統計学といくつかのよしなしごと

つれづれなるままに、日くらし統計学にむかひて、心にうつりゆくよしなし事を、そこはかとなく書きつくれば、あやしうこそものぐるほしけれ。

時系列分析の理論と実装を概観する①ARIMA過程まで

時系列 統計学 回帰

1.時系列分析における個人的諸混乱

　実務ではしょっちゅう時系列データに出くわすが、これを使って動向の予測やメカニズムの解釈をしたいと思ったときのデータの扱いが難しくて混乱しがちだ。ただ正しく使えば実務上大いに役に立つし、改めて回帰の枠組みを見直し理解するうえでも良い教材になる。ということでこの記事では時系列分析の標準的な手法であるARIMAモデルの利用*1を前提に、分析の基礎的な理論部分と実装手順を概観する。主な参考文献は以下。

経済・ファイナンスデータの計量時系列分析 (統計ライブラリー)

経済・ファイナンスデータの計量時系列分析 (統計ライブラリー)

時系列分析と状態空間モデルの基礎: RとStanで学ぶ理論と実装

時系列分析と状態空間モデルの基礎: RとStanで学ぶ理論と実装

現場ですぐ使える時系列データ分析 ~データサイエンティストのための基礎知識~

現場ですぐ使える時系列データ分析 ~データサイエンティストのための基礎知識~

2.1変量時系列データの分析概要

　時系列データに眺めた際には、「今日の株価は下落したから明日は買いが増えて株価は上昇するだろう*2」とか「直近5年間で売上は順調に伸びているから来年も伸びるだろう」とか言った具合に、過去の実現値が未来の値に影響する(時間依存する)ことを期待する。これを『自己相関』と言って、通常の回帰と枠組みが異なってくる*3。分析においてはこの自己相関をうまくモデル化してやる必要があるが、基本的なアプローチに過去の自分自身に回帰させる方法、 {\displaystyle y_t }と {\displaystyle y_{t-1} }の間に共通の成分を含めてやる方法の2通りがある。前者を『自己回帰(AR)過程(autoregressive process)』、後者を『移動平均(MA)過程(moving average process)』という。それぞれの式の形を以下に示しつつ、基本になるいくつかのモデルと押さえるべき諸条件を確認する。

AR(p)過程

