【統計】時系列分析（AR、MA、ARMA、ARIMA、SARIMA、ARIMAXモデル）

本記事は、AR、MA、ARMA、ARIMA、SARIMA、ARIMAXモデル等の特徴について。

【目次】

計算式等
ノート
参考

計算式等

＜前提＞

■ 弱定常性（Weakly stationary）
　時系列の期待値、分散および自己共分散が常に一定である状態（時点 t および k に依存しない状態）。単に定常性と言ったらこちらを指す。

$\quad E[y_{t}]=\mu_{t}=\mu_{t-k}=\mu$
$\quad Var[y_{t}]=\sigma^{2} _{t}=\sigma^{2} _{t-k}=\sigma^{2}$
$\quad Cov \left(y_{t},\ y_{t-s} \right) = Cov \left(y_{t-k},\ y_{t-s-k} \right)$

■ 強定常性（Strongly stationary）
　任意の時点の集合 {t1, t2, ..., tm} と k に対する時系列の同時密度関数 f が一定の状態（同時分布が同一となる場合）。

$\quad f \left(y_{t_{1}} \ _{,}\ y_{t_{2}} \ _{,}\ ...\ _{,}\ y_{t_{m}} \right) = f \left( y_{t_{1-k}}\ _{,}\ y_{t_{2-k}}\ _{,}\ ...\ _{,}\ y_{t_{m-k}} \right)$

■ ホワイトノイズ（White noise、白色雑音）
　時系列の誤差項。正規分布であり、他の時点と自己相関がない。

$\quad E[\varepsilon _{t}]=0, \quad Var[\varepsilon _{t}]=\sigma^{2}, \quad Cov[\varepsilon _{t},\ \varepsilon _{t-s}]=0$

■ ランダムウォーク（Random walk）
　ホワイトノイズの累積和により得られる時系列。

$\quad y_{t} = \sum \varepsilon _{t}$

＜各モデル＞

■ AR(p)モデル（AutoRegressive model、自己回帰モデル）
　ラグpまでの過去時点（AR(p) 成分）を説明変数に投入して現時点を表すモデル。pは自己回帰の次数、Φは自己回帰係数。

$\quad y_{t} = c + \sum ^{p}_{s=1} \phi _{s} y_{t-s} + \varepsilon _{t}, \quad \varepsilon _{t}\sim W.N.\left( \sigma ^{2}\right)\quad\verb|※W.N.=ホワイトノイズ|$

　※AR(1) の時：
$\qquad y_{t} = c + \phi _{1} y_{t-1} + \varepsilon _{t}$

■ MA(q)モデル（Moving Average model、移動平均モデル）
　ラグqまでの過去のホワイトノイズ（MA(q) 成分）を説明変数に投入して現時点を表すモデル。qは移動平均の次数、θは移動平均係数。

$\quad y_{t} = c + \varepsilon _{t} + \sum ^{q}_{s=1}\theta _{s}\varepsilon_{t-s},\quad \varepsilon _{t}\sim W.N.\left( \sigma ^{2}\right)$

　※MA(1) の時：
$\qquad y_{t} = c + \varepsilon _{t} + \theta _{1}\varepsilon_{t-1}$

■ ARMA(p, q)モデル（AutoRegressive Moving Average model、自己回帰移動平均モデル）
　AR(p) 成分とMA(q) 成分の両方を説明変数に投入して現時点を表すモデル。

$\quad y_{t} = c + \sum ^{p}_{s=1}\phi _{s}y_{t-s} + \varepsilon _{t} + \sum ^{q}_{s=1}\theta _{s}\varepsilon_{t-s},\quad \varepsilon _{t}\sim W.N.\left( \sigma ^{2}\right)$

　※ARMA(1,1) の時：
$\qquad y_{t} = c + \phi _{1}y_{t-1} + \varepsilon _{t} + \theta _{1}\varepsilon_{t-1}$

■ ARIMA(p, d, q)モデル（AutoRegressive Integrated Moving Average model、自己回帰和分移動平均モデル）
　観測値の時系列に対して、近似的に定常な時系列が得られるまで d階の階差（和分プロセス）をとった後に、階差の時系列に対して、ARMA(p,q) モデルを当てはめるモデル。ARIMAモデルは前期と今期のラグ１の d階差をとる。

　和分プロセス：
　[ １階差 ]
$\quad \Delta_{1} y_{t} = y_{t} - y_{t-1}$

　[ ２階差 ]
$\quad \Delta _{1} \Delta _{1} y_{t} = \Delta _{1} y_{t} - \Delta _{1} y_{t-1} = \left(y_{t} - y_{t-1}\right) - \left(y_{t-1} - y_{t-1-1}\right) = y_{t} - 2 y_{t-1} + y_{t-2}$

