経営のための統計学II

第4回 時系列分析2

自己回帰推定の診断

自己回帰式は、変数の原系列を特定のモデルで変換し、ホワイト・ノイズに近い性質をもつ残差を導くという役割を果たす。

• たとえば、GDPの原系列（対数系列）は時間に関して独立ではなく互いに相関しているが、真のモデルがAR(1)ならAR(1)式から得た残差はホワイト・ノイズの性質をもつと考えられる。
• 自己回帰式は、標準的な条件を満たさない時系列変数を、ホワイト・ノイズに変換している。
  • この変換をフィルター（濾過）と呼ぶこともある。

実際の時系列分析においては真のフィルターは未知だから、自己回帰式の推定を行い、残差がホワイト・ノイズの条件を満たすかどうか検討する。

• これを診断検定（diagnostic test）という。
• ここでは、以下の方法を説明する。
  1. 残差のプロット
  2. 残差の標本AC関数
  3. ふろしき検定

残差のプロット

時系列の平均や誤差の分散が観測期間の途中で変化することを構造変化（structural change）という。

• 残差のプロットを検討することにより、構造変化と異常値（特定の観測点における異常に大きい値あるいは小さな値）を検出することができる。

鉱工業生産指数（IIP）の月次成長率についてAR(2)を推定した結果は次のようになる。

# IIPデータの読み込み
iip <- read.csv("iip.csv")
head(iip)

##   yyyymm  iip
## 1 197801 63.0
## 2 197802 62.4
## 3 197803 63.9
## 4 197804 63.9
## 5 197805 64.1
## 6 197806 64.4

# データをts型に変換
x <- ts(data = iip$iip, start = c(1978, 1), frequency = 12)

# 鉱工業生産指数の連続型成長率の標本AC関数と2シグマ
r <- diff(log(x)) #連続型成長率
acf(r, lag.max = 20, ci.type = "ma")

# AR(2)の推定
ar2 <- arima(r, order = c(2, 0, 0))

# 残差のプロット
par(family = "HiraginoSans-W3") #日本語が文字化けしないようにおまじない
ts.plot(ar2$residuals, gpars=list(xlab="年月", ylab="残差", lty=3, col=3))

• 平均の構造変化を見出すことはできない。
• 2009年と2011年に残差が大きな値を示している。

２シグマを超える残差を調べる。

# ２シグマ
twosigma <- 2*sqrt(ar2$sigma2)
twosigma

## [1] 0.03601337

# ２シグマを超える残差（異常値）のインデックスを抽出
i.outlier <- abs(ar2$residuals) > twosigma
i.outlier

##        Jan   Feb   Mar   Apr   May   Jun   Jul   Aug   Sep   Oct   Nov   Dec
## 1978       FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 1979 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 1980 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 1981 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 1982 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 1983 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 1984 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 1985 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 1986 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 1987 FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE
## 1988 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 1989 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 1990 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 1991 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 1992 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE
## 1993 FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 1994 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 1995 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 1996 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 1997 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 1998 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 1999 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 2000 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 2001  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 2002 FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 2003 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 2004 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 2005 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 2006 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 2007 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 2008 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE  TRUE
## 2009  TRUE  TRUE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 2010 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 2011 FALSE FALSE  TRUE FALSE  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE
## 2012 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 2013 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 2014 FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 2015 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 2016 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 2017 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE

# 年月ラベルを作成
date <- seq(as.Date("1978-02-01"), as.Date("2017-12-01"), by = "month")

# 年月と残差のデータフレーム
zansa <- data.frame(date = format(date, "%Y-%m"), residuals = ar2$residuals)
head(zansa)

##      date    residuals
## 1 1978-02 -0.010566482
## 2 1978-03  0.022383563
## 3 1978-04  0.000470380
## 4 1978-05 -0.001010678
## 5 1978-06  0.003734482
## 6 1978-07  0.004847762

#  ２シグマを超える年月と残差の値
zansa[i.outlier, ]

##        date   residuals
## 113 1987-06  0.03732241
## 176 1992-09  0.03943973
## 182 1993-03 -0.04436986
## 276 2001-01 -0.04475336
## 292 2002-05  0.04050402
## 370 2008-11 -0.07178681
## 371 2008-12 -0.08403957
## 372 2009-01 -0.08506241
## 373 2009-02 -0.07946025
## 375 2009-04  0.05445839
## 398 2011-03 -0.18140191
## 400 2011-05  0.08941582
## 401 2011-06  0.03767297
## 435 2014-04 -0.04458957

