景気予測と底打ち判断
西山 茂
2021年5月26日
5月公表の景気動向指数に基づいて予測分析を行ってみる。特に2020年3月時点から12か月後の2021年3月における値を予測する。月を追って最新値を情報に追加して、翌年3月値を予測するとき、時間経過とともにどのように予測誤差が縮小するかを視覚化する。
データ入力と要約
data0 <- read.csv(file="C:/Users/nisiy/Documents/1_Data/1_DICI/ci-20210512.csv")データファイルは標準的なCSVファイルで、第１列と第２列が年(year)と月(month)、第３列、第４列、第５列がそれぞれ先行指数（le）、一致指数（co）、遅行指数（lg）となっている。第１行は各列のデータ名である。
dim(data0)## [1] 435 5head(data0)tail(data0)各指数を個別に取り出したうえで日付をつけて時系列型（ts）にする。
le <- ts(data0[,3], freq=12, start=c(1985,1))
co <- ts(data0[,4], freq=12, start=c(1985,1))
lg <- ts(data0[,5], freq=12, start=c(1985,1))各指数の1985年1月以降の推移をグラフに描くと以下のようになる。
ts.plot(le,co,lg, lty=c(2,1,2), col=c(1,2,3))
legend(
"topleft",
legend=c("leading","coincident","lagging"),
lty=c(2,1,2),
col=c(1,2,3),
cex=0.8,
horiz=TRUE)
予測計算と精度評価
予測作業のために、パッケージ“forecast”をロードしておく — OBS授業『将来予測の技術』を履修した方たちは使い慣れていることと思います。
library(forecast)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zooちなみに、パッケージ“forecast”ではautoplotが使える。これはグラフ作成向けパッケージの定番である“ggplot2”を基礎にしているので、ggplot2をロードしておくとよい。
library(ggplot2)上では日付を与えて時系列型にしたが、下のように日付なしで時系列型にすることも可である。例えば、日次データや週次データの場合は、休日、祝祭日が混じり、あえて日付を与えない方がデータとしては使いやすい場合がある。
autoplot(ts(data0[,3:5]))
先行指数（le）の予測作業
モデル同定と推定
この資料では、先行指標“le”を予測対象にする。一致指数、遅行指数を予測する場合も、“le”を“co”、“lg”に変えるだけで全く同じである。
予測にあたっては、2020年3月までのデータをトレーニングデータ、2020年4月以降をテストデータとする。予測モデルはトレーニングデータに基づいて推定し、予測ではモデル推定では使わなかったテストデータを予測する。そして予測値とテストデータ実績値との誤差の大小によって予測精度を評価する。
le.trn <- window(le, end=c(2020,3))
le.tst <- window(le, start=c(2020,4))予測方法としてはボックス・ジェンキンズ法によるアプローチをとる — これまたOBS授業『将来予測の技術』を履修した方なら周知のことと思います。
まずトレーニングデータであるle.trnについて、自己相関図、偏自己相関図を吟味し、モデル同定を試みる。
acf(le.trn)
pacf(le.trn)
自己相関図は「高止まり型」、偏自己相関図のラグ1が概ね1である。明らかに非定常だと判断できるので、階差をとってから再描画する – 定常・非定常の判別には、ディッキー・フラーの単位根検定などが利用できるが、（偏）自己相関図を目視することで容易に判断できることが多いので、ここでは敢えて厳密な検定は行わない。
階差をとる場合、景気動向指数を作成するための基礎データは季節調整済みであるという建て前であるので、前年同期比や前年同期差などの季節階差はとらない。しかしながら、後でみるように、先行指数には（ないはずの）季節成分が残存していることが分かる。
acf(diff(le.trn))
pacf(diff(le.trn))
目視によってモデル同定を行うときの要点は、
自己相関図、偏自己相関図を「減衰型」、「カット型」のいずれかに分類する。
その際、カット型☓カット型という組み合わせはない。
減衰型☓カット型には自己回帰モデル、カット型☓減衰型には移動平均モデル、減衰型☓減衰型には自己回帰移動平均モデルが予測モデルとしては当てはまる。
モデル型を決めたあと、（偏）自己相関係数が有意に出ているスパイクの数を数え、モデルの次数を決める。自己回帰移動平均の混合型の場合は、先ずはそれぞれの次数を1から始めて、回帰診断図を吟味しながら改良する。
これを基本方針とするとよい。
得られた自己相関図を目視した判断としては
自己相関図は減衰型。偏自己相関図はカット型で有意なスパイク数は3本。モデルとしてはARIMA(3,1,0）を採る。
自己相関図はカット型で有意なスパイク数は左端のラグ0を除く4本、偏自己相関図は減衰型とする。モデルとしてはARIMA(0,1,4)を採る。
自己相関図、偏自己相関図とも減衰型とみる。モデルとしては単純なARIMA(1,1,1)とより複雑なARIMA(3,1,3)を候補として採る。
なお、景気動向指数算定の基礎データは季節調整されているので、モデル定義に季節部分は含めないことにする。
これらのモデルを推定したうえで、最小のAIC値を示すものを採用する。
mod1 <- Arima(le.trn,order=c(3,1,0))mod2 <- Arima(le.trn, order=c(0,1,4))mod3 <- Arima(le.trn, order=c(1,1,1))mod4 <- Arima(le.trn, order=c(3,1,3))summary(mod1)## Series: le.trn
## ARIMA(3,1,0)
##
## Coefficients:
## ar1 ar2 ar3
## 0.2607 0.1567 0.1709
## s.e. 0.0493 0.0507 0.0493
##
## sigma^2 estimated as 1.162: log likelihood=-629.12
## AIC=1266.24 AICc=1266.33 BIC=1282.42
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.001816579 1.072828 0.7817234 0.006183776 0.8574685 0.1424845
## ACF1
## Training set 0.01335762summary(mod2)## Series: le.trn
## ARIMA(0,1,4)
##
## Coefficients:
## ma1 ma2 ma3 ma4
## 0.2757 0.2496 0.2836 0.1073
## s.e. 0.0508 0.0498 0.0448 0.0523
##
## sigma^2 estimated as 1.159: log likelihood=-628.04
## AIC=1266.08 AICc=1266.23 BIC=1286.31
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.005656243 1.070076 0.790121 0.007148848 0.8656823 0.1440151
## ACF1
## Training set -0.0008348118summary(mod3)## Series: le.trn
## ARIMA(1,1,1)
##
## Coefficients:
## ar1 ma1
## 0.7745 -0.4735
## s.e. 0.0566 0.0748
##
## sigma^2 estimated as 1.181: log likelihood=-633.11
## AIC=1272.21 AICc=1272.27 BIC=1284.35
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.0023447 1.083105 0.7991381 0.006148584 0.8765713 0.1456587
## ACF1
## Training set -0.02785963summary(mod4)## Series: le.trn
## ARIMA(3,1,3)
##
## Coefficients:
## ar1 ar2 ar3 ma1 ma2 ma3
## 0.5652 -0.7606 0.6374 -0.2942 0.9274 -0.3412
## s.e. 0.0870 0.0428 0.0706 0.1046 0.0236 0.0979
##
## sigma^2 estimated as 1.123: log likelihood=-620.83
## AIC=1255.65 AICc=1255.92 BIC=1283.97
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.002646111 1.050731 0.793413 0.006199295 0.8698958 0.1446152
## ACF1
## Training set 0.001833419以上の結果から、目視による最適モデル選択としては、ARIMA(3,1,3)を採用する。
診断をすると、下図のように予測誤差はほぼホワイトノイズとなっていることが示唆され、問題はない。
tsdiag(mod4)
念のために、auto.arimaを利用して、最良のモデルを自動探索しておこう。
model.auto <- auto.arima(le.trn)summary(model.auto)## Series: le.trn
## ARIMA(1,1,4)(0,0,2)[12]
##
## Coefficients:
## ar1 ma1 ma2 ma3 ma4 sma1 sma2
## -0.9013 1.1653 0.4869 0.5063 0.3102 -0.1598 -0.1818
## s.e. 0.0768 0.0852 0.0773 0.0707 0.0440 0.0523 0.0496
##
## sigma^2 estimated as 1.106: log likelihood=-617.33
## AIC=1250.66 AICc=1251.01 BIC=1283.02
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.02135903 1.04168 0.7732326 0.0202143 0.8476153 0.1409369
## ACF1
## Training set 8.232076e-05景気動向指数算定の基礎データは季節調整済みであるのだが、自動探索の中でデータに残存する季節成分が検出されたことが分かる。実際、季節パートの係数は有意に推定されている。AIC値をみると、自動選択されたモデルの方が“mod4”よりも一層低くなっている。
かなり複雑なモデルが選択されており、オーバーフィッテイングの懸念があるものの、予測作業ではこの自動選択された結果も参考として使うことにする。
診断図にも問題はない。
tsdiag(model.auto)
12か月予測と予測精度の比較
fore <- forecast(mod4, h=12)plot(fore,xlim=c(2018.0, 2021.5))
points(le.tst)
2021年3月時点の予測経路は、4月以降、低下傾向が一貫して続くというものだった。
ところが、実際には、◯でプロットされているテストデータから分かるように、予想経路を大きく上回る急回復をとげたわけである。
この結果は、自動探索された“model.auto”によっても同じである。
fore.auto <- forecast(model.auto, h=12)
plot(fore.auto,xlim=c(2018.0, 2021.5))
points(le.tst)
自動探索から得られた“model.auto”による12か月予測では、比較的早期に先行指数は底を打ち、回復に向かうという予測になっている。この予測は、実際の動きに合致しており、予測モデルとしては目視による“mod4”を上回る性能をもっていることが分かるが、それでも急な回復をフォローできてはいない。
予測誤差の評価と比較
予測誤差の評価には“forecast”の“accuracy”を使う。
accuracy(mod4)## ME RMSE MAE MPE MAPE MASE
## Training set 0.002646111 1.050731 0.793413 0.006199295 0.8698958 0.1446152
## ACF1
## Training set 0.001833419accuracy(model.auto)## ME RMSE MAE MPE MAPE MASE
## Training set 0.02135903 1.04168 0.7732326 0.0202143 0.8476153 0.1409369
## ACF1
## Training set 8.232076e-05トレーニングデータに対する予測誤差の平均は、平均絶対誤差(MAE)、平均二乗誤差の平方根（Root Mean Squared Error: RMSE)のいずれをみても、自動探索で得られた予測モデルのほう誤差は小さく、その分予測精度が高いといえる。
しかし、トレーニングデータに対する予測誤差の小ささはオーバー・フィッテイングの危険を示唆していることもあるので、モデル推定には使用しなかったテストデータに対する予測誤差をみる必要がある。
accuracy(fore, le.tst)## ME RMSE MAE MPE MAPE MASE
## Training set 0.002646111 1.050731 0.793413 0.006199295 0.8698958 0.1446152
## Test set 11.284474521 14.537177 12.807002 11.493969573 13.4262382 2.3343288
## ACF1 Theil's U
## Training set 0.001833419 NA
## Test set 0.740553243 4.508357accuracy(fore.auto, le.tst)## ME RMSE MAE MPE MAPE MASE
## Training set 0.02135903 1.04168 0.7732326 0.0202143 0.8476153 0.1409369
## Test set 9.31625677 12.07668 10.7725045 9.4770847 11.3264628 1.9635014
## ACF1 Theil's U
## Training set 8.232076e-05 NA
## Test set 7.323769e-01 3.752208上に見るように、テストデータ期間（2020年4月～2021年3月）における予測誤差（＝実績値－予測値）をみると、モデル“mod4”では平均絶対誤差(MAE)が約12.8、モデル“model.auto”で約10.8となっている。RMSEでみても同様である。前に図示しているようにモデル“model.auto”の方が予測モデルとしては優れているが、それでもテスト期間中の予測誤差はいずれの予測モデルでも大きく、2020年4月以降のボトムアウト、ボトムアウト後の急速な回復は、2020年3月時点では予測不可能であったことが分かる。
トレーニングデータとテストデータとで予測誤差が極端に違う場合、オーバー・フィッテイングが示唆されるのであるが、自動探索で得られた複雑なモデルと、相関図の目視から得られた（比較的）単純なモデルと、いずれにおいてもボトムアウト後の回復を予測できていないことは共通している。
これは、線形の予測モデルの限界を示すもので、特に下方転換点、上方転換点、転換点前後の動きを既存データから予測することは、そもそも難しい課題であるということの現れだと解釈してもよい。
直近月情報を追加した逐次予測
2020年3月時点では、ボトムアウトとボトムアウト後の急速な回復を正確に予測できなかったとしても、それでは毎月新たな実績値が情報として追加される中で、2021年3月時点の景気予測がどのように修正されていったのかを、視える化してみよう。
トレーニングデータとテストデータの長さを確認しておく。
length(le.trn)## [1] 423length(le.tst)## [1] 122020年4月以降、指定した月数を推定期間に繰り入れ、その分だけ、テスト期間数を削って予測計算をする。最終月の予測値が2021年3月の予測値であるから、これを保存していく。
予測モデルとしては“mod4”と“model.auto”の二つを使うことにする。
まず、毎月新たな実績値をトレーニングデータに追加し、その分だけ予測期間数を減じて将来予測をするという、この反復計算を行う関数を定義しておく。
mkfore <- function(th){
mod <- Arima(subset(le, end=423+th), order=c(3,1,3))
f <- forecast(mod, h = 12-th)
return(f$mean[length(f$mean)])
}
sqfore <- function(iter){
kekka <- numeric()
for(i in 0:iter){
f <- mkfore(i)
kekka <- c(kekka,f)
}
return(kekka)
}2020年3月時点における21年3月予測値、4月時点における21年3月予測値、・・・2021年2月時点における3月予測値は、上の関数を使って、以下のように得られる。
(mar21 <- sqfore(11))## [1] 79.23860 69.23026 70.75584 88.56817 91.39827 91.15980 94.67811 97.50623
## [9] 99.77810 97.08162 97.45724 98.996992020年3月時点の当初予測値から1月分の情報が加わるごとに2021年3月予測値がどのように変化していったか。その推移を図にしよう。21年3月予測値を茶の〇で示し、情報が1月分追加されるごとに予測確度が上がる様子を色が濃くなることで表現している。
下の作図では、plot関数の中の色をRGB値で指定するためrgb関数を使っている。
plot(fore,xlim=c(2018.0, 2021.5))
points(le.tst)
for(i in 1:12){
points(2021+2/12,mar21[i],col = rgb(red = 0.5, green = 0.3, blue = 0.1, alpha = i/12),pch=16,cex=2)
}
month <- as.character(c(seq(3,12),1:2))
for(i in 1:12){
text(2021+2/12+0.25,mar21[i],month[i],cex=1)
}
上のグラフをみると幾つかのことが分かる。
2021年3月予測値は、4月、5月とも当初予測の95パーセント区間に含まれているものの、大きく下方に修正されており、今後１年間の「想定外の落ち込み」が心配されていた様子が見てとれる。ところが、早くも7月には実績値が当初予測（3月時点）の95パーセント区間（薄いグレーの範囲）の上方外側にはずれている。その後、急速な回復が一貫して続き、当初の予測区間に戻ることはなかった。
一方、21年3月予測値が当初予測の95パーセント区間の外に出たのは、ようやく9月時点になってからである。おそらく当初予測が可能性の範囲として頭に入っている状況では、当初の予測モデルが当てはまらくなったことを確信できるのは、2020年秋以降であったろうと推測できる。
他方、上のグラフの濃いグレーで示された範囲は80パーセント予測区間である。21年3月予測値がこの80パーセント区間の外に出たのは、2020年6月である。
図に示された様子から考えると、急速に経済状況が変化している場合、将来予測の範囲としては95パーセント区間をみるのは余りに保守的に過ぎ、むしろ80パーセント区間を基本として変化の方向を推測し、レジームが変わったことを早期に検知するほうが得策である。こう言えるかもしれない。６月実績値は当初予測値から大きく上振れしており、2021年3月時点の見通しには当初の見通しが当てはまらなくなった、こんな見方の方が、素直な判断だろう。
同じ作業を“model.auto”で行う。但し、最適モデルは新たに1月分の情報が加わるごとに、自動探索で得られるので、その都度、選択される予測モデルは異なっている可能性がある–というより、その都度、異なった予測モデルを採って21年3月予測値を計算するステップを反復計算の中で、いわば視えない化するわけである。
mkfore <- function(th){
mod <- auto.arima(subset(le, end=423+th))
f <- forecast(mod, h = 12-th)
return(f$mean[length(f$mean)])
}
sqfore <- function(iter){
kekka <- numeric()
for(i in 0:iter){
f <- mkfore(i)
kekka <- c(kekka,f)
}
return(kekka)
}(mar21auto <- sqfore(11))## [1] 83.59841 76.71637 76.63831 88.85926 93.39103 94.81019 99.25759
## [8] 99.50894 100.70001 99.18247 99.33764 100.11387得られた21年3月予測値の見方は前と同じである。
plot(fore.auto,xlim=c(2018.0, 2021.5))
points(le.tst)
for(i in 1:12){
points(2021+2/12,mar21auto[i],col = rgb(red = 0.5, green = 0.3, blue = 0.1, alpha = i/12),pch=16,cex=2)
}
month <- as.character(c(seq(3,12),1:2))
for(i in 1:12){
text(2021+2/12+0.25,mar21auto[i],month[i],cex=1)
}
右端に表示されている月が重なっているので、みづらいが、図から分かる点をまとめると以下のようになる：
テスト期間中の実績値が、7月に当初予測（3月時点）の95パーセント区間（薄いグレーの範囲）の外に出ている点は変わらない。但し、こちらの予測モデルでは8月実績値もそれほど当初予測値とは違いがなく、当初予測の範囲内と判断する可能性もある。逆に言えば、そもそも優れた予測モデルはより正確な予測ができるわけであり、それ故に当初の予測モデルが「もはや当てはまらない」と判断する時機が遅れる可能性もある。優れた予測モデルであるが故の落とし穴とみてもよいかもしれない。
21年3月予測値が当初予測の95パーセント区間の外に出たのは9月時点であり、これは前と共通している。
前のケースでは6月実績値が80パーセント予測区間から外れているが、上のグラフをみると6月も当初予測の80パーセント区間内に入っている。6月実績値はボトムアウトを示唆する重要な情報であるが、今回のケースでは「21年3月の当初予測を変更する必要はなし」と判断してしまう余地がある。
まとめ
予測モデルとして何を選ぶかにもよるが、いずれにしても、昨年3月時点に見通されていた景気予測は、早ければ6月に、遅くとも9月時点には、大きく修正されていたものと思われる。