　※ARIMA(p, 1, q) の時：
$\qquad z_{t} = \Delta _{1} y_{t}$
$\qquad z_{t} = c + \sum ^{p}_{s=1}\phi _{s} z_{t-s} + \varepsilon _{t} + \sum ^{q}_{s=1}\theta _{s}\varepsilon_{t-s},\quad \varepsilon _{t}\sim W.N.\left( \sigma ^{2}\right)$

　※ARIMA(1, 1, 1) の時：
$\qquad z_{t} = c + \phi _{1} z_{t-1} + \varepsilon _{t} + \theta _{1}\varepsilon_{t-1}$

※追記：SARIMAとARIMAXを加えます。
■ SARIMA(p, d, q)(P, D, Q)[m] モデル（Seasonal AutoRegressive Integrated Moving Average model、季節自己回帰和分移動平均モデル）
　ARIMAモデルと同じ要領でラグ m の季節階差（D階差）をとった後にARMAを当てはめるモデル。ラグ演算子（L^m）を用いて表現する。ラグ演算子の上付き文字はラグを表す（L^m y_t = y_t-m）。

$\quad \phi \left( L \right) \Phi \left( L^{m}\right) \Delta _{m}^{D} \Delta _{1}^{d} y_{t} = c + \theta \left( L\right) \Theta \left( L^{m}\right) \varepsilon _{t}$

$\qquad \phi \left( L \right) = 1 - \phi_{1} L^{1} - \phi_{2} L^{2} - ... - \phi_{p} L^{p} = 1 - \sum ^{p}_{k=1} \phi_{k} L^{k}$
$\qquad \Phi \left( L^{m} \right) = 1 - \Phi_{1} L^{m} - \Phi_{2} L^{2m} - ... - \Phi_{P} L^{P\cdot m} = 1 - \sum ^{P}_{k=1} \Phi_{k} L^{k\cdot m}$
$\qquad \theta \left( L \right) = 1 + \theta_{1} L^{1} + \theta_{2} L^{2} + ... + \theta_{q} L^{q} = 1 + \sum ^{q}_{k=1} \theta_{k} L^{k}$
$\qquad \Theta \left( L^{m} \right) = 1 + \Theta_{1} L^{m} + \Theta_{2} L^{2m} + ... + \Theta_{Q} L^{Q\cdot m} = 1 + \sum ^{Q}_{k=1} \Theta_{k} L^{k\cdot m}$

　 $\theta \left( L \right)$ と $\Theta \left( L^{m} \right)$ について、演算子をプラスで示している文献とマイナスで示している文献の両方ある。ホワイトノイズは正規分布を仮定するためどちらでも良いということだろうか。

　※SARIMA(1, 1, 1) (1, 1, 1) [12] の時（階差と12周期季節階差の２階差）：
$\qquad z_{t}=\Delta _{12}^{1} \Delta _{1}^{1} y_{t} = \Delta _{12}^{1} y_{t} - \Delta _{12}^{1} y_{t-1} = \left(y_{t} - y_{t-12}\right) - \left(y_{t-1} - y_{t-1-12}\right) = y_{t} - y_{t-12} - y_{t-1} + y_{t-13}$

$\qquad \phi \left( L \right) \Phi \left( L^{12}\right) z_{t} = c + \theta \left( L\right) \Theta \left( L^{12}\right) \varepsilon _{t}$
$\qquad \Leftrightarrow \left( 1 - \phi_{1} L^{1} \right) \left( 1 - \Phi_{1} L^{12} \right) z_{t} = c + \left( 1 + \theta_{1} L^{1} \right) \left( 1 + \Theta_{1} L^{12} \right) \varepsilon _{t}$
$\qquad \Leftrightarrow \left( 1 - \phi_{1} L^{1} - \Phi_{1} L^{12} + \phi_{1} \Phi_{1} L^{13} \right) z_{t} = c + \left( 1 + \theta_{1} L^{1} + \Theta_{1} L^{12} + \theta_{1} \Theta_{1} L^{13} \right) \varepsilon _{t}$
$\qquad \Leftrightarrow z_{t} - \phi_{1} z_{t-1} - \Phi_{1} z_{t-12} + \phi_{1} \Phi_{1} z_{t-13} = c + \varepsilon _{t} + \theta_{1} \varepsilon _{t-1} + \Theta_{1} \varepsilon _{t-12} + \theta_{1} \Theta_{1} \varepsilon _{t-13}$
$\qquad \Leftrightarrow z_{t} = c + \phi_{1} z_{t-1} + \Phi_{1} z_{t-12} - \phi_{1} \Phi_{1} z_{t-13} + \varepsilon _{t} + \theta_{1} \varepsilon _{t-1} + \Theta_{1} \varepsilon _{t-12} + \theta_{1} \Theta_{1} \varepsilon _{t-13}$

