Data Analysis Using Statistical Packages: Time Series Analysis
Sho Kuroda / 黒田翔
2024年7月 (Last update: 2024年7月)
社会科学における時系列分析の初学者に最適な和書は通称「沖本本」と呼ばれる次の一冊． このページは沖本本に依拠している．
- 沖本竜義『経済・ファイナンスデータの計量時系列分析』朝倉書店，2010．
1 Data
quantmod
パッケージを使って1990年から2019年の米国における失業率を月次で取得し
GitHub に CSV でアップロードした． \(T = 12
\times 30 = 360.\)
library(tidyverse)
# install.packages("quantmod")
library(quantmod)
# unemp <- getSymbols(Symbols = "UNRATE", src = "FRED",
# from = "1990-01-01", to = "2019-12-31", auto.assign = FALSE)
# chart_Series(unemp)
# unemp_df <- as.data.frame(time = index(unemp), unemp)
# setwd("c://ws_stat"); write.csv(unemp_df, "unemp_1990_2019.csv")
unemp <- read.csv("https://raw.githubusercontent.com/kurodaecon/bs/main/data/unemp_1990_2019.csv")
unemp_xts <- xts(unemp[, 2], order.by = as.Date(unemp$time)) # 時系列オブジェクトとして定義
plot(unemp_xts)
ur <- unemp$UNRATE # unemployment rate in US
c(summary(ur), SD = sd(ur))## Min. 1st Qu. Median Mean 3rd Qu. Max. SD
## 3.500000 4.600000 5.500000 5.841667 6.800000 10.000000 1.5962312 Basics
2.1 Autocorrelation coefficient / 自己相関係数
自己相関係数 \(\rho_{kt}\) は同一の時系列データにおける異時点間の共分散である自己共分散（autocovariance） \(\gamma_{1t}\) を基準化したもの．
\(k\) 次の自己共分散：
\[ \gamma_{kt} = Cov(y_t, y_{t-k}) = E[(y_t - \mu_t)(y_{t-k} - \mu_{t-k})] , \quad \mu_{t-k} \equiv E(y_{t-k}) \]
\(k\) 次の自己相関係数：
\[ \rho_{kt} = Corr(y_t, y_{t-k}) = \frac{Cov(y_t, y_{t-1})}{\sqrt{Var(y_t) \cdot Var(y_{t-1})}} = \frac{\gamma_{kt}}{\sqrt{\gamma_{0t} \gamma_{0,t-k}}} \]
自己相関係数を \(k\) の関数として見たものは自己相関関数と呼ばれ，これをグラフにしたものはコレログラム（correlogram）と呼ばれる．
定常性やエルゴード性の仮定の下で次のように標本自己相関係数 \(\hat{\rho}_k\) を計算できる．
\[ \hat{\rho}_k = \frac{\hat{\gamma}_k}{\hat{\gamma}_0}, \quad \hat{\gamma}_k \equiv \frac{1}{T} \sum_{t = k+1}^T (y_t - \bar{y})(y_{t-k} - \bar{y}), \quad \bar{y} \equiv \frac{1}{T} \sum_{t=1}^T y_t \]
2.1.1 Calculate according to the definition
標本自己相関係数を計算する．
mu <- mean(ur)
n_ur <- length(ur)
gamma0 <- sum((ur - mu)^2) / n_ur
gamma1 <- sum((ur[-1] - mu) * (ur[-n_ur] - mu)) / n_ur
gamma2 <- sum((ur[-(1:2)] - mu) * (ur[-((n_ur-1):n_ur)] - mu)) / n_ur
gamma3 <- sum((ur[-(1:3)] - mu) * (ur[-((n_ur-2):n_ur)] - mu)) / n_ur
c(gamma0, gamma1, gamma2, gamma3) # 自己共分散## [1] 2.540875 2.521875 2.499946 2.472201c(gamma1, gamma2, gamma3) / gamma0 # 自己相関係数## [1] 0.9925222 0.9838919 0.97297222.1.2 Calculate using
acf function
acf 関数を使えば簡単に計算できる．
acf(ur, lag.max = 5, type = "covariance", plot = FALSE) # 自己共分散##
## Autocovariances of series 'ur', by lag
##
## 0 1 2 3 4 5
## 2.54 2.52 2.50 2.47 2.44 2.40acf(ur, lag.max = 5, type = "correlation", plot = FALSE) # 自己相関係数##
## Autocorrelations of series 'ur', by lag
##
## 0 1 2 3 4 5
## 1.000 0.993 0.984 0.973 0.960 0.944acf(ur, lag.max = 200, type = "correlation") # コレログラム
2.1.3 Test for \(\rho\)
\(H_0: \rho_k = 0\) を検定するためには \(H_0\) のもとでの \(\hat{\rho}_k\) の分布が必要．
\(y_t\) がiid系列の場合には漸近的に \(\hat{\rho}_k \sim N(0, 1/T)\) であることが知られているため，有意水準5%の場合は \(|\hat{\rho}_k| > 1.96 / \sqrt{T}\) を計算することで検定が可能となる．
1.96 / sqrt(n_ur)## [1] 0.1033011上で描いたコレログラムではこの値が青色の破線で表されている．
2.2 Stationarity / 定常性
ADF (Augmented Dickey-Fuller) 検定を行う．
library(tseries)
tseries::adf.test(x = ur)##
## Augmented Dickey-Fuller Test
##
## data: ur
## Dickey-Fuller = -2.5298, Lag order = 7, p-value = 0.3533
## alternative hypothesis: stationarytseries::adf.test(x = na.omit(diff(ur)))##
## Augmented Dickey-Fuller Test
##
## data: na.omit(diff(ur))
## Dickey-Fuller = -3.7087, Lag order = 7, p-value = 0.02385
## alternative hypothesis: stationary原系列 (level data) では \(p > 0.05\) であるため非定常である一方（原系列の図を見れば検定するまでもないが），差分を取ると \(p < 0.05\) であるため定常と言ってよいだろう．
この場合の自己相関係数は以下の通り．
acf(na.omit(diff(ur)), lag.max = 5, type = "correlation", plot = FALSE) # 自己相関係数##
## Autocorrelations of series 'na.omit(diff(ur))', by lag
##
## 0 1 2 3 4 5
## 1.000 0.119 0.242 0.249 0.212 0.256acf(na.omit(diff(ur)), lag.max = 5, type = "correlation") # コレログラム
2.3 White noise / ホワイトノイズ
\(W.N.(\sigma^2)\) は分散が \(\sigma^2\) の white noise を表す． White noise はすべての時点において期待値が 0，分散が一定，かつ自己相関を持たない． すなわち white noise \(\epsilon_t\) は以下を満たす：
\[ E(\epsilon_t) = 0, \quad E(\epsilon_t, \epsilon_{t-k}) = \begin{cases} \sigma^2 \quad \mbox{if } k = 0 \\ 0 \quad \mbox{if } k \ne 0 \end{cases} \]
3 Autoregressive (AR) process / 自己回帰過程
\(p\) 次のAR過程は「AR(p)」と表される：
\[ y_t = c + \sum_{s=1}^{p} \phi_s y_{t-s} + \epsilon_t, \quad \epsilon_t \sim W.N. (\sigma^2) \]
ここでは1次のAR過程 AR(1) を推定する方法を示す．
\[ y_t = c + \phi y_{t-1} + \epsilon_t, \quad \epsilon_t \sim W.N. (\sigma^2) \]
3.1 Estimation using OLS / 最小二乗法
\[ \left( \begin{matrix} \hat{c} \\ \hat{\phi} \end{matrix} \right) = \arg \min_{c, \phi} \sum_{t=2}^T \hat{\epsilon}_t^2 \]
ただし，ここでは初期値 \(y_0\) を無視して推定する．
上の通り差分を取ると定常になるため，差分を取った時系列を \(y_t\) として扱う． 以下同様．
unemp_2 <- unemp %>%
mutate(yt = c(NA, diff(UNRATE)), # y_{t}
yt1 = lag(yt)) # y_{t-1}
summary(lm(formula = yt ~ yt1, data = unemp_2))##
## Call:
## lm(formula = yt ~ yt1, data = unemp_2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.54325 -0.09585 0.00415 0.10415 0.46860
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.004153 0.008086 -0.514 0.6078
## yt1 0.118502 0.052598 2.253 0.0249 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1529 on 356 degrees of freedom
## ( 2 個の観測値が欠損のため削除されました )
## Multiple R-squared: 0.01406, Adjusted R-squared: 0.01129
## F-statistic: 5.076 on 1 and 356 DF, p-value: 0.024873.2 Estimation using ML / 最尤法
対数尤度関数は \(\Theta = (c, \phi, \sigma^2)'\) を用いて次のように表される：
\[ L (\Theta) = \sum_{t=1}^T \log f_{Y_t | Y_{t-1}} (y_t | y_{t-1} ; \Theta) = \sum_{t=1}^T \log \frac{1}{\sqrt{2 \pi \sigma^2}} \exp \left[ \frac{- (y_t - c - \phi y_{t-1})^2}{2 \sigma^2} \right] \] \[ = - \frac{T}{2} \log (2 \pi \sigma^2) - \sum_{t=1}^T \left[ \frac{(y_t - c - \phi y_{t-1})^2}{2 \sigma^2} \right] \]
2番目の等式は，\(\epsilon_t \sim \mbox{iid } N(0, \sigma^2)\) のとき \(y_t | y_{t-1} \sim N(c + \phi y_{t-1}, \sigma^2)\) より導かれる．
これを最大化する \(c, \phi\) は \(\sum_t (y_t - c - \phi y_{t-1})^2\) を最小化する \(c, \phi\) であり，これはOLS推定量に他ならない．
\(\epsilon_t\) の分散の最尤推定量は次のように与えられる：
\[ \hat{\sigma}^2 = \frac{1}{T} \sum_{t=1}^T (y_t - \hat{c} - \hat{\phi} y_{t-1})^2 \]
対数尤度関数を optim
関数を使って最大化してみると，概ね同じ値が得られることが分かる．
ただし，OLSの場合と同様に初期値 \(y_0\)
を無視しており上の尤度関数とは僅かに異なる．
logl <- function (param, Y) { # 対数尤度関数
c_hat <- param[1] # estimate of c
phi_hat <- param[2] # estimate of phi
sigma2_hat <- param[3] # estimate of sigma^2
yt <- Y$yt[-1:-2] # y_{t-1} の1-2行目は NA なのでデータから除外
yt1 <- Y$yt1[-1:-2]
Tmax <- nrow(Y) - 2
# 対数尤度
- Tmax/2*log(2 * pi * sigma2_hat) - sum((yt - c_hat - phi_hat * yt1)^2 / (2 * sigma2_hat))
}
optim(par = c(0, 0.1, 0.02), fn = logl, Y = unemp_2, control = list(fnscale = -1))## $par
## [1] -0.004157148 0.118496998 0.023249788
##
## $value
## [1] 165.3249
##
## $counts
## function gradient
## 124 NA
##
## $convergence
## [1] 0
##
## $message
## NULL3.3 Estimation using
arima function
AR と後述する MA を統合した ARIMA (AR integrated MA) model を推定する
arima 関数を使うのが便利．
order 引数は，AR order, degree of differencing, MA order
の順で指定する．
arima(x = unemp_2$yt, order = c(1, 0, 0))##
## Call:
## arima(x = unemp_2$yt, order = c(1, 0, 0))
##
## Coefficients:
## ar1 intercept
## 0.1183 -0.0050
## s.e. 0.0524 0.0091
##
## sigma^2 estimated as 0.02321: log likelihood = 166.09, aic = -326.183.4 Estimation using
forecast::Arima function
差分系列に AR(1) を適用することは原系列に ARIMA(1, 1, 0) を適用することに等しい．
forecast パッケージの Arima 関数を使えば
include.drift
引数を指定して定数項を含めた推定ができる．
# install.packages("forecast")
library(forecast)
forecast::Arima(y = ur, order = c(1, 1, 0), include.drift = TRUE)## Series: ur
## ARIMA(1,1,0) with drift
##
## Coefficients:
## ar1 drift
## 0.1183 -0.0050
## s.e. 0.0524 0.0091
##
## sigma^2 = 0.02334: log likelihood = 166.09
## AIC=-326.18 AICc=-326.11 BIC=-314.533.5 Checking the residuals after fitting an AR(1) model
AR(1) を適用した後の残差について自己相関係数を計算すると以下のようになるため，AR(1) ではこの時系列データを十分に記述できていないと判断される．
forecast::checkresiduals(arima(x = unemp_2$yt, order = c(1, 0, 0)))
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,0,0) with non-zero mean
## Q* = 93.458, df = 9, p-value = 3.331e-16
##
## Model df: 1. Total lags used: 10Ljung-Box 検定の帰無仮説は「データは独立に分布している」なので，これが棄却されることは残差に自己相関（系列相関）があることを意味している．
4 Moving average (MA) process / 移動平均過程
\(q\) 次のMA過程は「MA(q)」と表される：
\[ y_t = \mu + \epsilon_t + \sum_{s=1}^{q} \theta_s \epsilon_{t-s}, \quad \epsilon_t \sim W.N. (\sigma^2) \]
1次のMA過程 MA(1) は次の通り：
\[ y_t = \mu + \epsilon_t + \theta \epsilon_{t-1}, \quad \epsilon_t \sim W.N. (\sigma^2) \]
4.1 Estimation using ML / 最尤法
対数尤度関数は \(\Theta = (\mu, \theta, \sigma^2)'\) を用いて次のように表される：
\[ L (\Theta) = - \frac{T}{2} \log (2 \pi \sigma^2) - \sum_{t=1}^T \left[ \frac{\epsilon_t^2}{2 \sigma^2} \right], \] \[ \epsilon_t = y_t - \mu - \theta \epsilon_{t-1}, \quad \epsilon_1 = y_1 - \mu \]
対数尤度関数を optim 関数を使って最大化する．
logl_ma <- function (param, Y) { # 対数尤度関数
mu_hat <- param[1] # estimate of mu
theta_hat <- param[2] # estimate of theta
sigma2_hat <- param[3] # estimate of sigma^2
yt <- Y$yt[-1] # y_{t} の1行目は NA なのでデータから除外
Tmax <- nrow(Y) - 1
# epsilon を計算
epsilon_hat <- numeric(length = Tmax)
epsilon_hat[1] <- yt[1] - mu_hat # t = 1
for (i in 2:Tmax) { # t >= 2
epsilon_hat[i] <- yt[i] - mu_hat - theta_hat * epsilon_hat[i-1]
}
# 対数尤度
- Tmax/2*log(2 * pi * sigma2_hat) - sum(epsilon_hat^2 / (2 * sigma2_hat))
}
optim(par = c(0, 0.1, 0.02), fn = logl_ma, Y = unemp_2, control = list(fnscale = -1))## $par
## [1] -0.005027544 0.083652019 0.023310866
##
## $value
## [1] 165.3193
##
## $counts
## function gradient
## 114 NA
##
## $convergence
## [1] 0
##
## $message
## NULL4.2 Estimation using
arima function
差分系列に対する ARIMA(0, 0, 1) または原系列に対する ARIMA(0, 1, 1) で推定できる．
arima(x = unemp_2$yt, order = c(0, 0, 1))##
## Call:
## arima(x = unemp_2$yt, order = c(0, 0, 1))
##
## Coefficients:
## ma1 intercept
## 0.0835 -0.0050
## s.e. 0.0448 0.0087
##
## sigma^2 estimated as 0.02331: log likelihood = 165.32, aic = -324.63# forecast::Arima(y = ur, order = c(0, 1, 1), include.drift = TRUE) # same as above 4.3 Model selection
AIC などの情報量規準を使ってモデル選択をすることもできる．
forecast::auto.arima(y = ur, ic = "aic", trace = TRUE)##
## Fitting models using approximations to speed things up...
##
## ARIMA(2,1,2) with drift : -374.9525
## ARIMA(0,1,0) with drift : -319.3597
## ARIMA(1,1,0) with drift : -321.8254
## ARIMA(0,1,1) with drift : -320.8906
## ARIMA(0,1,0) : -320.9765
## ARIMA(1,1,2) with drift : -371.9543
## ARIMA(2,1,1) with drift : -368.25
## ARIMA(3,1,2) with drift : -371.6955
## ARIMA(2,1,3) with drift : -372.9533
## ARIMA(1,1,1) with drift : -359.9761
## ARIMA(1,1,3) with drift : -372.3223
## ARIMA(3,1,1) with drift : -373.6062
## ARIMA(3,1,3) with drift : -371.5262
## ARIMA(2,1,2) : -376.9169
## ARIMA(1,1,2) : -373.9322
## ARIMA(2,1,1) : -370.2166
## ARIMA(3,1,2) : -373.5588
## ARIMA(2,1,3) : -374.9183
## ARIMA(1,1,1) : -361.943
## ARIMA(1,1,3) : -374.304
## ARIMA(3,1,1) : -375.4581
## ARIMA(3,1,3) : Inf
##
## Now re-fitting the best model(s) without approximations...
##
## ARIMA(2,1,2) : -382.106
##
## Best model: ARIMA(2,1,2)## Series: ur
## ARIMA(2,1,2)
##
## Coefficients:
## ar1 ar2 ma1 ma2
## 1.4601 -0.5242 -1.4994 0.6695
## s.e. 0.1519 0.1565 0.1393 0.1377
##
## sigma^2 = 0.01983: log likelihood = 196.05
## AIC=-382.11 AICc=-381.94 BIC=-362.69ARIMA(2, 1, 2) が選択される：
\[ \Delta y_{t} = \phi_1 \Delta y_{t-1} + \phi \Delta y_{t-2} + \epsilon_{t} + \theta_1 \epsilon_{t-1} + \theta_2 \epsilon_{t-2} \]
4.4 Incorporating seasonality / 季節性
データが12か月周期であることを宣言してARIMAモデルで季節性を考慮する（SARIMAモデルとして推定する）こともできる．
ur_ts <- ts(data = ur, start = c(1990, 1), frequency = 12)
# arima(ur_ts, order = c(1, 1, 1), seasonal = list(order = c(0, 0, 1)))
auto.arima(ur_ts)## Series: ur_ts
## ARIMA(1,1,3)(0,0,2)[12]
##
## Coefficients:
## ar1 ma1 ma2 ma3 sma1 sma2
## 0.9524 -1.0094 0.1503 0.0905 -0.2964 -0.2365
## s.e. 0.0224 0.0563 0.0731 0.0565 0.0516 0.0467
##
## sigma^2 = 0.01793: log likelihood = 213.62
## AIC=-413.25 AICc=-412.93 BIC=-386.06ARIMA(1, 1, 3)(0, 0, 2)[12] が選択される：
\[ \Delta y_{t} = \phi_1 \Delta y_{t-1} + \epsilon_{t} + \theta_1 \epsilon_{t-1} + \theta_2 \epsilon_{t-2} + \theta_3 \epsilon_{t-3} + \vartheta_1 \epsilon_{t-12} + \vartheta_2 \epsilon_{t-24} \]
4.5 Forecasting / 予測
次のように予測することもできる．
sarima_forecast <- forecast::forecast(auto.arima(ur_ts), h = 12)
# summary(sarima_forecast)
ggplot2::autoplot(sarima_forecast, ylab = "Unemployment rate [%]")