2008年07月16日
2008年7月16日 ARIMAモデル
時系列データの解析と予測
時系列データを料理するために身につけるべきこと
① データの生成メカニズムを知る
② 予測する
データの発生メカニズム(DGP：data generating process)を解き明かすとは、すなわち何らかの方程式などのモデルに当てはめることと思えばよい。式が分かれば、予測も簡単である。しかし、統計的手法は万能ではないので、特別な前提をおく。
ここでは
●平均が一定（通常は標準化して０とおくと、やり易い）
●分散も一定（ずっと大きくなっていくとか、小さくなっていくとかしない）
この条件をみたすものを、とりあえず「定常時系列」とよぶ。
ちなみに、そうじゃないのを「非定常」という。
なので、実際のデータがこの条件に合わないとき、何らかの変換を加えて定常化すればよい。
１．データの生成機構を調べる方法
最初にトレンド(Ｔ)があるかをみる。
次は周期性(Ｃ)
そして季節性(Ｓ)などなど
→ プロットしたデータを見て特徴を記述する。
データが非定常なら、定常にする。
色々な方法がある。今回は差分をする。
２．予測する
今回は汎用性の高いARIMA/SARIMAをつかう。
ARFIMAやGARCH、カルマンフィルターなど時系列データの予測手法はデータの発生メカニズムに関連する。
ちなみに原典はBox=Jenkins(1976)『TIME SERIES ANALYSIS ,forecasting and control』です。
モデルの推計結果を使って、予測用の系列を作成し、描画出力する。
実際の予測作業はモデルの同定と診断に多くの労力を割き、時間もかかるが、今回は一通りの出力を得る練習と割り切って
ARIMAモデルの推計は
AR過程を１から３次程度
MA過程も１から３次程度
和分は１にする。
SARIMAの場合の
季節和分も１にする。
例えば
ARIMA(1,1,1)はAR1次、和分１階、MA1次となる。
ARIMA(2,1,1)はAR2次、和分１階、MA1次となる。
ARIMA(1,1,2)はAR1次、和分１階、MA2次となる。
ARIMA(2,1,2)はAR2次、和分１階、MA2次となる。
とりあえず、この４つ(SARIMAを除いて)を試して、いちばん良さげなものを採用し、予測系列を作成する。
※ モデルの診断に興味がある学生さんはAIC、SBICといったモデル選択のための統計量や、自己相関、編自己相関の詳しい見方を身につける必要がある。株で一儲け目論んでいる方は計量経済学の講義の時系列分析のパートを気合いをいれて勉強しましょう。
データはこれです。
例によって、ワークシート名[trianing]B1：B229を選択、コピーして
x<-read.table("clipboard",header=T)
でRに取り込む
このデータは1987年1月からなのでｔｓオブジェクトに変換して、表示
ty1<-ts(x,start=c(1987,1),freq=12)
plot(ty1)
TCSI分離で変動特性の目視
plot(decompose(ty1))
トレンドと周期成分をスペクトルで目視
spectrum(ty1)
ついでに
spectrum(ty1,method="ar")
最初に自己相関と編自己相関をみる
par(mfrow=c(2,1))
acf(ty1,lag.max=100)
pacf(ty1,lag.max=100)
par(mfrow=c(1,1))
トレンドを除去するか、１階差分でみる
dty1<-diff(ty1)
par(mfrow=c(3,1))
plot(dty1)
acf(dty1,lag.max=60)
pacf(dty1,lag.max=60)
par(mfrow=c(1,1))
とりあえず2次のAR、2次のMA、1次の差分でARIMAを推計する。
arima01<-arima(ty1,order=c(2,1,2))
arima01
この式は以下のように定式化できる
Yt=0.977Y(t-1)－0.1273Y(t-2)－0.1511m(t-1)+0.2913m(t-2)
Rのarima関数は残差系列しか出力されないので、原系列から残差を引いて、予測値系列とする
par(mfrow=c(2,1))
tres<-arima01$resid
ty1_hat<-ty1-tres
plot(ty1)
lines(ty1_hat,lty=2,col=2)
plot(tres,type="l",col=4)
abline(h=mean(tres))
par(mfrow=c(1,1)）
予測する。xtとした期間にちゅういすること
pred1<-predict(arima01,n.ahead=12)
xt<-c(2005,2007)
y1es<-pred1$pred
plot(ty1,xlim=xt)
lines(y1es,col=2)
以下はやってもやらなくてもよいけど、参考にしてくだされ
① 信頼限界付きかっこいい予測結果の出力
pred1<-predict(arima01,n.ahead=12)
xt<-c(2005,2007)
y1es<-pred1$pred
sig1<-pred1$se
tyl<-y1es-2*sig1
tyu<-y1es+2*sig1
yl<-c(min(ty1,tyl),max(ty1,tyu))
plot(ty1,xlim=xt,ylim=yl)
lines(y1es)
lines(tyl,lty=2,col=4)
lines(tyu,lty=2,col=4)
② 出力期間を変更してプロットする
pred1<-predict(arima01,n.ahead=24)
xt<-c(1987,2008)
y1es<-pred1$pred
sig1<-pred1$se
tyl<-y1es-2*sig1
tyu<-y1es+2*sig1
yl<-c(min(ty1,tyl),max(ty1,tyu))
plot(ty1,xlim=xt,ylim=yl)
lines(y1es)
lines(tyl,lty=2,col=4)
lines(tyu,lty=2,col=4)
③ 参考までに出力タイプ別でやってみた
xt<-c(2005,2007)
plot(ty1,type="p",xlim=xt,ylim=yl)
lines(y1es,type="p")
lines(tyl,lty=2,col=4)
lines(tyu,lty=2,col=4)
++++++++++++++++++++++++++++
季節階差付きSARIMAモデルの練習
++++++++++++++++++++++++++++
セルC1:C241を選択してコピー
x<-read.table("clipboard",header=T)
で読み込み、ｔｓオブジェクトに変換して、プロット
ty1<-ts(x,start=c(1987,1),freq=12)
plot(ty1)
par(mfrow=c(2,1))
acf(ty1,lag.max=100)
pacf(ty1,lag.max=100)
par(mfrow=c(1,1))
SARIMAの場合12か月毎の周期性があるためseasonal=list(order=c(1,1,1),period=12)というオプションを追加する。詳しくはｈｅｌｐ（ａｒｉｍａ）としてヘルプをみてみること。
sarima01<-arima(ty1,order=c(2,1,2),seasonal=list(order=c(1,1,1),period=12))
sarima01
推計結果の出力
par(mfrow=c(2,1))
tres<-sarima01$resid
ty1_hat<-ty1-tres
plot(ty1)
lines(ty1_hat,lty=2,col=2)
plot(tres,type="l",col=4)
abline(h=mean(tres))
par(mfrow=c(1,1))
12か月先の予測
pred1<-predict(sarima01,n.ahead=12)
xt<-c(2005,2009)
y1es<-pred1$pred
sig1<-pred1$se
tyl<-y1es-2*sig1
tyu<-y1es+2*sig1
yl<-c(min(ty1,tyl),max(ty1,tyu))
plot(ty1,xlim=xt,ylim=yl)
lines(y1es)
lines(tyl,lty=2,col=4)
lines(tyu,lty=2,col=4)
24か月先の予測(上との違いを確認してください)
pred1<-predict(sarima01,n.ahead=24)
xt<-c(1987,2009)
y1es<-pred1$pred
sig1<-pred1$se
tyl<-y1es-2*sig1
tyu<-y1es+2*sig1
yl<-c(min(ty1,tyl),max(ty1,tyu))
plot(ty1,xlim=xt,ylim=yl)
lines(y1es)
lines(tyl,lty=2,col=4)
lines(tyu,lty=2,col=4)
ここまで出来たらＩＩＰ(鉱工業生産指数)で練習してみます。
時系列データを料理するために身につけるべきこと
① データの生成メカニズムを知る
② 予測する
データの発生メカニズム(DGP：data generating process)を解き明かすとは、すなわち何らかの方程式などのモデルに当てはめることと思えばよい。式が分かれば、予測も簡単である。しかし、統計的手法は万能ではないので、特別な前提をおく。
ここでは
●平均が一定（通常は標準化して０とおくと、やり易い）
●分散も一定（ずっと大きくなっていくとか、小さくなっていくとかしない）
この条件をみたすものを、とりあえず「定常時系列」とよぶ。
ちなみに、そうじゃないのを「非定常」という。
なので、実際のデータがこの条件に合わないとき、何らかの変換を加えて定常化すればよい。
１．データの生成機構を調べる方法
最初にトレンド(Ｔ)があるかをみる。
次は周期性(Ｃ)
そして季節性(Ｓ)などなど
→ プロットしたデータを見て特徴を記述する。
データが非定常なら、定常にする。
色々な方法がある。今回は差分をする。
２．予測する
今回は汎用性の高いARIMA/SARIMAをつかう。
ARFIMAやGARCH、カルマンフィルターなど時系列データの予測手法はデータの発生メカニズムに関連する。
ちなみに原典はBox=Jenkins(1976)『TIME SERIES ANALYSIS ,forecasting and control』です。
モデルの推計結果を使って、予測用の系列を作成し、描画出力する。
実際の予測作業はモデルの同定と診断に多くの労力を割き、時間もかかるが、今回は一通りの出力を得る練習と割り切って
ARIMAモデルの推計は
AR過程を１から３次程度
MA過程も１から３次程度
和分は１にする。
SARIMAの場合の
季節和分も１にする。
例えば
ARIMA(1,1,1)はAR1次、和分１階、MA1次となる。
ARIMA(2,1,1)はAR2次、和分１階、MA1次となる。
ARIMA(1,1,2)はAR1次、和分１階、MA2次となる。
ARIMA(2,1,2)はAR2次、和分１階、MA2次となる。
とりあえず、この４つ(SARIMAを除いて)を試して、いちばん良さげなものを採用し、予測系列を作成する。
※ モデルの診断に興味がある学生さんはAIC、SBICといったモデル選択のための統計量や、自己相関、編自己相関の詳しい見方を身につける必要がある。株で一儲け目論んでいる方は計量経済学の講義の時系列分析のパートを気合いをいれて勉強しましょう。
データはこれです。
例によって、ワークシート名[trianing]B1：B229を選択、コピーして
x<-read.table("clipboard",header=T)
でRに取り込む
このデータは1987年1月からなのでｔｓオブジェクトに変換して、表示
ty1<-ts(x,start=c(1987,1),freq=12)
plot(ty1)
TCSI分離で変動特性の目視
plot(decompose(ty1))
トレンドと周期成分をスペクトルで目視
spectrum(ty1)
ついでに
spectrum(ty1,method="ar")
最初に自己相関と編自己相関をみる
par(mfrow=c(2,1))
acf(ty1,lag.max=100)
pacf(ty1,lag.max=100)
par(mfrow=c(1,1))
トレンドを除去するか、１階差分でみる
dty1<-diff(ty1)
par(mfrow=c(3,1))
plot(dty1)
acf(dty1,lag.max=60)
pacf(dty1,lag.max=60)
par(mfrow=c(1,1))
とりあえず2次のAR、2次のMA、1次の差分でARIMAを推計する。
arima01<-arima(ty1,order=c(2,1,2))
arima01
この式は以下のように定式化できる
Yt=0.977Y(t-1)－0.1273Y(t-2)－0.1511m(t-1)+0.2913m(t-2)
Rのarima関数は残差系列しか出力されないので、原系列から残差を引いて、予測値系列とする
par(mfrow=c(2,1))
tres<-arima01$resid
ty1_hat<-ty1-tres
plot(ty1)
lines(ty1_hat,lty=2,col=2)
plot(tres,type="l",col=4)
abline(h=mean(tres))
par(mfrow=c(1,1)）
予測する。xtとした期間にちゅういすること
pred1<-predict(arima01,n.ahead=12)
xt<-c(2005,2007)
y1es<-pred1$pred
plot(ty1,xlim=xt)
lines(y1es,col=2)
以下はやってもやらなくてもよいけど、参考にしてくだされ
① 信頼限界付きかっこいい予測結果の出力
pred1<-predict(arima01,n.ahead=12)
xt<-c(2005,2007)
y1es<-pred1$pred
sig1<-pred1$se
tyl<-y1es-2*sig1
tyu<-y1es+2*sig1
yl<-c(min(ty1,tyl),max(ty1,tyu))
plot(ty1,xlim=xt,ylim=yl)
lines(y1es)
lines(tyl,lty=2,col=4)
lines(tyu,lty=2,col=4)
② 出力期間を変更してプロットする
pred1<-predict(arima01,n.ahead=24)
xt<-c(1987,2008)
y1es<-pred1$pred
sig1<-pred1$se
tyl<-y1es-2*sig1
tyu<-y1es+2*sig1
yl<-c(min(ty1,tyl),max(ty1,tyu))
plot(ty1,xlim=xt,ylim=yl)
lines(y1es)
lines(tyl,lty=2,col=4)
lines(tyu,lty=2,col=4)
③ 参考までに出力タイプ別でやってみた
xt<-c(2005,2007)
plot(ty1,type="p",xlim=xt,ylim=yl)
lines(y1es,type="p")
lines(tyl,lty=2,col=4)
lines(tyu,lty=2,col=4)
++++++++++++++++++++++++++++
季節階差付きSARIMAモデルの練習
++++++++++++++++++++++++++++
セルC1:C241を選択してコピー
x<-read.table("clipboard",header=T)
で読み込み、ｔｓオブジェクトに変換して、プロット
ty1<-ts(x,start=c(1987,1),freq=12)
plot(ty1)
par(mfrow=c(2,1))
acf(ty1,lag.max=100)
pacf(ty1,lag.max=100)
par(mfrow=c(1,1))
SARIMAの場合12か月毎の周期性があるためseasonal=list(order=c(1,1,1),period=12)というオプションを追加する。詳しくはｈｅｌｐ（ａｒｉｍａ）としてヘルプをみてみること。
sarima01<-arima(ty1,order=c(2,1,2),seasonal=list(order=c(1,1,1),period=12))
sarima01
推計結果の出力
par(mfrow=c(2,1))
tres<-sarima01$resid
ty1_hat<-ty1-tres
plot(ty1)
lines(ty1_hat,lty=2,col=2)
plot(tres,type="l",col=4)
abline(h=mean(tres))
par(mfrow=c(1,1))
12か月先の予測
pred1<-predict(sarima01,n.ahead=12)
xt<-c(2005,2009)
y1es<-pred1$pred
sig1<-pred1$se
tyl<-y1es-2*sig1
tyu<-y1es+2*sig1
yl<-c(min(ty1,tyl),max(ty1,tyu))
plot(ty1,xlim=xt,ylim=yl)
lines(y1es)
lines(tyl,lty=2,col=4)
lines(tyu,lty=2,col=4)
24か月先の予測(上との違いを確認してください)
pred1<-predict(sarima01,n.ahead=24)
xt<-c(1987,2009)
y1es<-pred1$pred
sig1<-pred1$se
tyl<-y1es-2*sig1
tyu<-y1es+2*sig1
yl<-c(min(ty1,tyl),max(ty1,tyu))
plot(ty1,xlim=xt,ylim=yl)
lines(y1es)
lines(tyl,lty=2,col=4)
lines(tyu,lty=2,col=4)
ここまで出来たらＩＩＰ(鉱工業生産指数)で練習してみます。