• 2008年末から2009年はじめは世界金融危機の影響と考えられる。
• 2011年3月以降は大震災の影響と考えられる。

残差の標本AC関数

残差系列から計算した標本AC関数について統計的な検定を行う。

• ホワイト・ノイズのAC関数はすべての次数で０に近い値を示すので、これが時系列構造を検出する基準になる。
• 何らかの時系列構造が残差に残っていれば、AC関数はスパイク（尖り）を示し、０とは有意に異なる値を示すと考えられる。
• 残差がホワイト・ノイズであるという帰無仮説のもとでは、\(\rho(k) = 0\) \((k=1,2,\ldots)\)で、ACの分散は \[ V(AC) = \frac{1}{n} \]
• 有意水準5%の検定の境界値は\(\pm \frac{1.96}{\sqrt{n}} \approx \frac{2}{\sqrt{n}}\)で、この境界を超えるACは有意であると判定される。
  • 標本AC係数の\(t\)値は、標本AC係数に\(\sqrt{n}\)をかけて求められる。
  • \(t\)値による検定は、ACプロットにおいてACが上下の有意性境界を超えるか否かの判定と同じ。

上記の四半期GDPの対数系列のAR(1)から計算された残差のACプロットは次のようになる。

# 残差のACプロット
acf(ar2$residuals, lag.max = 20)

• 有意なAC値はまったく見られない。

自己回帰式の選択においては、原系列に関する標本AC関数のみならず標本PAC関数が分析に用いられたが、診断検定では残差に関する標本AC関数飲みが使われる。

• これは残差系列に何らかの相関が残っていれば、残差の標本AC関数に反映され、有意なACが見つかるからである。
• 残差に関しては、標本AC関数はPAC関数と非常によく似た変動を示すことが知られている。

# 残差のPACプロット
pacf(ar2$residuals, lag.max = 20)

ふろしき検定

診断検定では、推定された時系列過程の残差から計算された標本AC係数の平方和として計算されるふろしき統計量 が重要である。

• 平方和を取る範囲は、検定に先立って例えば\(k\)と決めておく。
• ふろしき統計量は、\(k\)次までのAC係数の推定値を使って次のように定義される。 \[ Q(k) = n \sum_{i=1}^{k}AC(i)^2 \quad \mbox{または} \quad Q^{*}(k) = n(n+2)\sum_{i=1}^{k}\frac{1}{n-i}AC(i)^{2} \]
  • \(Q\)はBox-Pierce（ボックス=ピアース）検定、\(Q^{*}\)はLjung-Box（リュン=ボックス）検定と呼ばれる。
  • 実際の計算では両者に大きな差異は見出せない。
• \(Q\)は「すべての母自己相関係数が０」という帰無仮説のもとで漸近的にカイ２乗分布に従う。
  • \(n\)が十分大きい場合にカイ２乗分布に従うと考えれば良い。
• 自由度は\(k\)から自己回帰の次数を引いた値である。
  • 定数項がモデルに含まれていても自由度には影響しない。
  • 残差ではなく、原系列に関して求めた\(Q\)統計量では、係数を推定する必要がないので自由度は\(k\)となる。
• ふろしき検定では対立仮説は特に指定されていない。
  • 何らかの原因（異常値、不均一分散、構造変化など）により、あるACが大きな値をもつと、検定統計量の値が大きくなり、帰無仮説は棄却される可能性が高くなる。
  • この検定では帰無仮説のモデルから何らかの乖離を含むモデルを対立仮説として考える。
  • 乖離の内容は自由に考えて良いから、「ふろしき」検定と呼ばれる。

Rでは関数Box.testを使ってふろしき検定を行うことができる。

# ふろしき検定
k <- 8 #次数
Box.test(ar2$residuals, k, type = "Box-Pierce", fitdf = 2)

## 
##  Box-Pierce test
## 
## data:  ar2$residuals
## X-squared = 4.2202, df = 6, p-value = 0.6469

Box.test(ar2$residuals, k, type = "Ljung-Box", fitdf = 2)

## 
##  Box-Ljung test
## 
## data:  ar2$residuals
## X-squared = 4.2899, df = 6, p-value = 0.6375