{\displaystyle y_t = c + φ_1 y_{t-1} + φ_2 y_{t-2} + φ_3 y_{t-3} + ... + φ_p y_{t-p} + ε_t, \ \ \ ε_t \ ～ \ W.N.(σ^2) }
　
　過去p時点までの自分自身の線形和で表現するため、p次のAR過程と言ったりする。ここで {\displaystyle ε_t }は時刻t時点における撹乱項(誤差項)であり、ホワイトノイズ( {\displaystyle } W.N. )と呼ばれる以下を満たす確率分布である。
{\displaystyle E(ε_t) = 0 \\ \\ γ_k = E(ε_t ε_{t-k}) = \left\{ \begin{array}{ll} σ^2, & k = 0 \\ 0, & k ≠ 0 \end{array} \right. }
　一見通常の回帰で誤差項に対して置く仮定と同様のようだが、独立で同一の分布に従うことを要求していない点がポイントになっている。これは実務上この仮定を必ずしも必要としないことや、実際のデータがこれを満たすと考えにくいことが背景にあるよう。
　AR(p)過程の性質についてここでは詳しく述べないが、ポイントとしては自己相関が指数的に減衰する、つまりラグが大きければ大きい程自己相関が弱くなることと、偏自己相関*4がp+1以降(理論上は)0になることである*5。

MA(q)過程

{\displaystyle y_t = μ + ε_t + θ_1 ε_{t-1} + θ_2 ε_{t-2} + θ_3 ε_{t-3} + ... + θ_q ε_{t-q}, \ \ \ ε_t \ ～ \ W.N.(σ^2) }

　こちらのアプローチでは {\displaystyle y_t } を過去の撹乱項の線形和で表現する。時点に限らず撹乱項がホワイトノイズであることに注意すると、q+1以降の自己相関が0になることが分かる。また、後述するがMA(q)過程は特定条件を満たすときAR(∞)過程で書き直すことができ、これにより偏自己相関はラグの大きさに応じて減衰していく。

ARMA(p, q)過程

　AR(p)過程とMA(q)過程を組み合わせたものを『自己回帰移動平均(ARMA)過程(autoregressive moving average process)』という。
{\displaystyle y_t = c + φ_1 y_{t-1} + ... + φ_p y_{t-p} + θ_1 ε_{t-1} + ... + θ_q ε_{t-q}, \ \ \ ε_t \ ～ \ W.N.(σ^2) }

　ARMA(p, q)はその名の通り両方の過程の性質を併せ持つ過程で、いずれか強いほうの性質に従う。そのため自己相関、偏自己相関はともにラグの大きさに応じて減衰していく。

定常性の仮定

　ところで、これらのモデルは『定常性』の仮定を満たす必要がある。定常性には『弱定常性』と『強定常性』があるが、ここでは弱定常性のみ記載する。強定常性は弱定常性により強い条件を課した概念である(断りなく定常性という場合は普通弱定常性を指す)。

{\displaystyle 弱定常性：任意のtとkに対して、 \\ E(y_t) = μ \\ Cov(y_t, \ y_{t-k}) = E [(y_t - μ)(y_{t-k} - μ) ] = γ_k }

　要は弱定常仮定の下では時点を問わず {\displaystyle y_t } の期待値が一定で、かつ自己共分散(よって自己相関も)時点に依存せず、時間差 {\displaystyle k } のみに依存することを課す条件である。この定常性を満たした時系列モデルを基礎として構築したうえで、季節性とかトレンドとか各種の性質をもった実際の時系列データにアプローチしていくことになる。
　定常性を条件に課す背景には、時系列データの特殊性がある。通常の場合ではある確率変数について複数の実現値があるため、それらを用いて平均なり分散なりといった母数を推定するが、時系列データの場合、1つの時点に対して1つの実現値しか存在しない。どういうことかというと、「2018年8月31日の店舗Aの売上は100万円から300万円の可能性があったが(確率変数)、実現値としてはただ1つ120万円が得られた」というようなイメージ。1つの確率変数に対して1つの値しかないため、分散からもう推定できなくなる。そこで、この定常性を課すことで得られた実現値全体からなる一つの過程を考えパラメータの推定を可能にしようというわけである。

定常性と反転可能性の条件

　時系列モデルではまずこの定常性を満たすかどうかの確認から始まるが、①MA過程は常に定常になるため、②AR過程(またはARMA過程のAR過程部分)について定常性を満たすかどうかを検討すればよい。
　①について。MA(q)過程のモデル式は以下であった。
{\displaystyle y_t = μ + ε_t + θ_1 ε_{t-1} + θ_2 ε_{t-2} + θ_3 ε_{t-3} + ... + θ_q ε_{t-q}, \ \ \ ε_t \ ～ \ W.N.(σ^2) }

　これについて、
{\displaystyle \begin{eqnarray} E(y_t) &=& E(μ + ε_t + θ_1 ε_{t-1} + θ_2 ε_{t-2} + θ_3 ε_{t-3} + ... + θ_q ε_{t-q}) \\ &=& E(μ) + E(ε_t) + θ_1 E(ε_{t-1}) + ... + θ_q E(ε_{t-q}) \\ &=& E(μ) \\ &=& μ \end{eqnarray} }
{\displaystyle \begin{eqnarray} γ_k = Cov(y_t, \ y_{t-k}) &=& E [ (y_t - μ)(y_{t-k} - μ) ] \\ &=& E [ (ε_t + θ_1 ε_{t-1} + ... + θ_q ε_{t-q})(ε_{t-k} + θ_1 ε_{t-1-k} + ... + θ_q ε_{t-q-k}) ] \\ &=& E [ (ε_t + θ_1 ε_{t-1} + ... + θ_k ε_{t-k} + ... + θ_q ε_{t-q})(ε_{t-k} + θ_1 ε_{t-1-k} + ... + θ_{q-k} ε_{t-q} + ... + θ_q ε_{t-q-k}) ] \\ &=& (θ_k + θ_1 θ_{k+1} + ... + θ_{q-k} θ_q)σ^2 \\ \end{eqnarray} }
{\displaystyle (ただし1 \leq k \leq q, \ \ \ k \geq q+1 ではγ_k = 0) }

　となって、定常性の条件を満たしていることが分かる。

　②について。AR(p)過程のモデル式は以下であった。
{\displaystyle y_t = c + φ_1 y_{t-1} + φ_2 y_{t-2} + φ_3 y_{t-3} + ... + φ_p y_{t-p} + ε_t, \ \ \ ε_t \ ～ \ W.N.(σ^2) }

　これについて、以下で表される『AR特性方程式』のすべての解の絶対値が1より大きくなる( {\displaystyle |z| > 1 } )ときAR過程は定常となる。ここの数理はしっかりとは追っていないが、 {\displaystyle y_t } の分散を有限にするために必要な制約っぽい。

{\displaystyle 1 - φ_1 z - ... - φ_p z^p = 0 }

　この下でAR(p)過程の期待値と自己共分散を確認しておく。定常であれば {\displaystyle μ = E(y_t) = E(y_{t-1} = ... = E(y_{t-p})) } であることに注意して、期待値は
{\displaystyle \begin{eqnarray} μ = E(y_t) &=& E(c + φ_1 y_{t-1} + φ_2 y_{t-2} + φ_3 y_{t-3} + ... + φ_p y_{t-p} + ε_t) \\ &=& E(c) + φ_1 E(y_{t-1}) + φ_2 E(y_{t-2}) + ... + φ_p E(y_{t-p}) + E(ε_t) \\ &=& c + (φ_1 + φ_2 + ... + φ_p)μ \\ \end{eqnarray} }
{\displaystyle \therefore μ = E(y_t) = \frac {c}{1 - φ_1 - φ_2 - ... - φ_p} }

一方自己共分散は、
{\displaystyle \begin{eqnarray} γ_k &=& Cov(y_t, \ y_{t-k}) \\ &=& Cov(c + φ_1 y_{t-1} + φ_2 y_{t-2} + φ_3 y_{t-3} + ... + φ_p y_{t-p} + ε_t, \ y_{t-k}) \\ &=& Cov(c, \ y_{t-k}) + Cov(φ_1 y_{t-1}, \ y_{t-k}) + Cov(φ_2 y_{t-2}, \ y_{t-k}) + ... + Cov(φ_p y_{t-p}, \ y_{t-k}) + Cov(ε_t, \ y_{t-k}) \\ &=& φ_1 γ_1 + φ_2 γ_2 + ... + φ_p γ_p \end{eqnarray} }

　ちなみに、この両辺を分散 {\displaystyle γ_0 = Var(y_t) } で割って自己相関の方程式に変換したもの*6を『ユール・ウォーカー方程式』という。グラフなど描くと自己相関の絶対値が減衰していくことが確認できる。

　ここまで定常性について見てきたが、パラメータの推定という観点でもう一つ気にしなければならないことにMA過程の『反転可能性』というものがある。これはMA過程をAR(∞)に書き直せるかどうかという話で、MA過程においては同じ期待値と自己相関構造を持つパラメータの組が複数ありうるので、その同定のために課す条件である。反転可能、つまり撹乱項 {\displaystyle ε_t } を {\displaystyle y_t } と関数として表現できるならば、過去の {\displaystyle y } で {\displaystyle y_t } を予測したときの予測誤差と解釈できて自然なので、反転可能性を満たすパラメータの組を選択しようという判断基準が与えられる。反転可能条件は以下のMA特性方程式のすべての解の絶対値が1より大きい( {\displaystyle |z| > 1 } )ことである。

{\displaystyle 1 + θ_1 z + ... + θ_q z^q = 0 }

　このあたりを満たすようなパラメータを導いてやろうというのが、時系列データに対するARMAモデルの推定の考え方となる。

非定常過程と単位根とARIMA過程

　ARMAはデータ系列の定常性の下で推定や予測を行うが、現実のデータは非定常な過程が多く、その場合定常過程とは性質が大きく異なってくる。有名なものに『単位根過程』があり、以下で定義される(「経済・ファイナンスデータの計量時系列分析」P105)。

原系列*7 {\displaystyle y_t } が非定常過程であり、差分系列 {\displaystyle Δy_t = y_t - y_{t-1} } が定常過程であるとき、過程は単位根過程(unit root process)といわれる。

　今回は単位根過程の性質については詳しく述べないが、予測のMSE(平均二乗誤差)が長期的に増大していったり、二変数間のモデル化において『見せかけの回帰*8』を引き起こしたりといろいろ厄介な挙動を示すので追加の処理を検討する必要が出てくる。単位根過程は『1次和分過程』とも言われ、和分過程が何かというと(「経済・ファイナンスデータの計量時系列分析」P106)、

d-1階差分をとった系列は非定常過程であるが、d階差分をとった系列が定常過程に従う過程は、d次和分過程もしくはI(d)過程と呼ばれる。また、I(0)は定常過程で定義される。

　差分系列が定常になるならそれはARMA過程として表現できるので、何階差分をとれば定常になるかを付与したものが『自己回帰和分移動平均(ARIMA)過程』と呼ばれる過程である。定義は以下の通り(「経済・ファイナンスデータの計量時系列分析」P106)。

d階差分をとった系列が定常かつ反転可能なARMA(p, q)過程に従う過程は次数(p, d, q)の自己回帰和分移動平均過程もしくはARIMA(p, d, q)過程と呼ばれる。

　次に、こうした基礎を頭に入れつつ分析の手順とそれに伴う結果を確認してみる。

3.Rによる分析の手順と結果の確認

　Rに用意されているNileデータを使って分析の手順を確認する。Nileは1871年から1970年までの100年間におけるナイル川の流量推移のデータ。原系列のプロットはこんな感じ。

> plot(Nile)

f:id:yanbow221:20180830161535j:plain

　実際の予測時に見立てるために最後の20年分は予測期間として、最初80年分を手元にあるデータとして以降分析していく。

> library(tidyverse)
> library(tseries)
> library(forecast)

> Nile_train <- window(Nile, end = 1950)     #学習データ
> Nile_test <- window(Nile, start = 1951)    #テストデータ
モデルの推定

　はじめに単位根検定を行い原系列が単位根過程かそうでないかを確認する。単位根を持っている場合、例えば差分系列をとってやって定常過程にするといった追加処理を検討する。単位根検定には結構種類があってその詳細は調べていないが、オーソドックスな方法としてここではADF検定を用いる。ADF検定は真のモデルをAR(p)過程として、帰無仮説に「単位根を持つ」、対立仮説に「単位根を持たない」を置いた検定である。p値は0.087で有意水準5％下で帰無仮説を棄却できない、すなわち単位根を持ってないとはいえないなぁという結果になる。

> adf.test(Nile_train)　#ADF検定

	Augmented Dickey-Fuller Test

data:  Nile_train
Dickey-Fuller = -3.2443, Lag order = 4, p-value = 0.08675
alternative hypothesis: stationary

　ということで差分系列を作成し、これに対して再び単位根検定を行ってみると、

> Nile_train_diff <- diff(Nile_train, lag = 1)  #差分系列をとる
> adf.test(Nile_train_diff)                     #差分系列に対してADF検定

	Augmented Dickey-Fuller Test

data:  Nile_train_diff
Dickey-Fuller = -5.8527, Lag order = 4, p-value = 0.01
alternative hypothesis: stationary

Warning message:
In adf.test(Nile_train_diff) : p-value smaller than printed p-value

　今度は有意な結果を示したので、この1階差分の系列に対して分析を進める。真の過程が不明な状態なので、各過程を真としたときの自己相関、偏自己相関の性質から仮定するモデルを選択する。「経済・ファイナンスデータの計量時系列分析」P49にある表を抜粋しておく。
f:id:yanbow221:20180830165637p:plain

　これを判断するために差分系列の自己相関、偏自己相関をプロットするとこんな感じ。

> par(mfrow = c(2, 1))
> acf(Nile_train_diff)                     #自己相関プロット
> acf(Nile_train_diff, type = "partial")   #偏自己相関プロット

f:id:yanbow221:20180830175425j:plain

　自己相関(ACF)については1次の自己相関が有意であり、2次以降は減衰というよりは切断されている(0とみなす)気もするが、とりあえずどちらの可能性も持たせておく。偏自己相関(PACF)については1次が有意で2次が微妙だが帰無仮説を棄却できておらず、3次以降は(10次時点が有意になっているが)減衰ないし切断ととらえて、1次までをモデルに盛り込むことを検討する。というわけで過程の候補をAR(1)、MA(1)、ARMA(1, 1)とする。これは差分系列に対してのものなので、それぞれARIMA(1, 1, 0)、ARIMA(0, 1, 1)、ARIMA(1, 1, 1)と同等となる。
　上で選択した候補においてパラメータの推定を行っていく。それぞれ定数項のありなし2種類において推定するので、計6つの過程において推定を行う。どのモデルが優れているかについて今回はAIC(赤池情報量基準)を用いる。

> # パラメータ推定
> model_ar1 <- arima(Nile_train_diff, order = c(1, 0, 0), include.mean = F)      #定数項なしAR(1)
> model_ar1_intercept <- arima(Nile_train_diff, order = c(1, 0, 0))              #定数項ありAR(1)
> 
> model_ma1 <- arima(Nile_train_diff, order = c(0, 0, 1), include.mean = F)      #定数項なしMA(1)
> model_ma1_intercept <- arima(Nile_train_diff, order = c(0, 0, 1))              #定数項ありMA(1)
> 
> model_arma1.1 <- arima(Nile_train_diff, order = c(1, 0, 1), include.mean = F)  #定数項なしARMA(1, 1)
> model_arma1.1_intercept <- arima(Nile_train_diff, order = c(1, 0, 1))          #定数項ありARMA(1, 1)

> # AICの比較
> result_aic <- as.matrix(c(model_ar1$aic, model_ar1_intercept$aic,
+                           model_ma1$aic, model_ma1_intercept$aic,
+                           model_arma1.1$aic, model_arma1.1_intercept$aic), ncol = 1)
> rownames(result_aic) <- c("AR(1)定数項なし", "AR(1)定数項あり",
+                           "MA(1)定数項なし", "MA(1)定数項あり",
+                           "ARMA(1,1)定数項なし", "ARMA(1,1)定数項あり")
> colnames(result_aic) <- "AIC"
> result_aic
                         AIC
AR(1)定数項なし     1027.300
AR(1)定数項あり     1029.236
MA(1)定数項なし     1018.066
MA(1)定数項あり     1019.439
ARMA(1,1)定数項なし 1016.574
ARMA(1,1)定数項あり 1015.936

　AICが最小になるモデルはARMA(1, 1)過程となった。このモデルの下での推定パラメータは以下の通り。

> model_arma1.1_intercept

Call:
arima(x = Nile_train_diff, order = c(1, 0, 1))

Coefficients:
         ar1      ma1  intercept
      0.3396  -0.9997    -4.2092
s.e.  0.1074   0.0664     1.0029

sigma^2 estimated as 19426:  log likelihood = -503.97,  aic = 1015.94

　実はここまでの流れは、Rにおいてはパッケージforecastに入っているauto.arimaが一発でやってくれる。auto.arimaを使うと総当たり的にモデルを推定し、指定した情報量基準の下で最適なモデルを選択してくれる。そのまま原系列に対してやってみるとこんな感じ。

> model_nile <- auto.arima(Nile_train, ic = "aic", trace = T, stepwise = F, approximation = F)

 ARIMA(0,1,0)                    : 1038.336
 ARIMA(0,1,0) with drift         : 1040.313
 ARIMA(0,1,1)                    : 1018.066
 ARIMA(0,1,1) with drift         : 1019.439
 ARIMA(0,1,2)                    : 1017.266
 ARIMA(0,1,2) with drift         : 1017.486
 ARIMA(0,1,3)                    : 1018.316
 ARIMA(0,1,3) with drift         : Inf
 ARIMA(0,1,4)                    : 1020.194
 ARIMA(0,1,4) with drift         : Inf
 ARIMA(0,1,5)                    : 1022.184
 ARIMA(0,1,5) with drift         : Inf
 ARIMA(1,1,0)                    : 1027.3
 ARIMA(1,1,0) with drift         : 1029.236
 ARIMA(1,1,1)                    : 1016.574
 ARIMA(1,1,1) with drift         : Inf
 ARIMA(1,1,2)                    : 1018.371
 ARIMA(1,1,2) with drift         : Inf
 ARIMA(1,1,3)                    : 1020.233
 ARIMA(1,1,3) with drift         : Inf
 ARIMA(1,1,4)                    : Inf
 ARIMA(1,1,4) with drift         : Inf
 ARIMA(2,1,0)                    : 1025.79
 ARIMA(2,1,0) with drift         : 1027.705
 ARIMA(2,1,1)                    : 1018.308
 ARIMA(2,1,1) with drift         : Inf
 ARIMA(2,1,2)                    : 1020.333
 ARIMA(2,1,2) with drift         : Inf
 ARIMA(2,1,3)                    : Inf
 ARIMA(2,1,3) with drift         : Inf
 ARIMA(3,1,0)                    : 1025.755
 ARIMA(3,1,0) with drift         : 1027.631
 ARIMA(3,1,1)                    : 1020.253
 ARIMA(3,1,1) with drift         : Inf
 ARIMA(3,1,2)                    : Inf
 ARIMA(3,1,2) with drift         : Inf
 ARIMA(4,1,0)                    : 1025.598
 ARIMA(4,1,0) with drift         : 1027.434
 ARIMA(4,1,1)                    : 1022.102
 ARIMA(4,1,1) with drift         : Inf
 ARIMA(5,1,0)                    : 1025.951
 ARIMA(5,1,0) with drift         : 1027.72



 Best model: ARIMA(1,1,1)

　こちらではドリフト率(先ほどまで述べていた定数項に当たるもの)なしのARIMA(1, 1, 1)が提案された。先ほどは定数項ありのARMA(1, 1)、すなわちドリフト率ありのARIMA(1, 1, 1)を選択したが、auto.arimaでは結果が「Inf」となっている。これはどうもAR過程の定常性やMA過程の反転可能性が危うい、つまり特性方程式の解が1より大きくなるかどうか怪しいときに出るようで、実際過程における特性方程式の解を確認してみると、反転可能性がちょっと怪しいという結果が確認できる。

> abs(polyroot(c(1, -coef(model_arma1.1_intercept)["ar1"])))  #差分系列における定数項ありARMA(1,1)のAR特性方程式の解
[1] 2.944355
> abs(polyroot(c(1, coef(model_arma1.1_intercept)["ma1"])))   #差分系列における定数項ありARMA(1,1)のMA特性方程式の解
[1] 1.000263
> abs(polyroot(c(1, -coef(model_nile)["ar1"])))               #原系列におけるドリフト率なしARIMA(1,1,1)のAR特性方程式の解
[1] 3.768194
> abs(polyroot(c(1, coef(model_nile)["ma1"])))                #原系列におけるドリフト率なしARIMA(1,1,1)のMA特性方程式の解
[1] 1.150311

　実務ではauto.arimaに頼ることが多いと思うので、こちらで選択されたモデル、原系列におけるARIMA(1, 1, 1)過程について以降のプロセスを進めていく。モデルの定常性、反転可能性は既に確認したので、次に残差について確認する。選択したモデルで原系列の自己相関が十分に説明されているか？という意味で残差 {\displaystyle e_t = y_t - \hat{y_t} } の自己相関が0になっていることを確認する。また、モデルの撹乱項に置いたホワイトノイズ自体は正規分布を仮定しないが、一般に撹乱項は正規性を仮定することが多い(「時系列分析と状態空間モデルの基礎」p30, 68)ので、併せて {\displaystyle e_t } の正規性も検定する。

> checkresiduals(model_nile)             #残差の自己相関の検定

	Ljung-Box test

data:  Residuals from ARIMA(1,1,1)
Q* = 9.6295, df = 8, p-value = 0.292

Model df: 2.   Total lags used: 10

> jarque.bera.test(resid(model_nile))    #残差の正規性の検定

	Jarque Bera Test

data:  resid(model_nile)
X-squared = 0.42784, df = 2, p-value = 0.8074

f:id:yanbow221:20180830182732j:plain

　10期の自己相関に怪しい値が出ているが、 {\displaystyle p > 0.05 } となって有意な自己相関はないことからひとまずOKとする。同様に正規性も問題なさそうなので、今回はこのモデルを採用する。

モデルによる予測

　作成したモデルでテスト期間20期分について予測すると、こんな感じになる。

> # テストデータ期間の予測
> forecast_nile <- forecast(
+   model_nile,
+   h = 20,
+   level = c(95, 80)
+ )
> autoplot(forecast_nile) +
+   theme_bw()

f:id:yanbow221:20180831132507j:plain

　差分系列が定常過程であり、AR過程、MA過程ともに次数が1であることもあって、3期後あたりから横ばいの予測値になっている。定常過程は平均回帰性を持っており、差分系列が定常過程であり定数項がないことによって、差分は徐々に0に近づいていく。そのため原系列も一定値で推移するようになる。一方信頼区間は期が増えるに従い線形に増大しており、これは原系列が単位根を持っていることによる特徴である(定常過程では分散は一定値に収束していく)。この2点から、素朴なARIMAが長期的な予測には向かないと言える。非定常過程のその他の手法として状態空間モデル等あるが、ARIMAを基礎にした工夫としては例えば以下が参考になる。先に紹介した「時系列分析と状態空間モデルの基礎」の著者の方のブログで、時系列予測に関して各記事とても勉強になる。本記事を書くにあたっても書籍と併せて大いに参考にした。今回は分析の大まかな流れを追うため、以降の工夫は省略する。

logics-of-blue.com

　さて、次に予測の良し悪しを評価するが、そのためには比較対象が必要になる。上回るべき最初の比較対象としては「学習データの平均値」や「学習データの最終時点の値」などがあり、こうした原系列からの素朴な予測値を『ナイーブな予測』という(「時系列分析と状態空間モデルの基礎」P113)。評価指標としてはMSE(平均二乗誤差)がよく用いられる。パッケージforecastに入っているaccuracyを使えばこれの平方根をとったRMSE(平均平方二乗誤差)をはじめいくつかの指標を計算してくれる。可視化と併せてこんな感じ。

> plot(Nile_train, xlim = c(1871, 1970))
> lines(Nile_test, lty = "dashed")
> lines(forecast_nile$mean, col = "steelblue", lwd = 2)
> lines(naive_f_mean$mean, col = "purple", lwd = 2)
> lines(naive_f_latest$mean, col = "rosybrown", lwd = 2)
> legend("topright", 
+        legend = c("ARIMA(1,1,1)", "学習データ平均", "学習データ最終点"),
+        col = c("steelblue", "purple", "rosybrown"),
+        lty = "solid", lwd = 2)

f:id:yanbow221:20180831140332j:plain

> accuracy(forecast_nile, x = Nile_test)    #ARIMA(1,1)による予測の評価
                    ME     RMSE      MAE        MPE     MAPE      MASE        ACF1 Theil's U
Training set -18.24538 143.2728 111.4261 -4.3359424 12.97822 0.8310673 -0.03829202        NA
Test set      22.29061 125.0452 106.1632  0.6760674 11.96723 0.7918141  0.19619479  0.869674
> accuracy(naive_f_mean, x = Nile_test)     #ナイーブ予測①(学習データの平均)の評価
                        ME     RMSE      MAE       MPE     MAPE      MASE      ACF1 Theil's U
Training set  4.546749e-14 176.4544 148.3694 -3.840484 16.70686 1.1066069 0.5301332        NA
Test set     -5.287500e+01 133.3132 108.0125 -8.046549 13.15615 0.8056068 0.1922875 0.9354761
> accuracy(naive_f_latest, x = Nile_test)   #ナイーブ予測②(学習データの最終点)の評価
                     ME     RMSE      MAE       MPE     MAPE      MASE       ACF1 Theil's U
Training set  -2.911392 170.7412 134.0759 -2.198271 15.14585 1.0000000 -0.3925129        NA
Test set     -12.950000 123.0624 101.9500 -3.407725 11.96265 0.7603899  0.1922875 0.8623788

　RMSEについて、ナイーブ予測①よりは小さいもののナイーブ予測②よりも大きい、つまり素朴な予測よりもイマイチな予測になっていることが分かった。ここからは前述したような工夫もしつつ、精度向上の手を考えていくことになる。

4.終わりに

　今回は最もよく用いられるだろうARIMAによる時系列分析について概観した。その他季節性や外因性を考慮したモデルや多変量への拡張、状態空間モデル等、整理してそのうちまとめたい。

*1:AR、MA、ARMAを内包する

*2:株のことはよく知らないのでこのような判断があるのかどうかわからないです。完全にたとえとしての妄想

*3:通常の回帰では同一変量内の値の間に独立を仮定する

*4:時点tと時点t-kにおいて、その間の時点の影響を取り除いた自己相関のこと

*5:p時点前までの自分自身で回帰するのでp+1時点以降影響がなくなるのは直感的に分かりやすいと思う

*6:自己相関 {\displaystyle ρ_k = γ_k / γ_0 }

*7:関心を持っている時系列データそのもののこと

*8:明らかに関係のない変数間に有意な関係を示してしまうこと