ボックス・ジェンキンズ法（２）
パッケージ“forecast”を使う
前回の授業で使用したデータファイル“FiveSeries.RData”に含まれていた系列“data1”をとりあげて、今回は予測用パッケージ“forecast”をまず使ってみよう。
パッケージを使えるようにlibraryコマンドで読み込んでおく。
library(forecast)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric## Loading required package: timeDate## This is forecast 7.3データファイルを読み込む。
load(file="FiveSeries.RData")
ls()## [1] "data.ima" "data1" "data2" "data3" "data4" "data5"
## [7] "sdp" "varp"もう一度、data1の変化を図で確認して変化のパターンを（偏）自己相関図でみておく。
plot(data1)
par(mfrow=c(2,1))
acf(data1)
pacf(data1)
増加トレンドを無視できない。実際、自己相関図は高止まり型であり、偏自己相関図の左端は概ね1である。
そこで前期差をとって、より明瞭に変化の特徴をみる。
plot(diff(data1))
par(mfrow=c(2,1))
acf(diff(data1))
pacf(diff(data1))
自己相関図は減衰型、偏自己相関図はスパイク2本のカット型である。故に、data1の変化には2次の自己回帰型が当てはまると推測できる。
前期差を変数“d.data1”という名前で保存し、それを予測の対象にしよう。
d.data1 <- diff(data1)
kekka1 <- Arima(d.data1,order=c(2,0,0))
kekka1## Series: d.data1
## ARIMA(2,0,0) with non-zero mean
##
## Coefficients:
## ar1 ar2 intercept
## 0.7500 -0.3325 0.3105
## s.e. 0.0938 0.0937 0.1705
##
## sigma^2 estimated as 1.02: log likelihood=-141.68
## AIC=291.36 AICc=291.78 BIC=301.78yosoku1 <- forecast(kekka1,h=6)
plot(yosoku1)
図に見るように、前期差は低下するものの、data1の元の値でいえば、増加はし続けるという予測になっている。しかし、前期差を直接の予測対象にするのは、（前期差そのものに関心がある場合は別として）いかにも分かりづらい。
そこでコマンドを次のように変える。
kekka2 <- Arima(data1,order=c(2,1,0))
kekka2## Series: data1
## ARIMA(2,1,0)
##
## Coefficients:
## ar1 ar2
## 0.7689 -0.3109
## s.e. 0.0945 0.0943
##
## sigma^2 estimated as 1.041: log likelihood=-143.2
## AIC=292.4 AICc=292.65 BIC=300.21yosoku2 <- forecast(kekka2,h=6)
plot(yosoku2)
Arimaコマンドのカッコ内で対象となる変数をd.data1から元の値であるdata1に修正した。それに応じて、orderの2番目の値を1とし、階差を1回とってから2次の自己回帰を当てはめると指定した。この結果に基づいて予測計算をした。ところが、得られた予測はほぼ横ばいである。
（新規追加： 前期差を予測した際、予測値はプラスの値だった。これは元の値が増加し続けるという意味である。このことと元も値を予測対象にしたときの上の結果は異なっている。）
この理由は、階差をとる場合は前期差の平均をゼロと前提するという点にある。足元の前期差がプラスであっても（＝増加していても）、一定の平均（＝ゼロ）に戻っていくと予測されるので、元のデータは横ばい状態に収束していくと前提されてしまう。これが理由である。
しかし、“data1”に関しては全期間にわたって前期差の平均値がプラスであり、全体として増加してきている。
mean(diff(data1))## [1] 0.305121全期間にわたって前期差の平均値がプラス（マイナス）であるとき、元のデータには「プラスの（マイナスの）ドリフトがある」という言い方をする。「トレンド」という用語と紛らわしいが、トレンドはもう少し短期の意味合いで、増加局面と減少局面が後退している場合でも「〇〇から△△までの期間には増加トレンド」がある。という言い方をする。
ゼロではないドリフトを含める場合、Arimaコマンドにおいてinclude.driftをTRUEに指定する。
kekka3 <- Arima(data1,order=c(2,1,0),include.drift=TRUE)
kekka3## Series: data1
## ARIMA(2,1,0) with drift
##
## Coefficients:
## ar1 ar2 drift
## 0.7500 -0.3325 0.3105
## s.e. 0.0938 0.0937 0.1705
##
## sigma^2 estimated as 1.02: log likelihood=-141.68
## AIC=291.36 AICc=291.78 BIC=301.78yosoku3 <- forecast(kekka3,h=6)
plot(yosoku3)
予測“yosoku3”なら納得できるだろう。ちなみに“yosoku1”と“yosoku3”は計算対象が前期差か、元のデータかの違いがあるだけで、得られた結果は同じである。要するに、前期差を自己回帰から予測して、元のデータを逆算して予測してくれるのがyosoku3であり、この点がArimaコマンドの便利な所である。
さて、前回の授業では省略した診断ステップをしておく。本当は、予測計算の前に、予測モデルは使えるかどうかを検証する意味で行っておくべきである。
tsdiag(kekka3)
診断図は推定した予測モデルからどのような予測誤差が発生したかを図にしたものである。
中段は予測誤差についての自己相関図である。予測誤差に何かの規則性、つまり自己相関でとらえられるはずの変動パターンが残っているとすれば、理屈からいっておかしい。なぜなら、予測誤差に規則性が残っているなら、その規則性を利用して予測の精度を上げることができるからである。予測誤差をもたらす要素は不規則成分であるはずなので、予測誤差は本来的にはランダムのはずである。故に、診断図の中段ではどの自己相関も有意であってはならない — 但し、有意水準５％で判定しているので、20回に1回程度は青い線の外に出てもおかしくない。
次に、最下段をみる。これはカバン検定といわれる結果を図にしているもので、横軸の値で示されるラグ数まで予測誤差の全ての自己相関係数がゼロとあると見なしてよいかどうかを調べている。ゼロとは判断できない場合は、青い線で示された５％ラインよりの下にドットが打点される。今回の図をみると、すべて青線のうえに打点されている。したがって、予測誤差は完全にランダムと見てもよい。
最上段にあるのは予測誤差の大きさを図にしている — 但し標準偏差で割って標準値にしている。大きさが2を超える誤差が頻繁に発生していれば不規則要素は（いわゆる）正規分布には従ってはいない可能性が強い — ARIMA分析では不規則成分は正規分布に従って発生するという前提を置いて予測式を推定しているため、この正規分布の仮定が当てはまらないとなると、「非正規型ショックモデル」をとることになるので結構面倒になってくる。また、大きさが4を超えるような箇所があれば、その時点で何らかの突発的で規模の大きい攪乱的なショックがあったと推測される。突発的な不規則変動は過去の実績からは予測できないものである。
今回の“data1”については診断図に問題はない。
自動探索機能を活用する
ボックス・ジェンキンズ法という将来予測における作業ルーティンは、（本来は）自己相関図と偏自己相関図を検討することによってデータに当てはまる（と思われる）予測モデルを識別する、つまり「同定」ステップが勘所である。
しかし、二つの図から適切な予測モデルを識別するという「読図」作業は、相当の分析経験をもつまでは中々できるものではない。そのためボックス・ジェンキンズの方法が（特に経済経営分野において）広く利用されるようになったのはそれほど古いことではない。普及の決定的要因になったのはソフトウェアの急速な進化に他ならない。
現在は（偏）自己相関図を読みとって予測モデルを同定するステップはソフトウェアにまかせるのが（作業現場では）普通である。パッケージ“forecast”でも自動選択機能がある。
kekka.auto <- auto.arima(data1)
yosoku.auto <- forecast(kekka.auto,h=6)
plot(yosoku.auto)
最適な予測モデルの選択から推定、予測、予測図の作成まで — 診断ステップを省略しているものの — わずか3行のコマンドで済むのである。
図のタイトルからも分かるように、予測モデルはARIMA(5,1,0)、つまり前期差に対して5次の自己回帰モデルを当てはめている。また“with drift”とあるので、「データの前期差は平均がゼロではない」つまり増加傾向が認められる。ここまでソフトウェアが自動的に分析している。
将来予測に利用する有力なソフトでは最適モデルの自動探索機能が当たり前になってきている。
クイズ
具体的にどんな予測式が求められたのだろうか？式に書いてみよ。
季節性のある時系列（実質GDP原系列）
上で分析したデータは実際のデータではなく季節成分も含まれていなかった。
次に分析するのは日本の実質GDP原系列である。データはRData形式になっているので（既にRで作業しているなら）loadコマンド、もしくはファイルをダブルクリックしてRをたちあげる。
load(file="Gdp-2014Q3.RData")ファイルには個別データとして“gdp”（1994年第1期から2014年第3期）と“gdp.s”が含まれており、gdp.sはgdpの最終四半期を除いたものになっている。
ここではデータ“gdp.s”を分析して2014年第3期を予測し、その予測値をデータgdpに含まれている実績と比較することにしよう。
まずデータを図に描いてみる。
plot(gdp.s)
増加トレンドが明瞭であり、かつ季節変動がある。しかし、季節成分は時間の経過とともに拡大しているとハッキリとは言えず、加法モデルを適用してもよさそうである。それ故、乗法モデルで必要な対数化はここでは採らないことにする。
だとすればgdp.sに基づく将来1期予測は非常に簡単である。
gdp.auto <- auto.arima(gdp.s)
gdp.auto## Series: gdp.s
## ARIMA(1,1,0)(1,0,0)[4]
##
## Coefficients:
## ar1 sar1
## 0.2140 0.8977
## s.e. 0.1136 0.0446
##
## sigma^2 estimated as 4220053: log likelihood=-735.07
## AIC=1476.13 AICc=1476.44 BIC=1483.32というモデルが選ばれたことがわかる。
ARIMAの後に二つの括弧がついている。これは季節ARIMAモデルを当てはめたという意味である。
季節ARIMA(アリマ)モデルとは、これまでのARIMAモデルに季節成分を含めるための修正を加えたものである。
上のARIMA記号の最初の括弧はこれまでと同じ意味をもつ。それぞれに数値は以下の意味である。
これに対して、二番目の括弧は季節要素に対応する部分である。それぞれの数値の意味合いは以下のようである。
（注）より詳細はパワーポイントおよび授業ノートを参照されたい。
上の推定結果をみると、自動推定機能が選んだ予測モデルは以下のようである。まず前期差を1回とっている（前の括弧）、前期差について1期前の実績で自己回帰を行っている。更に、前年同期の前期差も含めている。
実際には、季節性のあるデータを分析するとき、前期差をとるよりも前年同期との差をみるのが普通である。このことをRに指示することにしよう。
gdp.auto <- auto.arima(gdp.s,D=1)
gdp.auto## Series: gdp.s
## ARIMA(2,0,0)(0,1,1)[4] with drift
##
## Coefficients:
## ar1 ar2 sma1 drift
## 1.0109 -0.2547 -0.5636 233.8636
## s.e. 0.1142 0.1196 0.1029 85.1468
##
## sigma^2 estimated as 2733469: log likelihood=-687.92
## AIC=1385.83 AICc=1386.66 BIC=1397.61yosoku.gdp.auto <- forecast(gdp.auto,h=4)
plot(yosoku.gdp.auto)
points(window(gdp,start=c(2014,3),end=c(2014,3)),pch=17,col="red")
2014年第3期以降の1年間の予測ラインをグラフにした後、実際の第3期実績を赤い三角で書き加えている。
これをみると、第2期までの実績から事前に予測されていた値より第3期の公表値は下回っていたことがみてとれる。実績が事前予測値を下回る、いわゆる「下ぶれ」である。この下ブレの度合が大きければ「ネガティブ・サプライズ」と呼ばれて、その後の景気予測に影響する。経済データが公表される時にはこんな議論があるわけだが、今回の2014年第3四半期公表値は、それほどの下振れではない。「GDPショック」という状況ではなかったというほうが正確であった。
クイズ
上では自動推定機能で簡単にすませたが、もし自分が利用しているソフトウェアに自動機能がなければ、まず定常化してから（偏）自己相関図を読みとり予測モデルを同定しなければならない。各自、やってみよ。
季節性がある（かもしれない）時系列（自動車販売台数）
次に分析するのは日本の自動車販売台数（四半期ベース）である。
data0 <- read.csv(file="Car-Japan-1993Q1-2015Q2.csv")
dim(data0)## [1] 90 3head(data0); tail(data0)## year quarter sale
## 1 1993 1 3101281
## 2 1993 2 2776345
## 3 1993 3 2738075
## 4 1993 4 2611844
## 5 1994 1 2724604
## 6 1994 2 2515763## year quarter sale
## 85 2014 1 2664129
## 86 2014 2 2402147
## 87 2014 3 2380856
## 88 2014 4 2327533
## 89 2015 1 2480197
## 90 2015 2 2170771データ期間は1993年第1期から2015年第2期までである。
まずは販売台数“sale”を時系列データにしておこう。日付情報はいま読んだファイルに含まれている。
car <- ts(data0$sale,freq=4,start=c(1993,1))
head(car)## [1] 3101281 2776345 2738075 2611844 2724604 2515763plot(car)
1990年代後半に減少傾向をたどった後、2000年代に入って回復したもののリーマン危機で急減し、更に東日本大震災でも大きく減少している。その後、販売台数は以前の高さに戻らず概ね横ばい傾向のまま推移している。
クイズ
販売台数という元データの特徴を調べよ。「定常化」のために差をとる必要はあるだろうか？
【演習課題】
自動車販売台数について必要な作業を続けて、以下の予測を行え。
a. 2013年第4期までのデータを学習データ、2014年第1期から2015年第2期までのデータをテストデータとする。予測モデルの作成は学習データに基づいて推定し、将来予測はテスト期間である2014年第1期から2015年第2期までについて行う。
b． 予測誤差（＝予測精度）を評価するため以下のように定義される平均二乗誤差を求めよ。\[ 平均二乗誤差＝(実績値－予測値)^2の平均 \] （注）リーマン危機や東日本大震災によるサプライチェーン寸断によって自動車販売が大きな影響を受けたことは明らかである。このような外部環境の激変を考慮して予測する方がより良い予測につながる可能性がある。この点については次回にとりあげる予定である。
航空旅客数データ“passenger.csv”はボックス・ジェンキンズが著した古典的教科書にも掲載されている有名なサンプルデータである。これを読み込み、最終月以降1年間を予測せよ。最後の1年間をあらかじめ除いたデータから予測モデルをつくり、予測の信頼性を確かめて見よ。
データファイル“wine-wool.RData”にはワインとウールの月次生産量データが保存されている。最終月以降、将来6か月について予測せよ。