　※季節の自己回帰係数 $\Phi_{1}$ と移動平均係数 $\Theta_{1}$ を０にすれば、ひとつ上のARIMAの例で取り上げたARIMA(1, 1, 1) と同じになる。
$\qquad z_{t} = c + \phi _{1} z_{t-1} + \varepsilon _{t} + \theta _{1}\varepsilon_{t-1}$

■ ARIMAXモデル（AutoRegressive Integrated Moving Average with eXogenous variables、外生変数付き自己回帰和分移動平均モデル）
　ARIMAモデルに回帰成分 β x （外生変数、eXogenous vriables）を組み込んだモデル。例として閏日や祝日などのカレンダー効果、気温・降水量などの天候、災害や制度変更での水準変動を表す干渉シフト（一定期間フラグを立てる）などを導入する。

$\quad z_{t} = c + \beta x_{t} + \sum ^{p}_{s=1}\phi _{s}z_{t-s} + \varepsilon _{t} + \sum ^{q}_{s=1}\theta _{s}\varepsilon_{t-s},\quad \varepsilon _{t}\sim W.N.\left( \sigma ^{2}\right)$

ノート

上に示した過程はARIMAモデルで表すことができる（SARIMA、ARIMAX除く）。

・ARIMA(0,0,0)：ホワイトノイズ
・ARIMA(0,1,0)：ランダムウォーク
・ARIMA(p,0,0)：AR(p)モデル
・ARIMA(0,0,q)：MA(q)モデル
・ARIMA(p,0,q)：ARMA(p,q)モデル
・ARIMA(p,d,q)：ARIMA(p,d,q)モデル

また、MA(q) モデルは時点 t に依存しないため常に定常だが、
AR(p) モデルは、次に示す特性方程式の解（複素解含む）の絶対値がすべて１未満であることが、定常性の必要十分条件と言われる。

$\quad 1 = \sum \phi_{s} z^{s} = \phi_{1} z^{1} + \phi_{2} z^{2} + ... + \phi_{s} z^{s}$

AR(1) モデルなら 1 ＜ |Φ1| の場合。つまり、Φ1 が１を超えると非定常になってしまう。

ここからは、自己回帰係数Φと平均移動係数θの違いによる振る舞いを見てみる。

次のように適当な定数項 c とホワイトノイズ（ =NORM.INV(RAND(),0,0.5) ）を準備し、n=150 の観測値の時系列を作成。
※Excelでの例示は雑なのでイメージ程度にとどめてほしい。また、ARモデルの初期値は０にしている。
f:id:cochineal19:20210307162430p:plain

この観測値の時系列に対して、AR(1) とMA(1) モデルを当てはめて、自己回帰係数Φと平均移動係数θによる変化を見てみる。

f:id:cochineal19:20210307171843p:plain
f:id:cochineal19:20210307171921p:plain
f:id:cochineal19:20210307171951p:plain
f:id:cochineal19:20210307172016p:plain
f:id:cochineal19:20210307172120p:plain
f:id:cochineal19:20210307172145p:plain
f:id:cochineal19:20210307172215p:plain
f:id:cochineal19:20210307172238p:plain

先ほどの AR(1) モデルの定常性条件（1 ＜ |Φ1| ）を頭に入れて見てみると、自己回帰係数Φの絶対値が１以上の Φ=1.1 と Φ=－1.1 で非定常になっていることが良く分かる（予測モデルとして使えない）。
また、MA(1) モデルは定数項 c 付近でバラつく定常な時系列であることが分かる。

次に、先ほどの観測値の時系列にトレンド（時点毎に 0.5 増える）を加えてみる。このトレンド付きの時系列について、ラグ１の１階差（d=1）をとったときの振る舞いを見てみる。
※こちらも例示は雑なのでイメージ程度にとどめてほしい。また、初期値は０にしている。
f:id:cochineal19:20210307165011p:plain

１階差をとった後、階差の時系列が定常になったことが分かる。

今回、AR(1) とMA(1) モデルでいくつかと１階差を１つ例示したが、ARIMA(p,d,q) でどのようなパラメタ設定が良いか（ベストモデルは何か）は、AICやRMSE等の評価指標を用いて判断する。
r では forecastパッケージでここらができるので次回記事にしたい。
SAS にも proc ARIMA があるので余裕があれば。

なお、ARIMA(p,d,q) に従う観測値の時系列は、r の stats::arima.sim関数で比較的簡単に作れるので備忘録がてらメモする。

stats::arima.sim(model=list(order=c(p,d,q), ar=c(x,...), ma=c(x,...), n=x))

AR(p)、MA(q) に従う観測値の時系列：