• 自己回帰の次数は2だからfitdf = 2で自由度を\(k\)から2引いている。
• Box-Pierce検定とLjung-Box検定について、両者の統計量の値と\(p\)値に大きな差はない。

残差を回帰診断するため、\(Q^{*}\)統計量とその\(p\)値を８ラグ分まで示す。

# k=1から8までのLjung-Box検定統計量とp値
Q_LB <- numeric(8)
p_LB <- numeric(8)
for (k in 1:8) {
  LB <- Box.test(ar2$residuals, k, type = "Ljung-Box", fitdf = 2)
  Q_LB[k] <- LB$statistic
  p_LB[k] <- LB$p.value
}

## Warning in pchisq(STATISTIC, lag - fitdf): 計算結果が NaN になりました

data.frame(k = 1:8, Q_LB, p_LB)

##   k       Q_LB      p_LB
## 1 1 0.02148020       NaN
## 2 2 0.02159477 0.0000000
## 3 3 1.21303727 0.2707314
## 4 4 1.23327120 0.5397573
## 5 5 1.49840263 0.6826391
## 6 6 2.08325407 0.7204495
## 7 7 3.28200962 0.6565971
## 8 8 4.28994128 0.6375034

• \(k=3\)で自由度が\(k-2=1\)になり、検定統計量として意味をもつ。
• \(p\)値は3以上のすべての次数について5%より大きい。
• 一般的に、残差に関するPACとACは類似するが、AR(2)の残差についても両者にほとんど差異がない。
• したがって、残差にはARおよび次節で説明するMA過程は残っていないと判断できる。

移動平均（MA）法

時系列分析においては、観測される変数\(X_{t}\)についての自己回帰のみでなく、ノイズの移動平均が分析に利用される。

• 移動平均（moving average：MA）過程という。
• 最も単純なMA過程は次のようになる。 \[ \begin{align} X_{t} &= \mu + u_{t} \\ u_{t} &= \epsilon_{t} + \theta \epsilon_{t-1} \end{align} \]
  • 1次のMA過程と呼ばれMA(1)と記される。
  • \(\epsilon_{t}\)は自己回帰式の\(\epsilon_{t}\)と同様、ホワイト・ノイズの条件を満たす確率変数である。
  • ノイズをまとめて\(u_{t}\)と表記する。
• MA過程では観測できないノイズに関して式を想定しているため、モデルの直感的な理解が難しい。
• しかし、このような表現から\(u_{t}\)に関する性質を導くことができる。
• まず、MA(1)を考える。
  • \(u_{t}\)と1期前の\(u_{t-1} = \epsilon_{t-1} + \theta \epsilon_{t-2}\)は共通な項\(\epsilon_{t-1}\)を含むから、隣のノイズ間の共分散は０にはならない。
  • ２つのノイズの時間差が２以上あれば、共通な\(\epsilon\)がなくなるので共分散は０になる。
  • 上記の式からわかるように、\(u_{t}\)の相関は、\(X_{t}\)の相関と変わらない。
  • したがって、MA(1)では、観測できる変数\(X_{t}\)と\(X_{t-1}\)は相関をもつが、ラグが2以上あれば相関がなくなる。
• MA(\(k\))では\(u_{t}\)は次のように定義される。 \[ u_{t} = \epsilon_{t} + \theta_{1}\epsilon_{t-1} + \cdots + \theta_{k}\epsilon_{t-k} \]
  • この場合は、\(k\)時点離れている\(u_{t-k}\)まで共分散が０にならないが、\(k+1\)時点以上離れれば共分散は０になる。
  • このように、\(u_{t}\)の有限次のラグが相関するのが、MA過程の特色である。
  • 変数\(X_{t}\)についても同じで、\(MA(k)\)により、\(k\)次の時間差まで相関がある時系列変数を表現することができる。
• このように考えると、AR過程は無限次のMA過程に変換できたから、観測できる変数\(X_{t}\)は、無限次のラグまで相関をもつことになる。
• AR過程が最小２乗法で推定できるのと比べると、MA過程の最尤推定は困難であるが、Rなどの分析ソフトを利用すれば低次のMA過程は字数を指定するだけで推定できる。
• MA過程では、経済変数の変動が、観測できない確率変数\(\epsilon_{t}\)の変動で説明されるため、経済分析において好まれて利用されることはない。
  • 言葉は悪いが、\(\epsilon_{t}\)はいわばゴミのようなもので、経済変数の変動をゴミの和によって説明しようとするのは、経済的な意味合いに欠けると考えられる。
  • このため、MA過程が経済分析で単独で利用されることはほとんどない。

