自己回帰分析、及びその他の簡単な予測法
将来予測の技術は、統計分析の中でも急速な進歩が今も続いている分野である。本コースの目的は、一つの時系列データ（＝一変量）の今後の変化について予測することである。
そのための定番的方法であるボックス・ジェンキンズ流の予測法を到達目標にしている。ボックス・ジェンキンズ法とは、ARIMA（アリマ）モデルと呼ばれる予測式から誤差が最も小さくなる最善の計算式を決めるルーティン化された一連のやり方をさすのだが、ボックス・ジェンキンズ法が登場する以前にも、機械的ではあるが簡単で（マアマアの）精度を有する予測法が提案されていた。。
その中でも自己回帰分析や指数平滑法、ホルト・ウィンタース法が利用される機会は（今でも）多い。特に自己回帰分析はボックス・ジェンキンズ法の土台にもなっているので、今回のテーマにした。
次の第４回（M4）以降は、ボックス・ジェンキンズ流の予測法に入っていく。
消費支出月次データの分析
まず家計調査ベースの1世帯当たり消費金額を分析する。『家計調査』は、毎月、総務省統計局が実施している統計調査で、概ね8千世帯をサンプルにとって、その月の収入金額、品目別の消費支出額等を公表している。 今回使用するデータファイル“fies-201001-201409.csv”は、2010年1月から2014年9月までの消費支出合計を月次系列として保存している。
データファイルの読み込み、実額と前年（同月）比
データファイル中に年月も含めている。読み込んだあと、消費支出を時系列データに変換しておく。時系列をグラフに描く際は折れ線グラフがデフォールトで選ばれる。
data0 <- read.csv(file="fies-201001-201409.csv",header=T)
head(data0)## year month syouhi
## 1 2010 1 291918
## 2 2010 2 261163
## 3 2010 3 319991
## 4 2010 4 299996
## 5 2010 5 280714
## 6 2010 6 276494syouhi <- data0[,3]
syouhi <- ts(syouhi,freq=12,start=c(2010,1))
plot(syouhi)
図から明らかなように毎年12月には消費支出が増えるというハッキリした季節性がある。そして、季節性を除けばデータ全体に増加（減少）傾向や循環的な変動は認められず、期間を通して概ね横ばいである。
自己回帰分析の勘所は『前期に平均を超えたとすれば、今期はどうか？』という見方である – これは一次の自己回帰である。季節性が混じっていると、12月の消費が年間平均を超えたからといって、それは当たり前であり、1月の消費を予測する情報にはならない。
自己回帰で予測を行う大前提は、データの中の一定部分に着目すること、つまり定常化を済ませておくことである。定常化とは平均（＝高さ）と分散（＝バラツキ）を一定にしておく準備作業のことだが、今回の場合は、消費金額には季節性があり年間の水準が一定にそろわないので、前年比でみる。こういう手順になる。
高さとばらつきを概ね一定にしておいてから、変化のパターンを探るという発想は、現在に至る将来予測の技術の出発点になっている。
log.syouhi <- log(syouhi)
dlog.syouhi <- diff(log.syouhi,12)
nobi <- dlog.syouhi*100
dlog.syouhi## Jan Feb Mar Apr May
## 2011 -0.009385571 -0.001417744 -0.091881593 -0.025102786 -0.016359571
## 2012 -0.021202405 0.026718800 0.040093260 0.031588461 0.041674618
## 2013 0.020313324 0.000910526 0.039762854 0.008028674 -0.019447303
## 2014 0.027769513 -0.005569399 0.088560259 -0.007389696 -0.039569825
## Jun Jul Aug Sep Oct
## 2011 -0.039418639 -0.018496247 -0.039468493 -0.019645758 -0.006380052
## 2012 0.014947526 0.011534882 0.014182238 -0.012315817 -0.004797822
## 2013 -0.001453931 0.009845652 -0.004871374 0.051114804 0.022397329
## 2014 0.012441857 -0.020498925 -0.008899612 -0.019665404
## Nov Dec
## 2011 -0.038682105 0.003278962
## 2012 0.001257310 -0.007919597
## 2013 0.020871215 0.027098677
## 2014nobi## Jan Feb Mar Apr May Jun
## 2011 -0.9385571 -0.1417744 -9.1881593 -2.5102786 -1.6359571 -3.9418639
## 2012 -2.1202405 2.6718800 4.0093260 3.1588461 4.1674618 1.4947526
## 2013 2.0313324 0.0910526 3.9762854 0.8028674 -1.9447303 -0.1453931
## 2014 2.7769513 -0.5569399 8.8560259 -0.7389696 -3.9569825 1.2441857
## Jul Aug Sep Oct Nov Dec
## 2011 -1.8496247 -3.9468493 -1.9645758 -0.6380052 -3.8682105 0.3278962
## 2012 1.1534882 1.4182238 -1.2315817 -0.4797822 0.1257310 -0.7919597
## 2013 0.9845652 -0.4871374 5.1114804 2.2397329 2.0871215 2.7098677
## 2014 -2.0498925 -0.8899612 -1.9665404par(mfrow=c(2,1))
plot(dlog.syouhi,main="対数の前年差")
plot(nobi,main="前年比（パーセント")
データの対数をとってから対数の差を求める計算は、そのデータの増加率（減少率）を出す計算とほぼ一致することが知られている — 但し、10%を超える大きな増加率の場合は二つの計算法で違いが目立ってくる（参照）田中孝文『Rによる時系列分析入門』、40～41ページ。
定義通りに前年比増減率を求めると以下のようになる。
(nobi.teigi <- diff(syouhi,12)/lag(syouhi,-12)*100)## Jan Feb Mar Apr May
## 2011 -0.93416644 -0.14167397 -8.77868440 -2.47903305 -1.62264796
## 2012 -2.09792144 2.70789477 4.09078452 3.20926719 4.25551947
## 2013 2.05210438 0.09109406 4.05639792 0.80609906 -1.92594239
## 2014 2.81586798 -0.55539185 9.26000898 -0.73624590 -3.87971640
## Jun Jul Aug Sep Oct
## 2011 -3.86518333 -1.83262407 -3.86997590 -1.94540377 -0.63597430
## 2012 1.50597990 1.16016654 1.42832827 -1.22402874 -0.47863308
## 2013 -0.14528742 0.98942798 -0.48595282 5.24437112 2.26500327
## 2014 1.25195792 -2.02902502 -0.88601280 -1.94733017
## Nov Dec
## 2011 -3.79435070 0.32843434
## 2012 0.12581009 -0.78883199
## 2013 2.10905425 2.74691851
## 2014ts.plot(nobi,nobi.teigi,lty=c(1,2),col=c(1,2))
分子は前年同月との差、分母は前年同月の実績値である。なお、Rで入力する式をカッコで囲むと、結果を保存する場合であっても保存と同時に画面に表示される。
対数の前年差で求めた値（図の黒実線）と定義通りに求めた前年比（図の赤点線）を比べてみると、消費税率引き上げの直前の駆け込みがあった2014年3月でやや違いが目立ったものの、両者の値はほぼ一致していることがわかる。
今後、前年比増減率や前期比増減率を求めるときは、対数の差を求めるほうがコマンドとしては簡単なので近似法を常用することを了解されたい。
自己回帰分析と予測
今回のデータでは、前年比をとることによって消費の季節性が消え去り、ほぼ横ばいになった。そこで以下の作業をする。
- 最終月（2014年9月）の前月（8月）までのデータを用いて9月の前年比“nobi”を予測する。
- 一次の自己回帰、2次の自己回帰を推定してみる。
- 最小予測誤差モデルを探るためAIC（赤池情報量基準）に基づき最善の自己回帰式を求める。
- 最善の自己回帰式を用いて2014年9月の前年比を予測する。
- 実際の結果との予測誤差を求める。
以上の作業を行い、その後『予測誤差を最小化するモデルとはどんな意味だろうか』という問いかけを考えることにしよう。
nobi.sub <- window(nobi,end=c(2014,8))
(heikin.sub <- mean(nobi.sub))## [1] 0.1232193kekka1 <- ar(nobi.sub,order.max=1,aic=F)
kekka1##
## Call:
## ar(x = nobi.sub, aic = F, order.max = 1)
##
## Coefficients:
## 1
## 0.28
##
## Order selected 1 sigma^2 estimated as 8.517kekka2 <- ar(nobi.sub,order.max=2,aic=F)
kekka2##
## Call:
## ar(x = nobi.sub, aic = F, order.max = 2)
##
## Coefficients:
## 1 2
## 0.2099 0.2503
##
## Order selected 2 sigma^2 estimated as 8.178まずは一次の自己回帰と二次の自己回帰を推定した。行った計算は通常の回帰分析と（原理的には）同じなのであるが、求められた係数をみると回帰分析でみた切片と傾きのような形になっていない。
上の自己回帰分析から得られた関係は下の式である。
まず一次の自己回帰は \[ X_t-0.1232193 = 0.2800215\times (X_{t-1}-0.1232193) \]
また二次の自己回帰結果は下式のとおり（四捨五入している）。 \[ X_t-0.12 = 0.21\times (X_{t-1}-0.12)+0.25\times (X_{t-2}-0.12) \]
上に示した二次の自己回帰式は、前月と前々月の前年比から今月の前年比を予測している点で、GDP推計において消費の最終月を推定する方法と似ているところがある。ところが、上の式から\(X_t\)の予測式を出すと以下のようになる（注：回帰分析における通常の計算である「最小二乗法」とは異なる方法で解を求めているので、Rのlmコマンドの結果とはやや違う）。 \[ X_t = 0.21\times X_{t-1} + 0.25\times X_{t-2} + 0.065 \] この予測式は \[ X_t = 0.5 \times X_{t-1} + 0.5 \times X_{t-2} \] という実際に採用されていた予測法とは異なる。この点で、実際に採られていた予測法は、データの変動パターンと合致した最適な方法ではなかったことになる。
それでは得られた二次の自己回帰式から2014年9月を予測してみよう。それには以下のように“predict”コマンドを使う。予測期間数は1か月だからn.ahead=1と指定する。
predict(kekka2,n.ahead=1)## $pred
## Sep
## 2014 -0.6334694
##
## $se
## Sep
## 2014 2.8597872014年9月の実際の前年比は▲1.97%だった。事前予測値は▲0.63%だったから予測誤差は1.34%（＝▲0.63-(▲1.97)）、逆に言うと予測より1.34%だけ下振れした値が出た。
これに対して、実際の方法では次のように9月を予測していたはずである。
mean(tail(nobi.sub,2))## [1] -1.4699277月と8月の前年比の平均を事前予測値（＝▲1.47%）としていたから予測誤差は0.5%。実際には0.5%だけ下振れしたことになる。これは上の予測値より（結果的には）正確だったことになる。
あとで考えてもらう問いかけは次の疑問である。
【クイズ】
今月の前年比と前月・前々月の前年比に当てはまっている関係が二次の自己回帰式であったはずだ。ところが、これとは違う予測式から計算した結果のほうが正確だった。おかしくはないか。
自己回帰分析では、2か月だけ遡った実績のみを考慮する必然性はない。3か月前、4か月前の実績も考慮して予測値を求めてもよい。そんな時に便利であるのがモデル選択基準のAICである。この指標は、データに当てはまっている度合ではなく、将来に向けた予測誤差を最小化するという目的にかなうよう作成された指標である（計算式はパワーポイント参照）。
kekka.best <- ar(nobi.sub,order.max=12,aic=T)
kekka.best##
## Call:
## ar(x = nobi.sub, aic = T, order.max = 12)
##
## Coefficients:
## 1 2 3
## 0.1442 0.1952 0.2625
##
## Order selected 3 sigma^2 estimated as 7.805predict(kekka.best,n.ahead=1)## $pred
## Sep
## 2014 -0.1529755
##
## $se
## Sep
## 2014 2.793818最大で1年前の実績（＝order.max=12）まで説明変数に含めることにして、その中で予測誤差が最小になると示唆される次数をAICに基づいて探索（aic=T）しているのが一番上のコマンドである。
表示された結果を見ると、4か月より以前の実績は予測誤差を小さくする目的には寄与せず、直近3カ月の結果が大事であるという結果になっている — その中でも3か月前の実績が来月を予測するうえで大きなウェートをしめる。
この最善の自己回帰式を用いて2014年9月を予測すると▲0.15%となり、二次の自己回帰式より更に高い値となる（9月の予測誤差が大きくなる）— 3か月前の6月には前年比が1.2%のプラスだったから、この実績が織り込まれて予測値が上方に修正されたわけである。
【質問】
上のクイズをここで考えよ。 統計分析、機械学習などで一般に強調されている問題に『過学習』という概念がある。予測式をつくるときに用いたデータと、新たに追加されてくる新しいデータとは違うという点に関係する概念でもある。このことと上のクイズは非常に関係がある。
ホルト・ウィンタース法
上で紹介した（素朴な）自己回帰分析は、データを定常化（＝概ね横ばい）しておかないとうまく行かない。もちろん前年比を予測できれば、元の消費金額も予測できるのだが、作業コマンドで直接に消費金額の予測値を出してくれるわけではない。
ホルト・ウィンタース法を利用すると、データ数が少なく、ボックス・ジェンキンズの方法を使いづらい場合でも、季節変動を含んだ元のデータを直接的に予測することができる。（注）もう一つの機械的な方法である指数平滑法はここでは割愛する。
予測に用いるデータは2014年8月までに限定しておこう。
syouhi.sub <- window(syouhi,end=c(2014,8))
syouhi.hw <- HoltWinters(syouhi.sub,seasonal="add")
names(syouhi.hw)## [1] "fitted" "x" "alpha" "beta"
## [5] "gamma" "coefficients" "seasonal" "SSE"
## [9] "call"syouhi.hw.yosoku <- predict(syouhi.hw,n.ahead=12)
ts.plot(syouhi.sub,syouhi.hw$fitted[,1],type="l",lty=c(1,2),col=c(1,2))
ts.plot(syouhi.sub,syouhi.hw.yosoku,type="l",lty=c(1,2),col=c(1,2))
コマンド“HoltWinters”でseasonal=addとしているのは、データの基本成分であるTCSIが乗法モデルではなく加法モデルによって構成されていることを指示するものである — この点は、元のデータをグラフに描いて慎重に吟味しておかなければならない所である。
【クイズ】
ホルトウィンタース法では2014年9月の消費をどの程度と予測していたことになるか。その予測は自己回帰による予測と比べてどうだったか？
【実習課題】
- データファイル”passenger.csv”を読み込んで、将来12か月予測をしてください。
- データファイル”elec.csv”を読み込んで、将来12か月予測をしてください。
- データファイル”住宅着工戸数-196501-201508-元ファイル.xls”はe-Statからダウンロードした住宅着工戸数データの元ファイルである。これを利用して、本年9月以降6か月間の予測を総戸数について行いなさい。