library(tidyverse)
set.seed(1)
ads1 <- data.frame()
ads1 <- union_all(ads1, data.frame(x=1:200, AVAL=arima.sim(model=list(order=c(1,0,0), ar=c(0.3)),         n=200), TYPE="01. AR(1)  Φ1 = 0.3"))
ads1 <- union_all(ads1, data.frame(x=1:200, AVAL=arima.sim(model=list(order=c(2,0,0), ar=c(0.3,0.5)),     n=200), TYPE="02. AR(2)  Φ1 = 0.3, Φ2 = 0.5"))
ads1 <- union_all(ads1, data.frame(x=1:200, AVAL=arima.sim(model=list(order=c(3,0,0), ar=c(0.3,0.5,0.1)), n=200), TYPE="03. AR(3)  Φ1 = 0.3, Φ2 = 0.5, Φ3 = 0.1"))
ads1 <- union_all(ads1, data.frame(x=1:200, AVAL=arima.sim(model=list(order=c(0,0,1), ma=c(0.3)),         n=200), TYPE="04. MA(1)  θ1 = 0.3"))
ads1 <- union_all(ads1, data.frame(x=1:200, AVAL=arima.sim(model=list(order=c(0,0,2), ma=c(0.3,0.5)),     n=200), TYPE="05. MA(2)  θ1 = 0.3, θ2 = 0.5"))
ads1 <- union_all(ads1, data.frame(x=1:200, AVAL=arima.sim(model=list(order=c(0,0,3), ma=c(0.3,0.5,0.1)), n=200), TYPE="06. MA(3)  θ1 = 0.3, θ2 = 0.5, θ3 = 0.1"))

ggplot(ads1, aes(x=x, y=AVAL)) +
  geom_line() + facet_wrap(. ~ TYPE, ncol=3)

f:id:cochineal19:20210307180214p:plain

ARMA(p,q) に従う観測値の時系列：

set.seed(1)
ads2 <- data.frame()
ads2 <- union_all(ads2, data.frame(x=1:200, AVAL=arima.sim(model=list(order=c(1,0,1), ar=c(0.3),         ma=c(0.3)),         n=200), TYPE="01. ARIMA(1,0,1) \n Φ1 = 0.3 \n θ1 = 0.3"))
ads2 <- union_all(ads2, data.frame(x=1:200, AVAL=arima.sim(model=list(order=c(2,0,2), ar=c(0.3,0.5),     ma=c(0.3,0.5)),     n=200), TYPE="02. ARIMA(2,0,2) \n Φ1 = 0.3, Φ2 = 0.5 \n θ1 = 0.3, θ2 = 0.5"))
ads2 <- union_all(ads2, data.frame(x=1:200, AVAL=arima.sim(model=list(order=c(3,0,3), ar=c(0.3,0.5,0.1), ma=c(0.3,0.5,0.1)), n=200), TYPE="03. ARIMA(3,0,3) \n Φ1 = 0.3, Φ2 = 0.5, Φ3 = 0.1 \n θ1 = 0.3, θ2 = 0.5, θ3 = 0.1"))
ggplot(ads2, aes(x=x, y=AVAL)) +
  geom_line() + facet_wrap(. ~ TYPE, ncol=3)

f:id:cochineal19:20210307181223p:plain
ARIMA(p,d,q) に従う観測値の時系列：

set.seed(1)
ARIMA212 <- data.frame(y=arima.sim(model=list(order=c(2,1,2), ar=c(0.3,0.5), ma=c(0.3,0.5)), n=200), id=1:201)
ggplot(data=ARIMA212, mapping=aes(x=id, y=y)) +
  geom_line() +
  labs(title="ARIMA(2,1,2)", subtitle="（AR: Φ1 = 0.3, Φ2 = 0.5,  MA: θ1 = 0.3, θ2 = 0.5）"
       ,x="時点", y="", caption="🄫Cochineal19") +
  theme(plot.title=element_text(size=30)
        ,plot.subtitle=element_text(size=25)
        ,axis.title=element_text(size=25)
        ,axis.text=element_text(size=20))
 
set.seed(1)
ARIMA222 <- data.frame(y=arima.sim(model=list(order=c(2,2,2), ar=c(0.3,0.5), ma=c(0.3,0.5)), n=200), id=1:202)
ggplot(data=ARIMA222, mapping=aes(x=id, y=y)) +
  geom_line() +
  labs(title="ARIMA(2,2,2)", subtitle="（AR: Φ1 = 0.3, Φ2 = 0.5,  MA: θ1 = 0.3, θ2 = 0.5）"
       ,x="時点", y="", caption="🄫Cochineal19") +
  theme(plot.title=element_text(size=30)
        ,plot.subtitle=element_text(size=25)
        ,axis.title=element_text(size=25)
        ,axis.text=element_text(size=20))