数量や金額や価格の推移の記録する時系列データは観察記録や社会統計で広く見られ、それしか無い場合も多いです。しかし、時系列データは統計解析をかけるには厄介な性質を持つことが多々あります。
実験データやアンケート調査の分析でよく使われる統計手法は通用しないことが多いです。非定常なデータは、定常なデータに変換する必要があります。定常であっても、往々にして自己回帰や系列相関といった性質があり、それらに応じた統計モデルが必要です。
厄介と言っても、これまでの統計学者の努力で時系列データの分析手法は整備され、教科書的なアプローチがあり、ソフトウェアの対応も進んでいます。Rにも時系列データを処理するためのビルトイン関数や様々なパッケージがあります。
1
tsオブジェクト
使い勝手が微妙で、季節調整をしないときは使わない方が良い気もしますが、Rには時系列データ用のts型があります。
ベクトルもしくは行列によるデータ系列に、開始時点（start）と周期（frequency）をつけ加えたオブジェクトで、関数tsで作成することができます。
YEAR MONTH JTAR UR
1 1998 1 0.63 3.6
2 1998 2 0.61 3.7
3 1998 3 0.57 4.1 Jan Feb Mar Apr May Jun Jul
1998 0.63 0.61 0.57 0.56 0.54 0.52 0.51
Aug Sep Oct Nov Dec
1998 0.50 0.49 0.48 0.47 0.47周期12は月次データとして解釈されます。周期が2以上の場合、startは何周期目かと周期内の位置を指定します。年と月を指定しました。
tsオブジェクトに設定したstartとfrequencyは同名のメソッドで確認可能です。
[1] 1998 1[1] 12tsオブジェクトの系列を年月など周期と季節を添字に取り出す時は、windowメソッドを持ちいます。
Jan Feb Mar Apr May Jun Jul
2010
2011 0.60 0.62 0.62 0.62 0.61 0.62 0.64
2012 0.74 0.75 0.77 0.78 0.79 0.80 0.81
Aug Sep Oct Nov Dec
2010 0.55 0.56 0.58 0.59
2011 0.65 0.67 0.69 0.71 0.72
2012 0.82 位置を指定するときは添字を使えますが、tsオブジェクトではなく、ベクトル扱いになります。
[1] 0.48 0.47行列をtsオブジェクトにすることもできます。
JTAR UR
Jan 1998 0.63 3.6
Feb 1998 0.61 3.7
Mar 1998 0.57 4.1四半期ごとに平均をとって、周期4のtsオブジェクトを作ってみましょう。
df_q <- aggregate(cbind(JTAR, UR) ~ I(rep(1:4, nrow(df01)/4)) + YEAR, mean, na.rm = TRUE, data = df01)
colnames(df_q)[1] <- "Q"
ts_q <- with(df_q, ts(cbind(JTAR, UR), start = c(YEAR[1], Q[1]), frequency = 4))モデル式内のI(...)は、内側のRのコードを実行し、戻り値を説明変数として扱うための関数です。I(rep(1:4, nrow(df01)/4))は、1, 2, 3, 4の繰り返しをつくって説明変数とします。
周期4は四半期データとして解釈されます。
JTAR UR
2010 Q1 0.5000000 5.066667
2010 Q2 0.5100000 5.066667
2010 Q3 0.5300000 5.033333
2010 Q4 0.5400000 5.033333
2011 Q1 0.6266667 4.566667
2011 Q2 0.6433333 4.566667
2011 Q3 0.6566667 4.633333
2011 Q4 0.6633333 4.566667周期は1でも53でも問題なく計算できますが、printによる表示は4と12のときのみカスタマイズされたものが使われます。
1.1 プロット
見栄えがよくないのでggplot2を使うことになる気もしますが、plot(JTAR)とすると散布図ではなくて折れ線グラフが出ます。系列が複数の場合は、
と複数の折れ線グラフを並べて描画してくれます。
tsオブジェクトをプロットしたときは、textなど低水準グラフィックス関数のx座標が年などstartで指定した基準になることに注意してください。周期は小数点以下の数字になります。
1.2 季節調節
季節調整をかけるのも容易です。
ts_y_s <- matrix(NA, nrow(ts_y), ncol(ts_y))
for(j in 1:2) ts_y_s[, j] <- apply(stl(ts_y[, j],
s.window = frequency(ts_y))$time.series[, c("trend", "remainder")],
1, sum)
plot(ts(ts_y_s, frequency = frequency(ts_y),
start = start(ts_y),
names = colnames(ts_y)), main = "季節調整値")
上のコードは、行列の列ごとに季節調整をかけて、戻り値の中のトレンドと残差を合計して季節変動分を捨てているだけです。
詳しくはRでお手軽季節調整を参照してください。
2 単変量時系列モデル
単一系列しかない時系列データで将来予測を行いたいときに、よく使われるのが自己回帰和分移動平均（ARIMA）モデルです。
ARIMAは、自己回帰（AR）モデルと移動平均（MA）モデルを統合した自己回帰移動平均（ARMA）モデルに、時系列に0以上の階差をとった変数のモデルをデータ生成過程（DGP）と捉えるものです。
2.1 ARモデル
前期までの値が、今期の値に影響するモデルです。ラグの数を括弧でくくりAR(1)のように表記することが多いです。
\[\begin{align} y_{t} = \alpha_0 + \sum^L_{\ell=1} \alpha_i y_{t - \ell} + \epsilon_t \end{align}\]
\(y\)が時系列データで、\(t\)は時点を表す添字、\(L\)はラグの数で、\(\epsilon\)が撹乱項です。
推定はOLSでもできますし、自己相関係数の方程式（Yule-Walker方程式）をLevinson-Durbinアルゴリズムで解く高速な方法もあります。欠損値に対応させるために、状態空間モデルの形式にして最尤推定する方法もあります。
2.2 MAモデル
\(t\)期の値が、\(t\)期の撹乱項だけではなく、\(t\)期以前の撹乱項の値にも影響されるモデルです。
ラグの数を括弧でくくりMA(2)のように表記することが多いです。
\[ y_t = \mu + \sum^L_{\ell=0} \theta_\ell \epsilon_{t - \ell} \]
\(\mu\)は平均です。
数学的にはラグ無限のMAモデルと定常なARモデルは相互に変換可能であるそうですが、実用上は別のものとして扱われます。
推定は、最尤法で行います。尤度は、パラメーター\(mu\)と\(\theta_{\ell}\)が与えられたとして、時点\(t=1\)から昇順に残差（誤差）を計算していき、算出します。
2.3 ARMAモデル
ARモデルとMAモデルを組み合わせたモデルです。よくARモデルのラグとMAモデルのラグを括弧内に入れてARMA(1, 2)のように表記されます。一番目の値がARモデルの、二番目の値がMAモデルのラグになります。
\[ y_{t} = \alpha_0 + \sum^{L_{AR}}_{\ell=1} \alpha_i y_{t - \ell} + \sum^{L_{MA}}_{\ell=0} \theta_\ell \epsilon_{t - \ell} \]
2.4 ARIMAモデル
時系列の階差をとった変数が、ARMAモデルに従うとしたモデルです1。ARのラグが2で、階差が1階、MAのラグが3のときは、ARIMA(2, 1, 3)のように表記されます。ARIMA(2, 0, 3)であればARMA(2, 3)と同じで、ARIMA(2, 0, 0)であればAR(2)と同じに、ARIMA(0, 0, 3)であればMA(3)と同等になります。
2.5 SARIMAモデル
SARIMAは、時系列を周期成分とそれ以外に分けて、それぞれARIMAモデルに従うと考えて推定する手法です。非時系列成分のARIMAのラグと階差、時系列成分のARIMAのラグと階差、時系列の周期をARIMA(1,0,1)(1,1,2)[12]と言うように表記されます。
3 時系列データの定常性
時系列データの分析では、定常過程であることが要求されます。非定常であると、多変量解析では見せかけの相関を拾うことが多く、単変量解析でも最尤法が収束しなかったり、推定結果が不安定になることが多いです2。
感覚的には、正の方向にも負の方向にもどんどん大きくならないのが定常です。統計学では数理的に扱えるように、同じ長さの任意の二つの期間の同時分布が等しい時系列データのことだと定義されています。\(\{X_{t}, X_{t+1}, \cdots, X_{t+l}\}\)と\(\{X_{t + s}, X_{t + s + 1}, \cdots, X_{t + s + l}\}\)の同時分布が同じと言うことです。強定常と呼ばれます。ただし、この強定常の定義は扱いづらいので、どの時点でも期待値と分散が同じで、同じ間隔の任意の2時点の自己共分散が同じである弱定常であるかどうかを調べることが多いです。強定常であれば必ず弱定常ですが、逆は成立しません。
具体的な例として、ARモデルが定常過程になる条件を整理しましょう。まず、ARモデルを、行列を用いて書きなおします3。
\[\begin{align} \begin{pmatrix} y_{t} \\ y_{t-1} \\ y_{t-2} \\ \vdots \\ y_{t-L+1} \end{pmatrix} = \begin{pmatrix} \alpha_0 \\ 0 \\ 0 \\ \vdots \\ 0 \end{pmatrix} + \begin{pmatrix} \alpha_1 & \alpha_1 & \cdots & \alpha_L \\ 1 & 0 & \cdots & 0 \\ 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 1 \end{pmatrix}\begin{pmatrix} y_{t-1} \\ y_{t-2} \\ y_{t-3} \\ \vdots \\ y_{t-L} \end{pmatrix} + \begin{pmatrix} \epsilon_t \\ 0 \\ 0 \\ \vdots \\ 0 \end{pmatrix} \end{align}\]
右辺第2項の正方行列を\(A\)、\(\mathbf{y}_t = {}^t\!(y_t, y_{t-1}, \cdots , y_{t-L})\)、\(\mathbf{c} = {}^t\!(\alpha_0, 0, \cdots , 0)\)、\(\mathbf{\varepsilon}_t = {}^t\!(\epsilon_t, 0, \cdots , 0)\)と置くと、
\[ \mathbf{y}_t = \mathbf{c} + A \mathbf{y}_{t-1} + \mathbf{\varepsilon}_t \]
と簡潔に書けます。
\(\mathbf{y}_t\)が\(t \rightarrow \infty\)で発散せずに、ARモデルが定常になる条件は、\(A\)の固有値の絶対値がすべて\(1\)未満になります。初期値\(\mathbf{y}_0\)、切片項\(\mathbf{c}\)、撹乱項\(\mathbf{\varepsilon}_t\)の影響が有限の値に収束しることを順番に確認しましょう。
3.1 初期値\(\mathbf{y}_0\)の影響の収束
まず、\(\mathbf{c}\)と\(\mathbf{\varepsilon}_t\)が\(0\)で、\(\mathbf{y}_t = A \mathbf{y}_{t-1}\)のケースを考えます。
\[ \mathbf{y}_t = A \mathbf{y}_{t-1} = A A \mathbf{y}_{t-2} = \cdots = A^t \mathbf{y}_0 \]
と書けます。\(A\)は正則なので、ある正則行列\(P\)と対角成分に固有値をとる対角行列\(\Lambda\)で、\(A=P \Lambda P^{-1}\)とかけるので
\[\begin{align} \mathbf{y}_t = (P \Lambda P^{-1})^t \mathbf{y}_0 = P \Lambda^t P^{-1} \mathbf{y}_0 \end{align}\]
となり、すべての固有値の絶対値が\(1\)未満であれば、\(t \rightarrow \infty\)で\(0\)に収束します。
3.2 切片項\(\mathbf{c}\)の影響の収束
\(\mathbf{c}\)が非ゼロの場合は、
\[\begin{align} \mathbf{y}_t = \mathbf{c} + A \mathbf{y}_{t-1} = \bigg( \sum^t_{k=0} A^k \bigg) \mathbf{c} + A^t \mathbf{y}_0 \end{align}\]
ここで\(I\)を単位行列として、\(A\)のすべての固有値が\(1\)未満であれば、\((I - A)(I + A + A^2 + A^3 + \cdots) = I\)から\(\sum^t_{k=0} A^k = (I - A)^{-1}\)になることから、\(t \rightarrow \infty\)で\(\mathbf{y}_t\)は有限の値\((I - A)^{-1} \mathbf{c}\)に収束します。
3.3 撹乱項\(\mathbf{\varepsilon}_t\)の影響の収束
撹乱項の影響が収束するかも検討しましょう。撹乱項の累積的な影響は\(\sum^{t-1}_{k=0} A^k \mathbf{\varepsilon}_{t-k}\)と書けます。\(E\big[\mathbf{\varepsilon}_t\big]={}^t\!(0, 0, \cdots , 0)\)なので、期待値はゼロです。分散共分散行列は\(\sum^{t-1}_{k=0} A^k E\big[ \mathbf{\varepsilon}_{t-k} {}^t\! \mathbf{\varepsilon}_{t-k} \big] ({}^t\!A)^k\)となります。\(E\big[ \mathbf{\varepsilon}_{t-k} {}^t\! \mathbf{\varepsilon}_{t-k} \big]\)は撹乱項の累積効果を考えない分散で、1行1列以外はゼロの定数になります。\(A\)のすべての固有値は\(1\)未満なので、\(t \rightarrow \infty\)で累積効果の分散は有限の値に収束します。\(A^k\ {}^t\!A^k < A^k\)なので、級数の内側が収束する列に上から抑えられていることに注意してください。
3.4 特性方程式
教科書的には\(A\)の固有方程式の固有値\(\lambda\)を\(x=1/\lambda\)として\(x\)で置き換えた特性方程式の解（根）の絶対値が、収束条件として示される事が多いです。
\[ 1 - \alpha_1 x - \alpha_2 x^2 - \cdots - \alpha_L x^L = 0 \]
\(|\lambda|<0\)が定常条件でしたから、\(|x|>1\)のときにARモデルは定常となります。
解が\(1\)の場合を単位根と呼びます。単位根がある場合は非定常になります。\(1\)より大きい解がある場合は爆発的と呼び、やはり非定常です。
3.5 単位根検定
時系列モデルが分かっていれば、数理的に定常か非定常かを識別することができますが、実データではそうは行きません。単位根があるかどうかを検定して調べます。複数の手法が提案されており決定版はないようですが、拡張ディッキーフラー（ADF）検定やKPSS検定が使われる事が多いようです。複数の検定を行い、すべてが単位根を棄却したときのみ単位根がないと見做す保守的な運用も見かけます。
4 単変量時系列モデルの推定
ラグと階差の次数が定まっている場合は、ビルトイン関数のarimaを使っても簡単です。
Call:
arima(x = ts_y[, "JTAR"], order = c(1, 1, 1))
Coefficients:
ar1 ma1
0.8833 -0.2641
s.e. 0.0310 0.0622
sigma^2 estimated as 0.0001238: log likelihood = 1031.11, aic = -2056.234.1 モデル（ラグと階差）の自動決定
ARIMAでも3、SARIMAだと6のパラメーターを決めた上で推定するわけですが、事前に定まっている事は少ないでしょう。階差は単位根検定で単位根が棄却できる次数にするのが通例で、パラメーターは統計量規準で求める事が多いです。
機械的に選ぶわけですが、ビルトイン関数のarimaを使うと推定を繰り返す事になるので煩雑です。便利なパッケージを使いましょう。例えば、forecast::auto.arimaを使うと簡単にモデル選択ができます。
library(forecast)
# JTAR系列とUR系列を順次推定
r_arima <- lapply(1:ncol(ts_y), \(j){
# 情報量規準はA.I.C.
# サーチ方法は高速なstepwise法
# サーチの経過は表示しない
auto.arima(ts_y[, j], ic = "aic", stepwise = TRUE, trace = FALSE)
})
# 列名をリストの要素名にする
names(r_arima) <- colnames(ts_y)
# リストの要素名を確認
names(r_arima)[1] "JTAR" "UR" [1] "coef" "sigma2" "var.coef"
[4] "mask" "loglik" "aic"
[7] "arma" "residuals" "call"
[10] "series" "code" "n.cond"
[13] "nobs" "model" "bic"
[16] "aicc" "x" "fitted" Series: ts_y[, j]
ARIMA(3,1,0)
Coefficients:
ar1 ar2 ar3
0.6120 0.1773 0.0503
s.e. 0.0545 0.0640 0.0552
sigma^2 = 0.0001248: log likelihood = 1031.24
AIC=-2054.49 AICc=-2054.37 BIC=-2039.234.2 予測値
結果オブジェクトから予測値を簡単に出せます。
Point Forecast Lo 95
Jan 2026 1.205617 1.1837207
Feb 2026 1.210828 1.1692912
Mar 2026 1.215515 1.1525067
Apr 2026 1.219589 1.1337814
May 2026 1.223176 1.1138659
Jun 2026 1.226328 1.0931683
Jul 2026 1.229098 1.0719939
Aug 2026 1.231533 1.0505740
Sep 2026 1.233672 1.0290804
Oct 2026 1.235552 1.0076428
Nov 2026 1.237205 0.9863585
Dec 2026 1.238656 0.9653003
Jan 2027 1.239932 0.9445217
Feb 2027 1.241054 0.9240618
Mar 2027 1.242039 0.9039476
Apr 2027 1.242905 0.8841973
May 2027 1.243666 0.8648221
Jun 2027 1.244335 0.8458274
Jul 2027 1.244923 0.8272144
Aug 2027 1.245439 0.8089810
Sep 2027 1.245893 0.7911224
Oct 2027 1.246292 0.7736320
Nov 2027 1.246643 0.7565017
Dec 2027 1.246951 0.7397223
Hi 95
Jan 2026 1.227513
Feb 2026 1.252364
Mar 2026 1.278523
Apr 2026 1.305397
May 2026 1.332485
Jun 2026 1.359489
Jul 2026 1.386203
Aug 2026 1.412492
Sep 2026 1.438264
Oct 2026 1.463462
Nov 2026 1.488051
Dec 2026 1.512013
Jan 2027 1.535343
Feb 2027 1.558046
Mar 2027 1.580131
Apr 2027 1.601613
May 2027 1.622510
Jun 2027 1.642843
Jul 2027 1.662631
Aug 2027 1.681898
Sep 2027 1.700664
Oct 2027 1.718952
Nov 2027 1.736784
Dec 2027 1.754179forecastの戻り値のオブジェクトの中身は予測値ですが、それに対応したplotメソッドがあるので、そのままプロットもできます。
4.3 欠損値処理
欠損値が含まれている場合は内部的には状態空間モデルを用いるので、問題なく計算できます。
ts_y_missing <- ts_y[, "JTAR"]
ts_y_missing[sample(1:(length(ts_y_missing)), 100)] <- NA
r_arima_missing <- auto.arima(ts_y_missing, ic="aic", stepwise = TRUE, trace = FALSE)
r_forecast_missing <- forecast(r_arima_missing, h = 24, level = 95)
plot(r_forecast_missing)
4.4 過剰適合
モデル選択と（ラグと階差以外の）パラメーターの推定を同時に行うので、データに過剰適合している蓋然性が高い事には注意してください。
4.5 外生変数
ARIMAは単一時系列分析ですが、外生変数を付け加えたARIMAX/SARIMAXモデルもあります。forecast::auto.arimaとforecast::forecastは、引数xregで行列形式の外部変数を与えることができます。
5 多変量時系列分析
多変量の時系列分析はベクトル自己回帰（VAR）モデルとそのバリエーションを用いることが多いです。説明変数の定常性を気にしつつ分析します。
5.1 単位根検定と階差
VARモデルにも単位根の問題はあるため、単位根検定で単位根が棄却できない場合は、階差をとったり、（データがすべて正であれば）Box-Cox変換をかけたりします。
単位根検定をする関数を提供するパッケージは幾つもありますが、urca::ur.dfを使ってみましょう。
###############################################
# Augmented Dickey-Fuller Test Unit Root Test #
###############################################
Test regression none
Call:
lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
Residuals:
Min 1Q Median 3Q
-0.060326 -0.005823 0.001114 0.007828
Max
0.028328
Coefficients:
Estimate Std. Error
z.lag.1 -0.0000446 0.0005913
z.diff.lag1 0.6187454 0.0550019
z.diff.lag2 0.1708425 0.0632770
z.diff.lag3 0.0496128 0.0546660
t value Pr(>|t|)
z.lag.1 -0.075 0.9399
z.diff.lag1 11.250 <2e-16 ***
z.diff.lag2 2.700 0.0073 **
z.diff.lag3 0.908 0.3648
---
Signif. codes:
0 '***' 0.001 '**' 0.01 '*' 0.05
'.' 0.1 ' ' 1
Residual standard error: 0.01109 on 328 degrees of freedom
Multiple R-squared: 0.639, Adjusted R-squared: 0.6346
F-statistic: 145.1 on 4 and 328 DF, p-value: < 2.2e-16
Value of test-statistic is: -0.0754
Critical values for test statistics:
1pct 5pct 10pct
tau1 -2.58 -1.95 -1.62
###############################################
# Augmented Dickey-Fuller Test Unit Root Test #
###############################################
Test regression none
Call:
lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
Residuals:
Min 1Q Median 3Q
-0.47396 -0.11693 -0.00575 0.11352
Max
0.56017
Coefficients:
Estimate Std. Error t value
z.lag.1 -0.00218 0.00245 -0.890
z.diff.lag1 -0.03058 0.05447 -0.561
z.diff.lag2 -0.09859 0.05391 -1.829
z.diff.lag3 -0.14028 0.05414 -2.591
Pr(>|t|)
z.lag.1 0.37431
z.diff.lag1 0.57489
z.diff.lag2 0.06838 .
z.diff.lag3 0.00999 **
---
Signif. codes:
0 '***' 0.001 '**' 0.01 '*' 0.05
'.' 0.1 ' ' 1
Residual standard error: 0.1764 on 328 degrees of freedom
Multiple R-squared: 0.03226, Adjusted R-squared: 0.02045
F-statistic: 2.733 on 4 and 328 DF, p-value: 0.02906
Value of test-statistic is: -0.8896
Critical values for test statistics:
1pct 5pct 10pct
tau1 -2.58 -1.95 -1.62lagsはラグです。厳密にはVARモデルのラグと合わせるべきですが、とりあえず3としました。先行研究や何かの理屈で正当化するか、情報量規準で機械的に選ぶかする方が良いです。なお、ドリフト項やトレンド項を仮定する場合は、typeで指定します。
どちらも単位根を棄却できませんでした。Value of test-statisticの絶対値が、Critical values for test statisticsの絶対値を超えていないです。
一階の差分をとって、もう一度単位根検定をしてみましょう。
###############################################
# Augmented Dickey-Fuller Test Unit Root Test #
###############################################
Test regression none
Call:
lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
Residuals:
Min 1Q Median 3Q
-0.060424 -0.005932 0.001104 0.007852
Max
0.027575
Coefficients:
Estimate Std. Error t value
z.lag.1 -0.16112 0.03724 -4.327
z.diff.lag1 -0.21776 0.05890 -3.697
z.diff.lag2 -0.05373 0.05868 -0.916
z.diff.lag3 -0.00846 0.05473 -0.155
Pr(>|t|)
z.lag.1 2.01e-05 ***
z.diff.lag1 0.000256 ***
z.diff.lag2 0.360523
z.diff.lag3 0.877243
---
Signif. codes:
0 '***' 0.001 '**' 0.01 '*' 0.05
'.' 0.1 ' ' 1
Residual standard error: 0.0111 on 327 degrees of freedom
Multiple R-squared: 0.1439, Adjusted R-squared: 0.1334
F-statistic: 13.74 on 4 and 327 DF, p-value: 2.292e-10
Value of test-statistic is: -4.3266
Critical values for test statistics:
1pct 5pct 10pct
tau1 -2.58 -1.95 -1.62
###############################################
# Augmented Dickey-Fuller Test Unit Root Test #
###############################################
Test regression none
Call:
lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
Residuals:
Min 1Q Median 3Q
-0.50922 -0.13104 -0.01825 0.09651
Max
0.58542
Coefficients:
Estimate Std. Error t value
z.lag.1 -1.52761 0.11774 -12.974
z.diff.lag1 0.46713 0.09668 4.832
z.diff.lag2 0.34580 0.07633 4.530
z.diff.lag3 0.20066 0.05370 3.737
Pr(>|t|)
z.lag.1 < 2e-16 ***
z.diff.lag1 2.08e-06 ***
z.diff.lag2 8.27e-06 ***
z.diff.lag3 0.00022 ***
---
Signif. codes:
0 '***' 0.001 '**' 0.01 '*' 0.05
'.' 0.1 ' ' 1
Residual standard error: 0.1732 on 327 degrees of freedom
Multiple R-squared: 0.5424, Adjusted R-squared: 0.5368
F-statistic: 96.89 on 4 and 327 DF, p-value: < 2.2e-16
Value of test-statistic is: -12.9745
Critical values for test statistics:
1pct 5pct 10pct
tau1 -2.58 -1.95 -1.62今度は単位根を棄却できました。
5.2 VARモデル
時系列が2本あるときのVARモデルは以下のようになります。
\[\begin{align} \begin{pmatrix} y_{1,t} \\ y_{2,t} \end{pmatrix} = \begin{pmatrix} \tau_{1} \\ \tau_{2} \end{pmatrix} + \sum^{L}_{\ell=1} \begin{pmatrix} \alpha_{1,1,\ell} & \alpha_{1,2,\ell} \\ \alpha_{2,1,\ell} & \alpha_{2,2,\ell} \end{pmatrix} \begin{pmatrix} y_{1,t - \ell} \\ y_{2,t - \ell} \end{pmatrix} + \begin{pmatrix} \epsilon_{1,t} \\ \epsilon_{2,t} \end{pmatrix} \end{align}\]
誤差項の共分散を考慮する表記もあるのですが、ここでは使用するパッケージにあわせて\(Cov(\epsilon_{1,t}, \epsilon_{2,t})=0\)を仮定します。また、誤差項のラグは入りません。
行列を展開すると、\(y_{1,t}\)が\(y_{2,t - \ell}\)に影響され、\(y_{2,t}\)が\(y_{1,t - \ell}\)に影響される系列間の相互作用があることがわかります。
内部的にはOLSなのですが、後々の事を考えてvarsパッケージを用いて推定します。
library(vars)
y <- with(df01, cbind(JTAR, UR))
dy <- diff(y, d)
(r_VAR <- VAR(dy, season = 12, p = L, type = "none"))
VAR Estimation Results:
=======================
Estimated coefficients for equation JTAR:
=========================================
Call:
JTAR = JTAR.l1 + UR.l1 + JTAR.l2 + UR.l2 + JTAR.l3 + UR.l3 + sd1 + sd2 + sd3 + sd4 + sd5 + sd6 + sd7 + sd8 + sd9 + sd10 + sd11
JTAR.l1 UR.l1
6.350202e-01 -7.683781e-04
JTAR.l2 UR.l2
1.850550e-01 -1.081519e-03
JTAR.l3 UR.l3
2.513386e-02 4.045717e-03
sd1 sd2
5.795690e-04 7.888610e-03
sd3 sd4
7.802413e-04 -6.596309e-05
sd5 sd6
5.396528e-03 4.708628e-03
sd7 sd8
1.879278e-03 5.199467e-03
sd9 sd10
5.821393e-03 2.106307e-03
sd11
5.383255e-03
Estimated coefficients for equation UR:
=======================================
Call:
UR = JTAR.l1 + UR.l1 + JTAR.l2 + UR.l2 + JTAR.l3 + UR.l3 + sd1 + sd2 + sd3 + sd4 + sd5 + sd6 + sd7 + sd8 + sd9 + sd10 + sd11
JTAR.l1 UR.l1 JTAR.l2
-0.24407047 -0.25872880 -2.32061641
UR.l2 JTAR.l3 UR.l3
-0.23332457 -0.41459354 -0.06044361
sd1 sd2 sd3
-0.09307099 0.07712959 -0.04889832
sd4 sd5 sd6
-0.21233242 -0.27812709 -0.26956426
sd7 sd8 sd9
-0.16107285 -0.15464012 -0.19255033
sd10 sd11
-0.28582985 -0.33773954 月次データなのでseasonは12とし、データ全体を見ると単調増加や単調減少ではないのでtypeはnoneでトレンドなどは無いものと見做します。
summary(r_VAR)とすると、係数の標準誤差とそのt検定の結果などが表示されますが、ラグや相互作用があるため係数間の関係は複雑で、直感的に結果の意味は把握はできません。インパルス応答列を見たり、グレンジャー因果性テスト、予測誤差分散分解をかけることになります。
5.2.1 情報量規準によるラグの選択
3としたラグですが、情報量規準で選ぶこともできます。
$selection
AIC(n) HQ(n) SC(n) FPE(n)
2 2 2 2
$criteria
1 2
AIC(n) -1.290409e+01 -1.299117e+01
HQ(n) -1.277336e+01 -1.284177e+01
SC(n) -1.257661e+01 -1.261692e+01
FPE(n) 2.488135e-06 2.280740e-06
3 4
AIC(n) -1.297312e+01 -1.297522e+01
HQ(n) -1.280505e+01 -1.278847e+01
SC(n) -1.255208e+01 -1.250740e+01
FPE(n) 2.322439e-06 2.317768e-06
5 6
AIC(n) -1.297310e+01 -1.295181e+01
HQ(n) -1.276768e+01 -1.272771e+01
SC(n) -1.245850e+01 -1.239042e+01
FPE(n) 2.322935e-06 2.373225e-06
7 8
AIC(n) -1.293016e+01 -1.292211e+01
HQ(n) -1.268739e+01 -1.266066e+01
SC(n) -1.232199e+01 -1.226716e+01
FPE(n) 2.425529e-06 2.445566e-06
9 10
AIC(n) -1.290882e+01 -1.290270e+01
HQ(n) -1.262869e+01 -1.260390e+01
SC(n) -1.220708e+01 -1.215419e+01
FPE(n) 2.478786e-06 2.494564e-06
11 12
AIC(n) -1.289539e+01 -1.290451e+01
HQ(n) -1.257792e+01 -1.256836e+01
SC(n) -1.210010e+01 -1.206243e+01
FPE(n) 2.513521e-06 2.491445e-06単位根検定にラグが必要なので、厳密には、この結果のラグ2で単位根検定をやり直す必要があります。もっとも単位根検定の結果はだいたい階差の次数で決まっており、ラグの長さはさほど影響しないです。あまり気にしなくて大丈夫です。もちろん、なるべく厳密にしたい時は、階差とラグごとに単位根を棄却できるか表をつくり、その中から情報量規準で選びましょう。
5.2.2 制約付きVARモデル
OLSなので制約付きモデルをつくる事も可能です。
# t値の閾値2.0でゼロ制約をつける係数を定める（これ自体はbad-practice）
r_restrict_auto <- restrict(r_VAR, method = "ser", thresh = 2.0)統計解析としてはbad-praticeなのですが、コーディング時の工夫になります。つまり、結果オブジェクトを見るとどういう構造の行列を渡せばよいかわかります。
JTAR.l1 UR.l1 JTAR.l2 UR.l2
JTAR 1 0 1 0
UR 0 1 1 1
JTAR.l3 UR.l3 sd1 sd2 sd3 sd4 sd5
JTAR 0 0 0 0 0 0 0
UR 0 0 1 0 1 1 1
sd6 sd7 sd8 sd9 sd10 sd11
JTAR 0 0 0 0 0 0
UR 1 1 1 1 1 1行が系列を表し、列が係数を表していますね。
ラグ3の係数を0にして、他は1とする制約をつくり推定しましょう。
restrictions[, ] <- 1
restrictions[, c("JTAR.l3", "UR.l3")] <- 0
# ゼロ制約を表す行列を指定してVARを行う
r_restrict_manual <- restrict(r_VAR, method = "man", resmat = restrictions)制約付きモデルの推定結果が計算できました。
なお、制約を表す行列は
restrictions <- matrix(1, length(coef(r_VAR)), nrow(coef(r_VAR)[[1]]), dimnames = list(names(coef(r_VAR)), rownames(coef(r_VAR)[[1]])))としてから、任意の場所に0を入れてもよいです。
5.2.3 グレンジャー因果
グレンジャー因果は、制約付きモデルと制約なしモデルの比較の応用です。 ある系列の予測値を、原因となっていそうな系列を説明変数に含めたモデルと、含めないモデルでそれぞれ出し、二つの差を検定にかけます。 予測精度が統計的に有意に上がるのであれば因果がありそうだと考えます。
Granger causality H0: JTAR do not
Granger-cause UR
data: VAR object r_VAR
F-Test = 14.32, df1 = 3, df2 = 630,
p-value = 4.9e-09
Granger causality H0: UR do not
Granger-cause JTAR
data: VAR object r_VAR
F-Test = 0.33845, df1 = 3, df2 =
630, p-value = 0.7975基本的なVARでは共分散構造は考えないので、Instantaneous causalityは参照していません。
変数JTARは有効求人倍率（jobs-to-applicants
ratio）は変数URは失業率（unemployment
ratio）で、有効求人倍率は景気一致指数、失業率は景気遅行指数として知られていますが、有効求人倍率は失業率のグレンジャー因果である一方、失業率は有効求人倍率とは言えない事が確認されました。
5.2.4 予測値
階差をとった値の予測になるので、直感的な把握はしづらいですが予測値が出ます。
diffinvで差分を足した方が把握しやすいかもです。
5.2.5 インパルス応答列
VARの推定結果の係数からは、ある変数が増えた時に、他の変数が増えるのか減るのか良く分かりません。インパルス応答列を見て解釈しましょう。
# 大きさ1のショックが生じたとして、18期先まで予測; 直交化は行わない
# 差分モデルなのでcumulative = TRUEをつけた方が、実データとの対応はわかりやすくなる
r_irf <- irf(r_VAR, n.ahead = n.ahead, ortho = FALSE)
par(mfrow = c(ncol(r_VAR$y), ncol(r_VAR$y)))
plot(r_irf, plot.type = "single")
plotの正体はplot.varirfで、?plot.varirfでヘルプを参照できます。
有効求人倍率に上昇ショックが起きると、完全失業率が減少することがわかりました。
5.3 構造化VARモデル
マクロ経済モデルでは、複数の変数が同時決定されることにより、撹乱項の共分散が非ゼロと予想される事があります。
これまでの有効求人倍率と完全失業率の分析では、異なる系列間の作用は一ヶ月後からと仮定してきましたが、もっと同時に生じて、撹乱項の共分散が非ゼロと言う事もあります。
\[\begin{align} \begin{pmatrix} y_{1,t} \\ y_{2,t} \end{pmatrix} &= \begin{pmatrix} \tau_{1} \\ \tau_{2} \end{pmatrix} - \begin{pmatrix} 0 & \alpha_{1,2,0} \\ \alpha_{2,1,0} & 0 \end{pmatrix} \begin{pmatrix} y_{1,t} \\ y_{2,t} \end{pmatrix} + \sum^{L}_{\ell=1} \begin{pmatrix} \alpha_{1,1,\ell} & \alpha_{1,2,\ell} \\ \alpha_{2,1,\ell} & \alpha_{2,2,\ell} \end{pmatrix} \begin{pmatrix} y_{1,t - \ell} \\ y_{2,t - \ell} \end{pmatrix} + \begin{pmatrix} \epsilon_{1,t} \\ \epsilon_{2,t} \end{pmatrix} \\ \begin{pmatrix} 1 & \alpha_{1,2,0} \\ \alpha_{2,1,0} & 1 \end{pmatrix} \begin{pmatrix} y_{1,t} \\ y_{2,t} \end{pmatrix} &= \begin{pmatrix} \tau_{1} \\ \tau_{2} \end{pmatrix} + \sum^{L}_{\ell=1} \begin{pmatrix} \alpha_{1,1,\ell} & \alpha_{1,2,\ell} \\ \alpha_{2,1,\ell} & \alpha_{2,2,\ell} \end{pmatrix} \begin{pmatrix} y_{1,t - \ell} \\ y_{2,t - \ell} \end{pmatrix} + \begin{pmatrix} \epsilon_{1,t} \\ \epsilon_{2,t} \end{pmatrix} \end{align}\]
煩雑なので、書き直します。
\[ A \mathbf{y}_t = \mathbf{\tau}^* + \sum^{L}_{\ell=1} A^*_\ell \mathbf{y}_{t - \ell} + \mathbf{\epsilon}_t \]
構造化VAR（SVAR）のAモデルと呼ばれる形式になりました。
\(B=A^{-1}\)を両辺に生じると、Bモデルになります。\(A^{-1}\mathbf{\tau}^*=\mathbf{\tau}\), \(A^{-1}A^*_\ell=A_\ell\)になっているので注意してください。
\[ \mathbf{y}_t = \mathbf{\tau} + \sum^{L}_{\ell=1} A_\ell \mathbf{y}_{t - \ell} + B \mathbf{\epsilon}_t \]
両方を組み合わせたABモデルもあります。
これらは最尤法で推定しますが、vars:SVARが計算してくれるので尤度関数を考える必要はありません。しかし、AやBは何でも良いと言うことはなく、推定可能なように指定する必要があります。
library(vars)
d <- 1; L <- 3; n.ahead <- 18
y <- with(df01, cbind(JTAR, UR))
dy <- diff(y, d)
Amat <- matrix(c(1, NA, 0, 1), 2, 2) # 推定する箇所にNAを入れる;NAの数はK*(K-1)以下
r_VAR <- VAR(dy, season = 12, p = L, type = "none")
# このデータではデフォルトのestmethod = "scoring"では収束しないためdirect/BFGSを指定した
(r_SVAR <- SVAR(r_VAR, estmethod = "direct", Amat = Amat, max.iter = 1000, method="BFGS"))
SVAR Estimation Results:
========================
Estimated A matrix:
JTAR UR
JTAR 1.000 0
UR 1.157 1# インパルス応答列
r_irf_svar <- irf(r_SVAR, n.ahead = n.ahead, ortho = FALSE)
par(mfrow = c(ncol(r_VAR$y), ncol(r_VAR$y)))
plot(r_irf_svar, plot.type = "single")
有効求人倍率の効果が大きく計測されるようになりました。
5.4 共和分検定
差分を取ると単位根がなくなる二つの時系列\(y_1\)と\(y_2\)の間に、定数\(k\)と定常過程\(z_t\)で表せる、\(y_{1t} - k y_{2t} = z_t\)という関係があるとき、\(y_1\)と\(y_2\)は共和分関係にあるといいます。また、共和分関係があるときは、VARではモデル特定化が誤っている事になります。
urca::ca.poで共和分テストができるのでやってみましょう。
########################################
# Phillips and Ouliaris Unit Root Test #
########################################
Test of type Pu
detrending of series none
Call:
lm(formula = z[, 1] ~ z[, -1] - 1)
Residuals:
Min 1Q Median 3Q Max
-0.7677 -0.4035 0.1230 0.6559 1.1291
Coefficients:
Estimate Std. Error t value
z[, -1] 0.217771 0.007964 27.34
Pr(>|t|)
z[, -1] <2e-16 ***
---
Signif. codes:
0 '***' 0.001 '**' 0.01 '*' 0.05
'.' 0.1 ' ' 1
Residual standard error: 0.5763 on 335 degrees of freedom
Multiple R-squared: 0.6906, Adjusted R-squared: 0.6897
F-statistic: 747.7 on 1 and 335 DF, p-value: < 2.2e-16
Value of test-statistic is: 0.989
Critical values of Pu are:
10pct 5pct 1pct
critical values 20.3933 25.9711 38.3413今回はトレンドなし過程を仮定しているので、demean = "none"にしています。
Value of test-statisticの絶対値がCritical valuesの10pct,
5pct, 1pct絶対値を超えたら、それぞれ10%, 5%,
1%有意で共和分が無い（\(z_t\)に単位根がない）という帰無仮説を棄却します。
今回は共和分がないと言う帰無仮説（\(z_t\)が単位根過程）は棄却されませんでした。
5.5 VECMの推定
共和分関係がある場合は、ベクトル誤差修正モデル（VECM）を使う方が望ましいと言うことになります。
VARをやや抽象的に表記すると以下のようになります。
\[ \mathbf{y}_t = \mathbf{\tau} + \sum^{L}_{\ell=1} A_\ell \mathbf{y}_{t - \ell} + \mathbf{\epsilon}_t \]
\(\mathbf{y}\)は時系列変数のベクトル、\(\mathbf{\tau}\)はトレンド項、\(A\)は係数行列、\(\mathbf{\epsilon}\)は撹乱項のベクトルです。
定常過程にするために、1期前の値を引きましょう。
\[\begin{align} \mathbf{y}_t - \mathbf{y}_{t-1} &= \mathbf{\tau} + \sum^{L}_{\ell=1} A_\ell \mathbf{y}_{t - \ell} + \mathbf{\epsilon}_t - \mathbf{y}_{t-1} \\ &= \mathbf{\tau} + (\sum^L_{\ell=1} A_{\ell} - I)\mathbf{y}_{t-1} + \sum^{L-1}_{i=1} \sum^{L}_{j=i+1} A_j (\mathbf{y}_{t - i} - \mathbf{y}_{t - i - 1}) + \mathbf{\epsilon}_t \end{align}\]
\(\Pi = \sum^L_{\ell=1} A_{\ell} - I\)、\(\Gamma_i = \sum^{L}_{j=i+1} A_j\)と置いてまとめると、
\[\begin{align} \Delta \mathbf{y}_t &= \mathbf{\tau} + \Pi \mathbf{y}_{t-1} + \sum^{L-1}_{i=1} \Gamma_i \Delta \mathbf{y}_{t - i} + \mathbf{\epsilon}_t \end{align}\]
となります。
階差をとれば定常過程だと分かっているので、左辺も右辺も定常です。右辺の第1項は定数、第3項は階差に行列を乗じており、第4項は撹乱項なので定常です。右辺の第2項は定常でなければなりません。
右辺の第2項が定常になるためには、\(\Pi\)がゼロ行列で階差をとったVARモデルと同様になるか、共和分関係があるかです。\(\Pi \mathbf{y}_{t-1}\)の\(i\)行目は、\(\pi_{ij}\)を\(\Pi\)の成分として、\(\pi_{i1} y_{1,t-1} + \pi_{i2} y_{2,t-1} + \cdots + \pi_{iL} y_{L,t-1}\)となっており、これが定常であることが、すなわち共和分関係があることになるからです。
\(\Pi\)は共和分ベクトルそのままではなく、推定するパラメーターである調整係数行列を乗じたものになります。時系列2つの例で、上の第2項を書き直すと
\[ \Pi \mathbf{y}_{t-1} = \begin{pmatrix} lv_1 \\ lv_2 \end{pmatrix} \begin{pmatrix} cv_1 & cv_2 \end{pmatrix} \begin{pmatrix} y_{1, t-1} \\ y_{2, t-1} \end{pmatrix} \]
となります。\(lv_i\)は調整係数行列の要素である送出ベクトル（loading vector）の成分で、\(ci_j\)は共和分ベクトル（cointegration vector）の成分です。共和分ベクトルの数は時系列の数より少なくなりますが、調整係数行列を乗じることで\(\Pi\)の行数を時系列の数と同じにしています。
VARに\(\Pi\)を付け加えたこのモデルがVECMになります。推定は\(\Pi\)があるために煩雑になります[^fn03]が、川崎 (1993)の説明が参考になると思います。
varsパッケージでは、urca::ca.jo関数で推定した後に、vars::vec2varで結果オブジェクトを変換して、インパルス応答列と予測誤差分散分解を見る事ができます。
r_ca.jo <- ca.jo(y, ecdet = "none", K = L, season = 12)
r_vecm <- vec2var(r_ca.jo)
r_vecm_irf <- irf(r_vecm, n.ahead = n.ahead, ortho = FALSE)
par(mfrow = c(ncol(r_vecm$y), ncol(r_vecm$y)))
plot(r_vecm_irf, plot.type = "single")
予測誤差分散分解（FEVD）の結果が大きく変わり、有効求人倍率が大きな予測因子になっています。共和分検定で帰無仮説を棄却できていないので、VARの方がそれらしいと言う事になりますが。
5.6 構造化VECM
共和分関係について制約を置く構造化VECM（SVECM）もあります。
6 状態空間モデル
被説明変数が離散もしくは非負である場合、観測値に予測区間をつける必要がある場合、推定を逐次更新したい場合は、状態空間モデルの利用を考えましょう。
7 終わりに
教科書的な手法を網羅的に見てきました。実務的には無目的に収集した時系列データがあって、ここから何か分からないかと言う話になることは多々あるのですが、そういうときに何かする助けになると思います。ただし、モデルの階差とラグの次数をデータから決めるようなことになりがちで、往々にしてデータに過剰適合した推定結果になることには注意してください。また、見ての通り予測値の幅は広くなるため、データが助けにならない事を確認して終わる事が多そうです。