ボックス・ジェンキンズ法（含む、説明変数）：他の要因を考慮する
まずパッケージ“forecast”をロードしておく。
library(forecast)データ例（１）：販売高と宣伝費
データファイル“sales-Advertising.csv”を読みこみ、“data0”という名前を付けて保存する。但し、このデータファイルにはヘッダーがついていない。第1列が商品の売上高、第2列が広告費である。単位は省略する。コマンドread.csvのカッコ内にheaderの有無を指定しないときはheader=TRUEになるので、ヘッダーがない場合は“header=FALSE”が必要である。
data0 <- read.csv(file="Sales-Advertising.csv",header=FALSE)
dim(data0)## [1] 36 2head(data0)## V1 V2
## 1 12.0 15
## 2 20.5 16
## 3 21.0 18
## 4 15.5 27
## 5 15.3 21
## 6 23.5 49ヘッダーがない場合は読んだファイルの第1列が“V1”、第2列が“V2”と自動的に名前がつく。これでは何のデータか分かりにくいので、第1列を“sale”、第2列を“promo”と名付ける。
colnames(data0) <- c("sale","promo")
head(data0)## sale promo
## 1 12.0 15
## 2 20.5 16
## 3 21.0 18
## 4 15.5 27
## 5 15.3 21
## 6 23.5 49データはまだ時系列になっていないのでtsコマンドで時系列データにしよう。
しかし、このデータファイルには — おそらく毎月の販売高と宣伝費を記録したものであろうが — 日付が含まれていない。そこでカレンダーの日付なしの、つまり時点番号のみの時系列にする。
sale <- ts(data0$sale)
promo <- ts(data0$promo)いま販売高を予測したいと考えている。そこで得られる予測モデルの精度を検証するため、「学習データ」と「テストデータ」に分けることにする。学習データは元のデータの時点1から時点30、テストデータは時点31から最終時点（＝36）までとする。
sale.train <- window(sale,end=30)
sale.test <- window(sale,start=31)
promo.train <- window(promo,end=30)
promo.test <- window(promo,start=31)単純なアリマ予測
まず学習データに基づいて販売高の予測式を推定する。 これはパッケージforecastのauto.arimaで簡単に行う。
kekka1 <- auto.arima(sale.train)
kekka1## Series: sale.train
## ARIMA(1,0,1) with non-zero mean
##
## Coefficients:
## ar1 ma1 intercept
## 0.4076 0.7093 23.3922
## s.e. 0.2057 0.1670 2.0704
##
## sigma^2 estimated as 18.26: log likelihood=-85.25
## AIC=178.51 AICc=180.11 BIC=184.11販売高の通期平均値(=23.4)との乖離に対して前期の実績及び前期の予測誤差を考慮した予測モデルARIMA(1,0,1)が最も当てはまっていることが分かる。
この結果に基づいて時点31から最終時点までを予測し、それを実績値と比較する。
yosoku1 <- forecast(kekka1,h=6)
plot(yosoku1)
points(sale.test,pch=20)
得られた予測値は全体として濃いグレーのゾーン（80%予測区間）には収まっているが、その変動パターンは全くフォローできていない — 6時点先にはいずれ下がるであろうという方向感だけは出ているが。
もしも販売高が宣伝費と相関があるならば、宣伝費の方は自社が将来期間にわたって予め決定することができるので、より正確な販売予測につながるのではないだろうか。これを確かめよう。
宣伝費を説明変数とするアリマ予測
そこで宣伝費を説明変数に含めた予測式を推定する。推定においては宣伝費についても学習データを用いる点に注意せよ。
kekka2 <- auto.arima(sale.train,xreg=promo.train)
kekka2## Series: sale.train
## ARIMA(0,1,0)
##
## Coefficients:
## promo.train
## 0.0926
## s.e. 0.0431
##
## sigma^2 estimated as 20.62: log likelihood=-84.52
## AIC=173.04 AICc=173.5 BIC=175.78結果を用いて予測する。予測期間の宣伝費はpromo.testであるのでこれをxregに指定している点を見よ。
yosoku2 <- forecast(kekka2,h=6,xreg=promo.test)
plot(yosoku2)
points(sale.test,pch=20)
今度は、変動パターンは（ある程度まで）フォローできているが水準が過大に過ぎるようである — 宣伝費を考慮して販売を予測すると変動パターンをフォローできるようになる、ということは宣伝と販売には関係があるということである。宣伝費を説明変数に含めておくことの効果はみてとれるのではないだろうか。おそらく宣伝費以外の要因、たとえば価格設定、景気等々から下振れしたものと思われる。
【補足】
説明変数“promo”を含めた予測モデルがARIMA(0,1,0)というのはどういう意味だろうか？
説明変数xregを含める時の予測モデルは以下の形を考えている。
\[ Y_t = a + b X_t + n_t \\ n_t \sim Arima(p,d,q) \]
要するに、回帰分析の残差についてARIMAモデルを当てはめる。これが基本になっている — 複数のデータを分析するときの基本型は（実は）これ以外にもある。この辺はパッケージ“forecast”を開発したHyndmanの考え方が反映されているところでもある。
上の式の前期差をとると
\[ \Delta Y_t = b \Delta X_t + \Delta n_t \]
こういう式になる。
上でxreg=promo.trainと指定してauto.arimaで推定したが、得られた推定結果を解釈すれば『1回差をとると自己相関がなくなるARIMA(0,1,0)』という意味になる。つまり前期差をとった時のモデルの残差\(\Delta n_t\)は完全ランダムでARIMA(0,0,0)が当てはまる。差をとる前の元のモデルでいえば、\(n_t\)にARIMA(0,1,0)を当てはめるということになる。
実際、差をとってから切片抜きの回帰分析を行うと、上と同じ結果が出てくる。
lm(diff(sale.train)~diff(promo.train)-1)##
## Call:
## lm(formula = diff(sale.train) ~ diff(promo.train) - 1)
##
## Coefficients:
## diff(promo.train)
## 0.09258実習
以下の計算を行え。
q1.kekka <- Arima(sale.train,order=c(1,1,1),xreg=promo.train)
q1.yosoku <- forecast(q1.kekka,h=6,xreg=promo.test)
plot(q1.yosoku)
points(sale.test,pch=20)データ例（２）：GDPと消費需要
国内市場の消費とGDPには関係がある。但し、今期の消費は今期のGDPによって影響されるという単純な関係ではない可能性がつよい。
data1 <- read.csv(file="real-gdp-cons.csv")
dim(data1)## [1] 185 4head(data1)## year quarter gdp cons
## 1 1955 1 10667.3 7178.7
## 2 1955 2 10685.2 7065.6
## 3 1955 3 11258.9 7404.5
## 4 1955 4 14463.6 9047.9
## 5 1956 1 11531.6 7831.9
## 6 1956 2 11769.6 7884.2時系列データにしておく。
gdp <- ts(data1$gdp,freq=4,start=c(1955,1))
cons <- ts(data1$cons,freq=4,start=c(1955,1))
ts.plot(gdp,cons,lty=c(1,2))
図から明らかにようにデータには季節成分が含まれている。その季節成分はGDPや消費の水準上昇に伴って次第に拡大しているので乗法モデルを適用するべきである。そのため元のデータの対数をとる。
log.gdp <- log(gdp)
log.cons <- log(cons)
plot(log.gdp)
plot(log.cons)
ここでの目的は消費需要の予測である。
GDPとの関係を考慮する方がより精度の高い予測になるかどうかを見ることにしよう。
今回も最後の1年間をテストデータとして留保しよう。
gdp.train <- window(log.gdp,end=c(2000,1))
gdp.test <- window(log.gdp,start=c(2000,2))
cons.train <- window(log.cons,end=c(2000,1))
cons.test <- window(log.cons,start=c(2000,2))まずGDPとの関係を考慮しない予測をする。
kekka.cons1 <- auto.arima(cons.train,D=1)
yosoku.cons1 <- forecast(kekka.cons1,h=4)
ts.plot(window(cons.train,start=c(1990,1)),yosoku.cons1$mean,lty=c(1,2))
points(cons.test,pch=20)
ほぼ正確に将来1年間（2000年第2期から2001年第1期）の消費需要を予測していることが分かる。この間の予測誤差の大きさを平均二乗誤差で計算するには以下のようにすればよい。
mean((cons.test-yosoku.cons1$mean)^2)## [1] 5.902755e-05次にGDPを説明変数に含めた予測式を求めることにしよう。
kekka.cons2 <- auto.arima(cons.train,D=1,xreg=gdp.train)
yosoku.cons2 <- forecast(kekka.cons2,xreg=gdp.test,h=4)
ts.plot(window(cons.train,start=c(1990,1)),yosoku.cons2$mean,lty=c(1,2))
points(cons.test,pch=20)
この時の予測誤差を評価しておくと次のようになる。
mean((cons.test-yosoku.cons2$mean)^2)## [1] 2.185019e-05消費の予測にGDPを考慮すると、考慮しない場合よりも平均二乗誤差が半分に減った。
そもそもGDPとの関係を考慮せずに予測した時の誤差が小さいこともあり、また対数にしていることも手伝って数値的には小さくなっているが、有力な影響要因があればそれを考慮するほうが予測精度の向上につながることが分かる。
ただし、将来4期間のGDPの値があらかじめ分かっていると仮定するのもは現実的ではない。そこで、説明変数に使っているGDPを単純なARIMAモデルで予測し、その予測値を考慮して消費を予測してみよう。
kekka.gdp <- auto.arima(gdp.train,D=1)
yosoku.gdp <- forecast(kekka.gdp,h=4)
yosoku.cons3 <- forecast(kekka.cons2,xreg=yosoku.gdp$mean,h=4)
mean((cons.test-yosoku.cons3$mean)^2)## [1] 3.404784e-05予測期間内のGDP実績値を使った際の消費予測に比べると誤差が拡大するのは当然である。しかし、単純なARIMA予測によった場合に比べて誤差は小さい点に注意せよ。
将来4期間のGDP実績をそのまま使った時に比べれば、GDP予測値を使う分、消費の予測力は落ちていることが分かる。しかし、それでも消費データだけから消費を将来予測するときに比べれば予測精度は高くなっていることが確認できる。
月並みな言い方だが、予測の際に考慮する材料を増やせば、それだけ予測力は高まる可能性がある。ボックスジェンキンズ法によるARIMAモデルは、単純にして、しかもかなり高い予測精度を提供する計算技術だが、これに甘んじない方法もあるわけであり、さらに言えば予測精度の向上は最近時点において最も活発な研究分野であることも付記しておきたい。
一層の勉強の余地は多々残されている。
実習
今回の消費予測ではデータに季節成分があるにもかかわらず、その処理はすべてauto.arimaに任せてしまっている。少なくともauto.arimaではどのような予測モデルが推定され採用されたのか、推定結果、診断図をみておくべきである。上の分析に関連してこの作業を行え。
【演習課題】
データファイル“HousingStart-196501-201508-ByCategory-ExclKyuuyo.csv”を用いて、最終月から将来1年間の「持家着工」を予測したい。作業を行え。分析対象期間は、各自、適切に選ぶこと（完璧を期す必要はない）。
持家着工の時系列データだけに基づいて予測せよ。但し、事前にデータ全体を学習データとテストデータに分けて予測力を検証しておけ。
持家以外の着工データを考慮することで持家の予測精度は高まるのか？この点を分析せよ。