MA過程のAR表現

前回、AR過程が定常性の条件を満たせば時系列変数をノイズ\(\epsilon_{t}\)だけの1次式として表現できることを示した。 \[ \begin{align} X_t &= \mu + \phi X_{t-1} + \epsilon_{t} \\ &= \epsilon_t + \phi \epsilon_{t-1} + \phi^2 \epsilon_{t-2} + \cdots \end{align} \] 同様に、ノイズの1次式であるMA過程を、時系列変数\(X_{t}\)だけの1次式として表現することができる。

• これを、MA過程の反転という。
• MA過程の反転とは、MA過程のAR表現を意味する。
• 反転可能性条件と呼ばれるAR過程の定常性に類した条件が満たされるならば、MA過程のAR表現は一意になる。
• 例としてMA(1)を考える。 \[ X_{t} = \mu + \epsilon_{t} + \theta \epsilon_{t-1} \]
  • ノイズは \[ \epsilon_{t} = X_{t} - \mu - \theta \epsilon_{t-1} \]
  • \(t-1\)に当てはまめたノイズ\(\epsilon_{t-1} = X_{t-1} - \mu - \theta \epsilon_{t-2}\)をもとの式に代入すると \[ \begin{align} X_{t} &= \mu + \epsilon_{t} + \theta (X_{t-1} - \mu - \theta \epsilon_{t-2}) \\ &= (1 - \theta)\mu + \theta X_{t-1} + \epsilon_{t} - \theta^{2}\epsilon_{t-2} \end{align} \]
  • 代入を繰り返すと \[ \begin{align} X_{t} &= (1 - \theta + \theta^{2} - \cdots )\mu + \theta X_{t-1} - \theta^{2} X_{t-2} + \theta^{3} X_{t-3} + \cdots + \epsilon_{t} \\ &= c' + \theta X_{t-1} - \theta^{2} X_{t-2} + \theta^{3} X_{t-3} + \cdots + \epsilon_{t} \end{align} \]
    • ここで、\(c' = \frac{\mu}{1 + \theta}\)と定義する。
    • \(|\theta|<1\)と仮定すると、\(\theta^{m}\epsilon_{t-m}\)は0、級数\((1 - \theta + \theta^{2} - \cdots)\)は\(\frac{1}{1+\theta}\)で近似できる。
  • \(X_{t}\)のラグは無限に続くので、MA過程はAR(\(\infty\))に書き換えられている。
  • MA過程は定義によりノイズの1次式であるから、常に定常である。
  • 条件\(|\theta|<1\)は、定常性の条件ではなく、MA過程をAR仮定で表現するために役立っている。

推定例

鉱工業生産指数の月次成長率に関する4次の移動平均式MA(4)の推定結果を示す。

• AR式の推定で用いたarima関数を使って簡単にMAモデルを最尤法により推定できる。
  • オプションのorderの最後の数値にMAの次数を指定する。

# MA(4)の推定
ma4 <- arima(r, order = c(0,0,4))
ma4

## 
## Call:
## arima(x = r, order = c(0, 0, 4))
## 
## Coefficients:
##           ma1     ma2     ma3     ma4  intercept
##       -0.0127  0.1325  0.0405  0.0104     0.0011
## s.e.   0.0457  0.0457  0.0447  0.0429     0.0010
## 
## sigma^2 estimated as 0.0003236:  log likelihood = 1244.93,  aic = -2477.85

最高次の項ma4が有意でないので、次数を下げて推定を繰り返す。

# MA(3)の推定
ma3 <- arima(r, order = c(0,0,3))
ma3

## 
## Call:
## arima(x = r, order = c(0, 0, 3))
## 
## Coefficients:
##           ma1     ma2     ma3  intercept
##       -0.0125  0.1316  0.0401     0.0011
## s.e.   0.0456  0.0450  0.0447     0.0010
## 
## sigma^2 estimated as 0.0003236:  log likelihood = 1244.9,  aic = -2479.8

# MA(2)の推定
ma2 <- arima(r, order = c(0,0,2))
ma2

## 
## Call:
## arima(x = r, order = c(0, 0, 2))
## 
## Coefficients:
##           ma1     ma2  intercept
##       -0.0172  0.1350     0.0011
## s.e.   0.0455  0.0451     0.0009
## 
## sigma^2 estimated as 0.0003242:  log likelihood = 1244.5,  aic = -2480.99

# MA(1)の推定
ma1 <- arima(r, order = c(0,0,1))
ma1

## 
## Call:
## arima(x = r, order = c(0, 0, 1))
## 
## Coefficients:
##           ma1  intercept
##       -0.0053     0.0011
## s.e.   0.0406     0.0008
## 
## sigma^2 estimated as 0.0003302:  log likelihood = 1240.13,  aic = -2474.25

• 対数尤度log likelihoodはMA(4)が最も大きい。

項数が増えることに対するペナルティを入れたAICを比較する。

# AICの比較
AIC(ma1, ma2, ma3, ma4)

##     df       AIC
## ma1  3 -2474.251
## ma2  4 -2480.994
## ma3  5 -2479.795
## ma4  6 -2477.854

• 最小AIC規準によれば、MA(2)が選ばれる。

移動平均過程の診断

推定結果を用いた診断は、自己回帰の場合と変わらない。

1. 残差のプロット
2. 残差の標本AC関数
3. ふろしき検定

まず、フィット値と残差のプロットを見てみる。

# 成長率とフィット値のプロット
par(family = "HiraginoSans-W3") #日本語が文字化けしないようにおまじない
ts.plot(r, r - ma2$residuals, gpars=list(xlab="年月", ylab="成長率", lty=c(1:2), col=c(1:2))) #成長率とフィット値

# 残差のプロット
par(family = "HiraginoSans-W3") #日本語が文字化けしないようにおまじない
ts.plot(ma2$residuals, gpars=list(xlab="年月", ylab="残差", lty=3, col=3))

• フィット値を見てもMA過程の特色は見えない

残差が２シグマを超える場合を調べる。

# ２シグマ
twosigma <- 2*sqrt(ma2$sigma2)
twosigma

## [1] 0.03601049

# ２シグマを超える残差（異常値）のインデックスを抽出
i.outlier <- abs(ma2$residuals) > twosigma

# 年月と残差のデータフレーム
zansa <- data.frame(date = format(date, "%Y-%m"), residuals = ma2$residuals)

#  ２シグマを超える年月と残差の値
zansa[i.outlier, ]

##        date   residuals
## 25  1980-02  0.03632681
## 113 1987-06  0.03695657
## 176 1992-09  0.03855896
## 182 1993-03 -0.04462793
## 276 2001-01 -0.04505106
## 292 2002-05  0.04024675
## 370 2008-11 -0.07199758
## 371 2008-12 -0.08545244
## 372 2009-01 -0.08569674
## 373 2009-02 -0.08049433
## 375 2009-04  0.05350490
## 398 2011-03 -0.18116492
## 400 2011-05  0.08980717
## 401 2011-06  0.03930056
## 435 2014-04 -0.04452816

• MA(2)の残差は1980年2月も２シグマを超えており、AR(2)と比べると精度は落ちると考えられる。

次に、残差の標本AC関数と標本PAC関数を見てみる。

# 残差のACプロット
acf(ma2$residuals, lag.max = 20)

# 残差のPACプロット
pacf(ma2$residuals, lag.max = 20)

• ２シグマを超えているAC係数はなく、自己相関は取り除かれていることがわかる。

最後に、８次までのふろしき検定の結果を見る。

# k=1から8までのLjung-Box検定統計量とp値
Q_LB <- numeric(8)
p_LB <- numeric(8)
for (k in 1:8) {
  LB <- Box.test(ar2$residuals, k, type = "Ljung-Box", fitdf = 2)
  Q_LB[k] <- LB$statistic
  p_LB[k] <- LB$p.value
}

## Warning in pchisq(STATISTIC, lag - fitdf): 計算結果が NaN になりました

data.frame(k = 1:8, Q_LB, p_LB)

##   k       Q_LB      p_LB
## 1 1 0.02148020       NaN
## 2 2 0.02159477 0.0000000
## 3 3 1.21303727 0.2707314
## 4 4 1.23327120 0.5397573
## 5 5 1.49840263 0.6826391
## 6 6 2.08325407 0.7204495
## 7 7 3.28200962 0.6565971
## 8 8 4.28994128 0.6375034

• AR(2)と同様、\(k=3\)で自由度が\(k-2=1\)になり、検定統計量として意味をもつ。
• \(p\)値は3以上のすべての次数について5%より大きい。
• したがって、残差系列に自己相関は残っていないと判断できる。

自己回帰移動平均（ARMA）法

AR過程において誤差項がMAであるモデルを自己回帰移動平均（auto-regressive moving average：ARMA）という。

• 経済分析においては系列の変動を説明する主要な部分は自己回帰であり、移動平均は付属部分である。

最も簡単な例として、１次のARおよび１次のMA部分をもつARMA過程は次のようになる。 \[ X_{t} + \phi X_{t-1} = \mu + \epsilon_{t} + \theta \epsilon_{t-1} \]

• 次数を明示するために、この系列はARMA(1,1)と表記される。

ARMA過程では、今期の変数値（\(X_{t}\)）は、過去の値（\(X_{t-1}, X_{t-2}, \ldots\)）の１次式と今期以前の誤差項（\(\epsilon_{t-1}, \epsilon_{t-2}, \ldots\)）の１次式として表現される。

• AR過程およびMA過程を特殊ケースとして含んでいる。
• 応用では、ARのラグ次数がMAの次数より長くなるのが自然である。

AR過程がMA化でき、MA過程がAR化できるように、ARMA過程をMA過程やAR過程に変換することもできる。

1. ARMA(1,1)のMA化
   • ARMA(1,1)の左辺の\(X_{t-1}\)を \[ X_{t-1} = - \phi X_{t-2} + \mu + \epsilon_{t-1} + \theta \epsilon_{t-2} \] の右辺に置き換える。
   • 新たに出てきた\(X_{t-2}\)もさらに同様に置き換え、この手続きを反復する。
   • このMA過程はエラーショック表現と呼ばれる。
2. ARMA(1,1)のAR化
   • ARMA(1,1)の右辺の\(\epsilon_{t-1}\)を \[ \epsilon_{t-1} = X_{t-1} + \phi X_{t-2} - \mu + - \theta \epsilon_{t-2} \] の右辺に置き換える。
   • 新たに出てきた\(\epsilon_{t-2}\)もさらに同様に置き換え、この手続きを反復する。

こうして求められた時系列過程は無限の項を含むから、MA(\(\infty\))やAR(\(\infty\))のように次数を明記して示すことが適切である。

• ARMA過程のモデルを選ぶとは、時系列変数の変動を、MA(\(\infty\))やAR(\(\infty\))のような複雑なモデルではなく、可能な限り項数の少ない簡潔なモデルで表現することである。

一般のARMA過程は次のようになる。 \[ X_{t} + \phi_{1} X_{t-1} + \phi_{2} X_{t-2} + \cdots + \phi_{p} X_{t-p} = \mu + \epsilon_{t} + \theta_{1} \epsilon_{t-1} + \cdots + \theta_{q} \epsilon_{t-q} \]

• 次数を明示してARMA(\(p,q\))と表記する。

ARMA(\(p,q\))を過去のノイズの１次式で表現するためには、AR部分の係数を用いた多項式 \[ 1 + \phi_{1}x^{1} + \phi_{2}x^{2} + \cdots + \phi_{p}x^{p} = 0 \] の根の絶対値が１より大でないといけない。

• これを定常性の条件という。
• 多項式は \[ x^{p} + \phi_{1}x^{p-1} + \phi_{2}x^{p-2} + \cdots + \phi_{p} = 0 \] となっていることもあるが、この場合は根の絶対値が１より小であることが定常性の条件となる。

一方、ARMA(\(p,q\))の自己回帰への書き直しが可能であるためには、MA部分の係数を用いた多項式 \[ 1 + \theta_{1}x^{1} + \theta_{2}x^{2} + \cdots + \theta_{q}x^{q} = 0 \] の根の絶対値が１より大でないといけない。

• これを反転可能性条件という。
• この条件も、多項式が \[ x^{q} + \theta_{1}x^{q-1} + \theta_{2}x^{q-2} + \cdots + \theta_{q} = 0 \] となっていれば、根の絶対値が１より小という条件にかわる。

推定例と診断

AR(2,2)の推定結果を示す。

# ARMA(2,2)の推定
arma22 <- arima(r, order = c(2, 0, 2))
arma22

## 
## Call:
## arima(x = r, order = c(2, 0, 2))
## 
## Coefficients:
##          ar1     ar2      ma1     ma2  intercept
##       0.3308  0.0897  -0.3437  0.0485     0.0011
## s.e.  0.3190  0.2701   0.3200  0.2685     0.0010
## 
## sigma^2 estimated as 0.0003234:  log likelihood = 1245.07,  aic = -2478.14

• 係数はすべて有意ではない。

AR部分の多項式\(1+\phi_{1}x + \phi_{2}x^{2} = 0\)の根を求める。

• 関数polyrootを用いる。
• 例えば、polyroot(c(1, 0.5, 0.6))で\(1 + 0.5x + 0.6x^{2}\)の根を求めることができる。

# AR多項式の根
root <- polyroot(c(1, arma22$coef[1], arma22$coef[2]))
root

## [1] -1.84332+2.783228i -1.84332-2.783228i

round(abs(root), digits = 3) #根の絶対値

## [1] 3.338 3.338

• 根の絶対値はどちらも1より大きいので、定常性の条件を満たす。

同様にMA部分の多項式\(1+\theta_{1}x + \theta_{2}x^{2} = 0\)の根を求め、反転可能性条件を確認する。

# MA多項式の根
round(abs(polyroot(c(1, arma22$coef[3], arma22$coef[4]))), digits = 3) #根の絶対値

## [1] 4.543 4.543

• 根の絶対値はどちらも1より大きいので、反転可能性条件を満たす。

ARMA(2,2)よりも少ない次数を推定し、AR(2)とMA(2)も含めてAICを比較する。

# ARMA(2,1)の推定
arma21 <- arima(r, order = c(2, 0, 1))
arma21

## 
## Call:
## arima(x = r, order = c(2, 0, 1))
## 
## Coefficients:
##          ar1     ar2      ma1  intercept
##       0.3166  0.1371  -0.3294     0.0011
## s.e.  0.2667  0.0467   0.2685     0.0010
## 
## sigma^2 estimated as 0.0003234:  log likelihood = 1245.05,  aic = -2480.11

# ARMA(1,2)の推定
arma12 <- arima(r, order = c(1, 0, 2))
arma12

## 
## Call:
## arima(x = r, order = c(1, 0, 2))
## 
## Coefficients:
##          ar1      ma1     ma2  intercept
##       0.3821  -0.3941  0.1348     0.0011
## s.e.  0.3530   0.3512  0.0472     0.0010
## 
## sigma^2 estimated as 0.0003235:  log likelihood = 1245.01,  aic = -2480.02

# ARMA(1,1)の推定
arma11 <- arima(r, order = c(1, 0, 1))
arma11

## 
## Call:
## arima(x = r, order = c(1, 0, 1))
## 
## Coefficients:
##           ar1     ma1  intercept
##       -0.8553  0.8181     0.0011
## s.e.   0.1080  0.1177     0.0008
## 
## sigma^2 estimated as 0.0003285:  log likelihood = 1241.35,  aic = -2474.7

# AR(2)の推定
ar2 <- arima(r, order = c(2, 0, 0))

# MA(2)の推定
ma2 <- arima(r, order = c(0, 0, 2))

# AICの比較
AIC(ar2, ma2, arma11, arma12, arma21, arma22)

##        df       AIC
## ar2     4 -2480.918
## ma2     4 -2480.994
## arma11  4 -2474.696
## arma12  5 -2480.023
## arma21  5 -2480.107
## arma22  6 -2478.135

• 最小AIC規準では、MA(2)が選ばれる。

MA(2)の定常性の条件と反転可能性条件を確認する。

# MA(2)のMA特性式の根
round(abs(polyroot(c(1, ma2$coef[1], ma2$coef[2]))), digits = 3) #根の絶対値

## [1] 2.722 2.722

• すべて１より大となり、定常性の条件と反転可能性条件をどちらも満たすことがわかる。

上記で行った残差の診断から、MA(2)は診断検定をパスする。

次数の選択

時系列モデルの古典的な分析では、次のように次数を選択する。

1. 標本PACおよび標本ACの形状からそれぞれAR過程とMA過程の次数を仮に定める。
2. 仮に定めた式を推定する。
3. 推定した式の残差を求め、残差に関して診断検定を行う。
4. 次数を変えて、2と3を繰り返す。

• これは、Box-Jenkins（ボックス=ジェンキンス）法と呼ばれる

1では、次の理論的な性質が背後にある。

• 真の時系列が\(p\)次の自己回帰であれば、\(p\)次を越える真のPAC（母偏自己相関関数）は０になる。
  • PACは、AR(\(p\))の最高次項の係数である。
• 真の時系列が\(q\)次の移動平均であれば、\(q\)次を越える真のAC（母自己相関関数）は０になる。

PACはAR過程の次数を決める際に、ACはMA過程の次数を決める際に、それぞれ重要な情報を与える。

• ある系列に関して、その標本PACが\(p\)次まで有意であり、他方、有意な標本PACがなければ、\(p\)次のAR過程であると判断する。
• ある系列に関して、その標本ACが\(q\)次まで有意であり、他方、有意な標本ACがなければ、\(q\)次のMA過程であると判断する。

経済変数に関するデータ分析でも、標本ACおよびPAC関数を次数を選ぶ際の基準に用いるが、上記のように明確には次数を選べない。 そこで、以下の手順に従って次数を選択することが推奨される。

1. 標本ACとPACの情報を用いて比較的次数の高いAR過程を選び、その式を推定する。
2. 推定されたAR過程の高次項の\(t\)値を調べ、有意でない項を除去し、簡潔な時系列過程の推定を繰り返す。
3. 望ましいと考えられる式について、残差に関する診断を行う。
4. 選ばれたAR過程にMA項を追加して、ARMA過程を推定し、望ましい式についての診断検定を行う。

時系列分析が統計ソフトを使って容易にできるようになった現状では、ARMA式の推定結果にもとづいて時系列過程の選択を進めていくことが実際的である。

• Rでは、forecastパッケージのauto.arima関数を用いて、最小AIC規準の次数選択を行うことができる。

# パッケージの読み込み
library(forecast)

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

# ARMA過程の次数選択
fit <- auto.arima(r, d=0, max.p=4, max.q=4, seasonal = F, ic = "aic", stepwise = T, trace = T)

## 
##  Fitting models using approximations to speed things up...
## 
##  ARIMA(2,0,2)            with non-zero mean : -2477.915
##  ARIMA(0,0,0)            with non-zero mean : -2476.234
##  ARIMA(1,0,0)            with non-zero mean : -2473.599
##  ARIMA(0,0,1)            with non-zero mean : -2474.251
##  ARIMA(0,0,0)            with zero mean     : -2476.538
##  ARIMA(1,0,2)            with non-zero mean : -2478.84
##  ARIMA(0,0,2)            with non-zero mean : -2480.994
##  ARIMA(0,0,3)            with non-zero mean : -2479.805
##  ARIMA(1,0,1)            with non-zero mean : -2474.434
##  ARIMA(1,0,3)            with non-zero mean : -2477.428
##  ARIMA(0,0,2)            with zero mean     : -2481.588
##  ARIMA(0,0,1)            with zero mean     : -2474.542
##  ARIMA(1,0,2)            with zero mean     : -2479.533
##  ARIMA(0,0,3)            with zero mean     : -2480.491
##  ARIMA(1,0,1)            with zero mean     : -2474.623
##  ARIMA(1,0,3)            with zero mean     : -2477.98
## 
##  Now re-fitting the best model(s) without approximations...
## 
##  ARIMA(0,0,2)            with zero mean     : -2481.59
## 
##  Best model: ARIMA(0,0,2)            with zero mean

• オプションで次のように指定している。
  • d = 0：差分を取らない
  • max.p = 4とmax.q = 4：ARとMAの最高次数を4
  • seasonal = F：季節調整は行わない（FはFALSEの略）
  • ic = "aic"：最小AIC規準
  • stepwise = F：選択するモデルが少なくなるように段階的に選択する
    • stepwise = Tとすると、指定した最高次数以下のすべてのモデルを比較するので、特に季節調整のあるモデルを含む場合に時間がかかる。
  • trace = T：比較されたARMAモデルのリストが表示される。

選択されたモデルを表示する。

# 選択されたモデルの表示
summary(fit)

## Series: r 
## ARIMA(0,0,2) with zero mean 
## 
## 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        ACF1
## Training set 0.0009701102 0.01803164 0.01166214 NaN  Inf 0.6433995 0.003341196

• 定数項なしの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/