時系列分析（ARIMA）
東京国際大学 データサイエンス教育研究所 竹田 恒
2023-12-05
データ
月間航空旅客数のデータを使用する。
y <- c(AirPassengers) # 月間航空旅客数
n <- length(y) # 月数
n.tr <- 120 # 訓練月数
n.te <- n - n.tr # テスト月数
t <- 1:n # 月
ii.tr <- 1:120 # 訓練期間インデックス
ii.te <- 121:n # テスト期間インデックス
matplot(y, type = 'o', pch = 1, col = 2)
fq <- 12グラフから周期（\(fq = 12\)）の周期性があることが分かる。
時間の経過とともに上昇するトレンドが有り振幅も徐々に大きくなっていることから，この時系列は非定常であることが分かる。 定常時系列にするためにはトレンド除去や変数変換が必要となる。
定常化（BoxCox変換）
BoxCox変換は次式で\(x\)を変換し\(y^{(\lambda)}\)を得る変数変換法である。 \(\lambda\)が0のときは定義により（自然）対数変換となる。
\[y^{(\lambda)} = \begin{cases}\frac{x^\lambda-1}{\lambda}\quad (\lambda \ne 0)\\ \ln(x)\quad (\lambda = 0)\end{cases}\]
library(MASS)
bc <- boxcox(y ~ t, data = data.frame(y = y[ii.tr], t = t[ii.tr]))
lambda <- round(bc$x[which.max(bc$y)], 1)
lambda## [1] 0このデータでは\(\lambda\)は0のとき対数尤度（たいすうゆうど；log-likelihood）が最大になっている。
library(forecast)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zooy.bc <- BoxCox(y[ii.tr], lambda)
matplot(y.bc, type = 'o', pch = 1, col = 2, ylab = 'y')
振幅がほぼ一定になった。
定常化（トレンド除去）
1次のトレンド（直線的傾向変動）があることから， 1次の階差をとりトレンド除去する。 2次のトレンド（曲線的傾向変動）がある場合は2次階差をとるなど， トレンドがなくなるまで階差を取り続ける（階差を取ることは微分することに相当）。
d <- 1
dy <- diff(y.bc, diff = d)
matplot(dy, type = 'o', pch = 1, col = 2, ylab = 'y')
1次の階差をとることで，トレンドが消えほぼ定常化できた。
定常化のためにBoxCox変換（\(\lambda = 0\)）と1次階差（\(d = 1\)）が必要となる ので，ARIMAモデルをフィッティングする関数（forecast::Arima）のオプションを設定する。 なお，1次階差をとる前の時系列は定常時系列から見ると和分（差分の反対） した時系列となる。
ACF・PACF
自己相関関数（ACF：auto-correlation function）
原系列とラグ系列の自己相関
偏自己相関関数（PACF：partial auto-correlation function）
原系列とラグ系列の本質的な自己相関（自己相関から他の短いラグ系列の影響を取り除いたもの）。
最大次数が分かる。
AR，MAの次数p，qをACF，PACFグラフを見て決める。
| モデル | ACF | PACF |
|---|---|---|
| AR | 減衰 | 有意なp次ラグ |
| MA | 有意なq次ラグ | 減衰 |
| ARMA | 減衰 | 減衰 |
ys <- ts(10 * sin(2 * pi * t / 10) + rnorm(n))
dir.create('acf_lag', showWarnings = F)
for (k in 1:25)
{
png(paste0('acf_lag/acf_lag', k, '.png'))
ys.lag <- lag(ys, k = -k)
ts2 <- cbind(ys, ys.lag)
par(mfrow = c(3, 1), mar = c(2, 4, 1, 1))
matplot(ts2, type = 'o', pch = 1, col = c(1, 4),
xlab = '', ylab = '', xlim = c(0, 120))
abline(v = c(0, k), col = 2)
arrows( 0, 0, k, 0, length = 0.1, col = 2)
mtext(paste(k), side = 1, at = k, line = 1, col = 2)
legend('topright', col = c(1, 4), bg = 'white', lty = 1,
legend = c('原系列（周期10の正弦波＋ノイズ）', paste0(k, '次ラグ系列')))
grid()
acf(ys, lag.max = k, xlim = c(0, 120), ylim = c(-1, 1), main = '')
pacf(ys, lag.max = k, xlim = c(0, 120), ylim = c(-1, 1), main = '')
dev.off()
}
#system('bash -c 'convert -delay 100 -loop 0 ./acf_lag/acf_lag{1..25}.png ./acf_lag/acf_lag.gif'')\[Y_t = 10\sin\left(2\pi \frac{t}{10}\right) + \epsilon_t, \quad \epsilon_t \sim \mathbf{N}(0, 1)\]
【注意】青点線で囲まれた範囲は有意ではない。
ACF
周期性がある場合は周期内で有意な最大ラグ次数を選ぶ。
a <- acf(dy, lag.max = 50) # MA(q)
grid()
abline(v = fq, col = 3, lty = 3)
q <- 12 PACF
pa <- pacf(dy, lag.max = 50) # AR(p)【注意】acfと違いグラフが1次ラグから始まる。
grid()
abline(v = fq, col = 3, lty = 3)
p <- 12フィッティング
Arima関数のオプション設定 order = c(p, d,
q)にAR，和分，MAの次数p，d，qを与える。
frequencyには季節周期を入力する。
seasonalには季節周期を単位時間とみなしたときの周期点系列のモデル次数を設定する。
変数変換（BoxCox変換）を自動で行うオプション設定（lambda =
‘auto’）を使用。
【注意】次数が大きいと計算時間が掛かる。また，訓練データは大きなサイズが必要となる。
fit.arima <- Arima(ts(y[ii.tr], frequency = fq), lambda = 'auto',
order = c(p, d, q),
seasonal = c(1, 1, 0)) # SARIMA
# c.f.
#fit.arima <- auto.arima(ts(y[ii.tr], frequency = fq), lambda = 'auto')
fit.arima## Series: ts(y[ii.tr], frequency = fq)
## ARIMA(12,1,12)(1,1,0)[12]
## Box Cox transformation: lambda= -0.3096643
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 ar6 ar7 ar8
## -0.2542 0.0940 -0.0545 0.0567 -0.0594 -0.0297 0.1077 -0.0028
## s.e. 0.1417 0.1441 0.1368 0.1352 0.1335 0.1306 0.1425 0.1334
## ar9 ar10 ar11 ar12 ma1 ma2 ma3 ma4
## 0.1405 0.0729 -0.2619 -0.2848 -0.0736 -0.1363 -0.0955 -0.2205
## s.e. 0.1335 0.1311 0.1333 0.1807 0.1472 0.1384 0.1264 0.1472
## ma5 ma6 ma7 ma8 ma9 ma10 ma11 ma12
## 0.1913 0.2118 -0.1829 -0.0367 -0.1446 0.0096 0.3593 -0.6336
## s.e. 0.1458 0.1479 0.1603 0.1385 0.1426 0.1516 0.1519 0.1503
## sar1
## 0.1202
## s.e. 0.2297
##
## sigma^2 = 5.446e-05: log likelihood = 381.33
## AIC=-710.67 AICc=-693.12 BIC=-641.17残差分析
ACFで有意なものがない（棒が青点線をでていない）こと， 残差分布が0点回りに左右対称（釣り鐘状）になっていることを確認する。
Ljung-Box（リュング・ボックス）検定は残差系列の自己相関の有無を判定する。
p値が0.05以下であると5％水準で有意な自己相関がある。
帰無仮説：残差は独立に分布している（ホワイトノイズ）。
対立仮設：残差には自己相関がある。
残差分析の結果が思わしくないときはモデルを見直す。 また，気温などの外的要因の影響があるときは，外生変数（説明変数；xregオプションで与える）の導入を検討する。
checkresiduals(fit.arima) # 残差分析グラフ
##
## Ljung-Box test
##
## data: Residuals from ARIMA(12,1,12)(1,1,0)[12]
## Q* = 6.3161, df = 3, p-value = 0.09721
##
## Model df: 25. Total lags used: 28予測
pred <- forecast(fit.arima, h = n.te, level = 95)グラフ
matplot(t, y, type = 'o', pch = 16, col = 2, main = pred$method)
grid()
abline(v = n.tr + 1, lty = 3, col = 3)
matlines(t[ii.tr], fit.arima$fitted, type = 'o', pch = 16, lty = 2, col = 4)
matlines(t[ii.te], pred$mean, type = 'o', pch = 1, lty = 2, col = 4)
matlines(t[ii.te], pred$lower, lty = 3, col = 'gray')
matlines(t[ii.te], pred$upper, lty = 3, col = 'gray')
legend('topleft', bg = 'white', lty = c(1, 2, 2, 3),
col = c(2, 4, 4, 'gray'), pch = c(16, 16, 1, NA),
legend = c('原系列', '1期先予測値', '予測値', '95%予測区間'))
精度評価
get.accuracy <- function(yhat, y)
{
data.frame(MBE = mean(yhat - y),
MAE = mean(abs(yhat - y)),
MAPE = mean(abs((yhat - y) / y)) * 100,
RMSE = sqrt(mean((yhat - y)^2))
)
}
get.accuracy(yhat = pred$mean, y = y[ii.te])## MBE MAE MAPE RMSE
## 1 -14.07091 14.92617 3.421171 17.18222周期性を有する時系列の長期予測モデル選定における注意点
Rの時系列解析ライブラリforecastには，auto.arimaと呼ばれる自動でARIMAモデル
の選択，フィッティングを行う関数がある。
この関数の利用は周期性を有する時系列の長期予測には以下の理由により推奨しない。
AICは期待平均対数尤度の不偏推定量であり， 最大対数尤度とパラメータ数を明示的に含んだ関数である。 AICではパラメータ数は罰則として機能しモデル選定でよく使用される。 重回帰モデルとは違い，時系列モデルにおいてはARの最尤計算などで 使用する誤差は1期先予測誤差であり， 本来時系列に内包する周期性を表現するため高次のモデルが必要な場合であっても， 低次のモデルでも1期先予測誤差は十分小さくなる場合がある。このため， モデル選択は比較的低次のモデルになりやすいという欠点がある。
周期性を有する時系列の長期予測のためには周期分の繰り返しを再現できる パラメータの数を十分確保できていることをよく確認すること。 長期予測のためには検証データ（validation data set）をテストデータとは 別に用意して，Arima関数を使って手動でパラメータ数を設定したほうが良い。 時系列データを訓練用，検証用，テスト用の3分割し，訓練と検証を繰り返し， 選定された最適なモデルで，最後にテストデータで精度評価する。 テストデータでパラメータ調整を行うと訓練データとテストデータに過適合するので， 必ず検証データで行うこと。
演習課題
次の時系列モデルで生成したデータを用いてARIMAモデルによる予測を行え。 （BoxCox変換やgif動画作成は実施しなくて良い。）
\[Y_t = 100 + 0.1t + 0.02t^2 + 100\sin\left(2\pi \frac{t}{24}\right) + \epsilon_t, \quad \epsilon_t \sim \mathbf{N}(0, 5^2)\]
n <- 24*7*2 # 時間数
t <- 1:n # 時間
n.tr <- floor(n * 0.8) # 訓練期間数
n.te <- n - n.tr # テスト期間数
ii.tr <- 1:n.tr # 訓練期間インデックス
ii.te <- (n.tr+1):n # テスト期間インデックス
set.seed(5)
y <- 100 + 0.1 * t + 0.02 * t^2 + 100 * sin(2 * pi * t / 24) + rnorm(n, sd = 5)
matplot(y, type = 'o', pch = 1, col = 2)