経営のための統計学II
第5回 時系列分析3
高橋 慎
2021/10/11
和文過程とARIMA過程
ここまでに紹介した分析手法は定常時系列に関するもので、定常性の条件を満たさない時系列については利用できない。
- 非定常時系列では、平均、分散、共分散が観測時点に依存して変化する。
非定常時系列は種種あるが、原系列の差分を取れば定常になる時系列を和文過程という。
- このような系列については、階差系列に関してARMA分析を行う。
1次の和文過程の最も簡単な系列は、\(\epsilon_{t}\)をホワイト・ノイズの系列として次のように定義される。 \[ X_{t} = X_{0} + \sum_{i=1}^{t} \epsilon_{i} \]
- \(X_{0}\)は時点0での初期値である。
- このような和文系列は定常性の性質を満たさないため、ARMA過程の分析ができない。
- 時系列が和文過程であるか否かを調べる必要があるが、そのための検定として、単位根検定が広く使われるようになった。
系列が和文過程である場合は、1次差分または階差を求めれば定常になる。
- 上記の\(X_{t}\)の階差を求めると \[ \Delta X_{t} = X_{t} - X_{t-1} = \epsilon_{t} \] となり、\(\Delta X_{t}\)は最も簡単な定常時系列であるホワイト・ノイズとなる。
ARMAなどの性質をもつ定常な時系列変数を\(u_{t}\)とすれば、和文過程は次のようになる。 \[ X_{t} = X_{0} + \sum_{i=1}^{t}u_{i} \]
- 階差を取ると、 \[ \Delta X_{t} = X_{t} - X_{t-1} = u_{t} \] となり、定常な時系列\(u_{t}\)になる。
- 次のステップとして、\(\Delta X_{t}\)に関して、すでに説明したようにARMAの次数を求める。
一般的には差分を一度でなく複数回取ることも可能である。
- \(X_{t}\)の階差を取ればARMA過程として分析できる時系列を、ARIMA(\(p,d,q\))過程（auto-regressive integrated moving average process：ARIMA）という。
- \(d\)は差分の回数であり、\(p\)と\(q\)は差分系列に関するARMAの次数である。
対数IIP（鉱工業生産指数）の推定
系列\(X_{t}\)の連続型成長率は \[ \Delta \log(X_t) = \log(X_t) - \log(X_{t-1}) \]
- これまでに行ってきたIIPの連続型成長率の推定は、対数IIPの階差に関する推定だと考えることができる。
もとの対数系列である対数IIPに関するARMA推定を行う。
- 対数IIP \(\log(IIP_{t})\)を省略して\(y_{t}\)と表記する。
# IIPデータの読み込み
iip <- read.csv("iip.csv")
# データをts型に変換
x <- ts(data = iip$iip, start = c(1978, 1), frequency = 12)
# 対数IIP
y <- log(x)まず、系列のプロットを調べ、次にACおよびPAC関数を計算する。
# 時系列プロット
plot(y)
# ACプロット
acf(y, lag.max = 20)
# PACプロット
pacf(y, lag.max = 20)
- AC関数が漸減していくのに比べ、PAC関数は1次しか有意でない。
- ACとPACの理論的な性質（前回の「次数の選択」を参照）から、仮に選ばれるモデルはARになる。
2次の自己回帰式AR(2)を推定する。
# AR(2)
ar2 <- arima(y, order = c(2, 0, 0))
ar2##
## Call:
## arima(x = y, order = c(2, 0, 0))
##
## Coefficients:
## ar1 ar2 intercept
## 0.9956 0.0005 4.4924
## s.e. 0.0456 0.0458 0.1503
##
## sigma^2 estimated as 0.000331: log likelihood = 1239.69, aic = -2471.38# 決定係数
1 - sum(ar2$residuals^2) / sum((y - mean(y))^2)## [1] 0.98428# AR多項式の根
root <- polyroot(c(1, ar2$coef[1], ar2$coef[2]))
root## [1] -1.004976+0i -1870.258347-0iround(abs(root), digits = 3) #根の絶対値## [1] 1.005 1870.258- AR(2)の2次のラグ係数は有意でない。
- 決定係数が異常に1に近いが、これがIIPあるいは対数IIPの原系列（つまり単位根過程が疑われる系列）に関する自己回帰式の特徴である。
- ARに関する多項式の根は1より大だが、そのうちの1つは1に非常に近い。
1次の自己回帰式AR(1)を推定する。
# AR(1)
ar1 <- arima(y, order = c(1, 0, 0))
ar1##
## Call:
## arima(x = y, order = c(1, 0, 0))
##
## Coefficients:
## ar1 intercept
## 0.9962 4.5034
## s.e. 0.0038 0.1513
##
## sigma^2 estimated as 0.000331: log likelihood = 1239.68, aic = -2473.37# 決定係数
1 - sum(ar1$residuals^2) / sum((y - mean(y))^2)## [1] 0.9842802# AR多項式の根
root <- polyroot(c(1, ar1$coef[1], ar1$coef[2]))
root## [1] -0.1106049+0.4580606i -0.1106049-0.4580606iround(abs(root), digits = 3) #根の絶対値## [1] 0.471 0.471- 決定係数は1に近くほぼ変わらないが、AICはAR(2)より小さい。
- AR多項式の根は1より大だが、1に非常に近い。
AR(1)の残差の診断を行う。
# 残差のプロット
par(family = "HiraginoSans-W3") #日本語が文字化けしないようにおまじない
ts.plot(ar1$residuals, gpars=list(xlab="年月", ylab="残差", lty=3, col=3))
# ２シグマ
twosigma <- 2*sqrt(ar1$sigma2)
twosigma## [1] 0.0363868# ２シグマを超える残差（異常値）のインデックスを抽出
i.outlier <- abs(ar1$residuals) > twosigma
# 年月ラベルを作成
date <- seq(as.Date("1978-01-01"), as.Date("2017-12-01"), by = "month")
# 年月と残差のデータフレーム
zansa <- data.frame(date = format(date, "%Y-%m"), residuals = ar1$residuals)
# ２シグマを超える年月と残差の値
zansa[i.outlier, ]## date residuals
## 114 1987-06 0.03661069
## 177 1992-09 0.04046167
## 183 1993-03 -0.04060958
## 277 2001-01 -0.04261244
## 293 2002-05 0.04285383
## 371 2008-11 -0.06839960
## 372 2008-12 -0.08548092
## 373 2009-01 -0.09272075
## 374 2009-02 -0.08968304
## 376 2009-04 0.04284169
## 399 2011-03 -0.17928316
## 401 2011-05 0.06610112
## 402 2011-06 0.04131928
## 436 2014-04 -0.04460743- 成長率のAR(2)やMA(2)と挙動はほとんど変わらない。
# 残差のACプロット
acf(ar1$residuals, lag.max = 20)
# 残差のPACプロット
pacf(ar1$residuals, lag.max = 20)
- AC(2)とPAC(2)が有意となっている。
# k=1から20までのLjung-Box検定統計量とp値
max.lag <- 20
Q_LB <- numeric(max.lag)
p_LB <- numeric(max.lag)
for (k in 1:max.lag) {
LB <- Box.test(ar1$residuals, k, type = "Ljung-Box", fitdf = 1)
Q_LB[k] <- LB$statistic
p_LB[k] <- LB$p.value
}
data.frame(k = 1:max.lag, Q_LB, p_LB)## k Q_LB p_LB
## 1 1 0.02690503 0.000000000
## 2 2 7.69832074 0.005527222
## 3 3 8.73319928 0.012694333
## 4 4 8.81862917 0.031802080
## 5 5 8.97638227 0.061692530
## 6 6 9.65628514 0.085581067
## 7 7 11.37075191 0.077571853
## 8 8 12.08831022 0.097691885
## 9 9 14.43406704 0.071129939
## 10 10 16.08337375 0.065160807
## 11 11 19.29231429 0.036703139
## 12 12 20.76660144 0.035871348
## 13 13 24.81121372 0.015743247
## 14 14 24.96282581 0.023344690
## 15 15 29.50449127 0.008924233
## 16 16 31.26380910 0.008097555
## 17 17 31.69484917 0.010954564
## 18 18 34.72818535 0.006752841
## 19 19 34.75518047 0.010145712
## 20 20 34.84077455 0.014596263- AR係数が1つなので、\(k=2\)で自由度が\(k-1=1\)になり、検定統計量として意味をもつ。
- \(p\)値はすべての次数で10%より小さいが、3次以降で1%より大きく、5次以降で5%より大きい。
- したがって、残差に系列相関がないという帰無仮説を強くは棄却できない。
最後に、forecastパッケージのauto.arima関数を用いて、ARIMA過程を含めて最小AIC規準の次数選択を行う。
# パッケージの読み込み
library(forecast)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo# ARMA過程の次数選択
fit <- auto.arima(y, max.p=2, max.q=2, seasonal = F, ic = "aic", stepwise = T, trace = T)##
## Fitting models using approximations to speed things up...
##
## ARIMA(2,1,2) with drift : -2469.902
## ARIMA(0,1,0) with drift : -2468.222
## ARIMA(1,1,0) with drift : -2465.587
## ARIMA(0,1,1) with drift : -2466.239
## ARIMA(0,1,0) : -2468.526
## ARIMA(1,1,2) with drift : -2470.828
## ARIMA(0,1,2) with drift : -2472.981
## ARIMA(1,1,1) with drift : -2466.421
## ARIMA(0,1,2) : -2473.576
## ARIMA(0,1,1) : -2466.53
## ARIMA(1,1,2) : -2471.521
## ARIMA(1,1,1) : -2466.611
##
## Now re-fitting the best model(s) without approximations...
##
## ARIMA(0,1,2) : -2481.59
##
## Best model: ARIMA(0,1,2)# 選択されたモデルの表示
summary(fit)## Series: y
## ARIMA(0,1,2)
##
## Coefficients:
## ma1 ma2
## -0.0148 0.1371
## s.e. 0.0455 0.0450
##
## sigma^2 estimated as 0.0003265: log likelihood=1243.79
## AIC=-2481.59 AICc=-2481.54 BIC=-2469.07
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.0009767195 0.01801383 0.01164648 0.02152651 0.2559383 0.2325809
## ACF1
## Training set 0.003128338- オプションで
d = 0と指定せずに、差分をとるARIMAモデルを含むように指定している。 - ARIMA(0,1,2)モデル、つまり連続型成長率についてのMA(2)モデルが選ばれる。
参考文献
- 計量経済学（有斐閣）http://www.yuhikaku.co.jp/books/detail/9784641053854
- 基礎からわかる時系列分析（技術評論社）https://gihyo.jp/book/2018/978-4-7741-9646-6
- 経済・ファイナンスデータの計量時系列分析（朝倉書店）https://www.asakura.co.jp/books/isbn/978-4-254-12792-8/