はじめに
こんにちは！
僕は研究を行いながら、長期インターンでデータサイエンティストとして働く大学院生です！
学部時代から長期インターンを始め、現在まで４社経験してきました。
この経験から、プログラミングの学習を始めたばかりの人や、長期インターンを行う勇気が出ない人に、学習サポートやデータ分析の実績作り支援、さらにはKaggle人材マッチングサービスなどを行わせてもらっています！
僕自身、プログラミングの習得や長期インターン探しに苦労したので、その経験をお伝えすることで、より多くの人が挫折せずデータサイエンティストになるまで成長して欲しいです！
以下でサポートを行なっているのでご興味ある方はご連絡ください！学生・社会人問わず専攻も問わずサポートいたします！
X(Twitter)
今回は、前回の「統計学実践ワークブック完全攻略5章〜8章」に引き続き、27章の時系列解析を解説します。
統計学実践ワークブック完全攻略5章~8章
他にもこの章を解説してほしいというのがあれば、XのDMにてお問合せください！
27章 時系列解析
この章では以下の重要な概念について詳しく解説していきます。また、ワークブックの内容に加えて、Prophetという時系列モデルについても解説します。
- 自己共分散
- 定常性
- 自己回帰過程
- 移動平均過程
- 自己回帰移動平均過程
- 次数の決定法
- 状態空間モデル
- ダービン・ワトソン検定
- Prophet
自己共分散
期間$t$における時系列データを$Y_t$とします。
一般に、時系列データにおける平均や分散は、観測時点$t$に依存するので、以下で表されます。
\displaylines{
E[Y_t]=μ_t\\
V[Y_t]=σ^2_t
}
また、時系列データは、同一変数の異なる時点($Y_t$と$Y_{t-1}$など)の共分散や相関係数で特徴づけられることがしばしばあります。これは直感的にもそうで、例えば、サラリーマン家庭の所得は、残業代などで前後しますが、前月と今月で所得額を比較すると、似た金額であることが多いです。このような場合は、$Y_t$と$Y_{t-1}$は正の相関となります。同一データの異なる時点の共分散を自己共分散、相関係数を自己相関係数と呼びます。
式で表すと以下になります。
\displaylines{
Cov[Y_t, Y_{t-h}]=γ_{t,h}\\
ρ[Y_t,Y_{t-h}]=\frac{Cov[Y_t,Y_{t-h}]}{\sqrt{V[Y_t]V[Y_{t-h}]}}=ρ_{t,h}
}
これはそれぞれ、時点$t$における$h$次の自己共分散、時点$t$における$h$次の自己相関係数です。
また、時間差$h$はラグと呼ばれています。
定常性
時系列データの平均と分散が有限で、平均が観測時点$t$には依存せず、また、自己共分散も観測時点$t$には依存せずに時間差$h$にのみ依存する場合、その系列は共分散定常過程と言います。ある値を中心に上下に動く時系列データなどがそれに当たります。
すなわち${Y_t}$が共分散定常過程であるとき、以下で表せます。
\displaylines{
E[Y_t]=μ\\
Cov[Y_t,Y_{t-h}]=γ_{|h|}
}
共分散定常過程の平均が0, $h\ne0$のすべての自己共分散が0である時、その系列をホワイトノイズと言います。ホワイトノイズは、時系列データを用いた回帰分析では、しばしば、誤差項が満たすべき仮定として想定されます。
たとえば以下の時系列重回帰モデルの$U_t$はホワイトノイズと仮定されることが多いです。
Y_t=β_0+β_1X_{1t}+β_2X_{2t}+...+β_KX_{Kt}+U_t
自己回帰過程
自己回帰過程(Autoregressive process)は時系列分析で広く用いられる統計モデルの一つです。このモデルは、現在の時系列値が過去のいくつかの時系列値の線形組み合わせとホワイトノイズによって生成されると仮定します。ARモデルは、特に経済学、金融学、気象学などの分野で時間に依存するデータをモデル化するためによく使用されます。
たとえば、家計の年間消費支出などはARモデルで説明され得ます。多くの家計では、年々の諸遠くはそれほど大きくは変化せず、そのため、消費行動も比較的似たものになる傾向があり、過去の消費額と現在の消費額は正の値を取りやすいからです。
最もシンプルな1次の自己回帰過程(AR(1)過程)であり、以下のように表されます。
Y_t=φ_1Y_{t-1}+U_t (t=1,2,...,T)\tag{1}
ここで${U_t}$は系列${Y_t}$を生成するショックと呼ばれ、ホワイトノイズと仮定されます。$φ_1$は未知パラメータであり、この値により${Y_t}$の性質が変わってきます。$(1)$を逐次的に過去へ向かって解くと、次のように表現できます。
\displaylines{
Y_t=φ_1Y_{t-1}+U_t\\
=φ_1(φ_1Y_{t-2}+U_{t-1})+U_t\\\
...\\
=φ_1^tY_0+\sum_{j=1}^tφ_1^{t-j}U_j
}
この式より、$Y_0$から$Y_t$への影響は$φ^t_1$であることがわかります。従って、$|φ|<1$であれば、過去から現在への影響は、時間差が開くほど弱くなることになります。一方で$|φ|=1$であれば過去からの影響は時間差に関係なく一定であり、$|φ|>1$であれば、過去に遡れば遡るほど、現在への影響は大きくなります。
AR(1)過程が共分散定常である時は以下のことが言えます。ここでは、AR(1)過程をより一般的な式で示します。
\displaylines{
Y_t=c+φ_1Y_{t-1}+U_t\\
E[Y_t]=E[c+φ_1Y_{t-1}+U_t]\\
=c+φ_1E[Y_{t-1}]\\
=c+φ_1E[Y_t](平均はtに依らず一定)\\
E[Y_t]=\frac{c}{1-φ_1}\tag{2}
}
AR(1)過程において、一般に$h$次の自己共分散については、$h=0,1,2,...$と逐次的に解くことにより、以下の式が得られます。
γ_h=φ_1^h\frac{σ^2}{1-φ_1^2}
また、$h$次の自己相関係数は以下で表されます。
ρ_h=φ_1^h
より一般的に$p$次のAR過程(AR(p))を考えてみます。式$(2)$と同様に考えれば、$Y_t$の平均は以下のように求まります。
\displaylines{
Y_t=c+φ_1Y_{t-1}+φ_2Y_{t-2}+...+φ_pY_{t-p}+U_t\\
μ=\frac{c}{1-φ_1-,...,-φ_p}
}
移動平均過程
移動平均過程は、AR過程とおこなり、$Y_t$が現在と過去のショックの加重和で生成されている系列です。具体的に1次の移動平均過程(MA(1))を以下に示します。
Y_t=μ+U_t+θ_1U_{t-1}
MA(1)の自己共分散は以下のように表されます。
γ_h =
\begin{cases}
(1+θ_1^2)σ^2 & \text{if } h = 0, \\
θ_1σ^2 & \text{if } h = 1, \\
0 & \text{else }
\end{cases}
式より、MA(1)は自己共分散がhに依らないので、共分散定常であると言えます。
これを一般的なMA(q)の場合にも拡張して考えると以下のように表せます。
\displaylines{
Y_t=μ+U_t+θ_1U_{t-1}+θ_2U_{t-2}+...+θ_qU_{t-q}\\
γ_h =
\begin{cases}
(1+θ_1^2+...+θ_q^2)σ^2 & \text{if } h = 0, \\
(θ_h+θ_1θ_{h+1}+...+θ_{q-h}θ_q)σ^2 & \text{if } 1 <= h <= q, \\
0 & \text{if } h>q
\end{cases}
}
こちらも$h$によらないため、共分散定常であることがわかります。
自己回帰移動平均過程
こちらはAR過程とMA過程を組み合わせたより複雑な過程です。式としては以下で表せます。
Y_t=c+φ_1Y_{t-1}+...+φ_pY_{t-p}+U_t+θ_1U_{t-1}+...+θ_qU_{t-q}
上記式は、(p,q)次の自己回帰移動平均過程ARMA(p,q)過程と呼ばれます。
ARMA(p,q)過程は比較的低次でも複雑な自己共分散構造を表現できます。
一方で、時系列データによっては、${Y_t}$が共分散定常ではなく、階差をとった${ΔY_t}$が強分散定常であるとみなされる場合もあります。ここでは、$ΔY_t=Y_t-Y_{t-1}$と定義されています。
つまり、時系列が時間によってその平均や分散が大きく変化しない場合、ARMAモデルが適切な可能性が高く、時系列がトレンド（時間とともに増加または減少する傾向）や季節性などの非定常性を示す場合ha,
ARIMAモデルが適している可能性が高いです。
式としては以下で表されます。
ΔY_t=φ_1ΔY_{t-1}+...+φ_pΔY_{t-p}+U_t+θ_1U_{t-1}+...+θ_qU_{t-q}
次数の決定法
ここまで様々な統計モデルを解説してきましたが、結局時系列データに対してどのモデルが適しているのかを検討することが大切です。その方法について説明します。
それでは、ある時系列データに対して、AR(1),AR(2),MA(1),MA(2),ARMA(1,1)過程のどれが、最も当てはまりが良いかを推定することを考えます。自己相関係数を考えると、AR(1),AR(2),ARMA(1,1)に関しては、徐々に0に近づいていきますが、MA(1),MA(2)に関しては、２次または3次以上の自己相関係数が0になります。
従って、得られたデータから標本自己相関係数を求めて2次もしくは3次以上で有意に0と異なるかどうかでMA過程を選択するか判断できます。
なお、横軸に次数、縦軸に自己相関係数をとった図をコレログラムと言います。
では、残ったAR系のモデルからどのように次数を決定すれば良いのでしょうか。そこで登場するのが、偏自己相関係数です。$h$次の偏自己相関係数とは、$Y_t$から$Y_{t-1},...,Y_{t-h+1}$の影響を取り除いたもの(回帰の残差)と$Y_{t-h}$から$Y_{t-1},...,Y_{t-h+1}$の影響を取り除いたもの(回帰の残差)との間の相関係数です。
この値を算出するためには、非常に複雑な計算が必要ですが、PythonやRなどの統計ソフトウェアを使えば、簡単に算出してくれます。
AR(p)過程の偏自己相関係数は、$p+1$次以上は全て0になりますが、MA(q)過程やARMA(p,q)過程の偏自己相関係数は、徐々に減衰することがわかっています。
従って、標本偏自己相関係数を求めて2次もしくは3次以上で有意に0と異なるかどうかで、AR過程稼働を選択する判断ができます。
以上をまとめると以下のようになります。
ただし、実際のデータを用いて標本自己相関係数や標本編自己相関係数を求めても、様々な次数で0から有意に離れたり離れなかったりするため、これだけで次数の選択を行うのは難しい場合もあります。
そのような場合には、AICやBICといった評価指標を用いて次数を選択できます。この場合は、ARの最高次数とMAの最高次数をあらかじめ設定して、それを上回らない範囲で全ての次数の組み合わせでモデルを推定して、最良の結果となる次数を選択すれば良いです。
状態空間モデル
ARIMA過程よりも複雑な時系列モデルの例として、状態空間モデルがあります。式としては次のように表せます。
Y_t=X_tβ_t+U_t\tag{1}
β_t=Φβ_{t-1}+V_t\tag{2}
ここで、$X_t$は$K$次元の説明変数であり、$β_t$は$K$次元回帰係数で時点$t$に依存し変動するパラメータです。式$(1)$は観測方程式、式$(2)$は状態方程式と呼びます。この２つの式を合わせて状態空間モデルと呼びます。
$β_t$の推定には、カルマンフィルターという手法が用いられます。
状態空間モデルを用いれば、様々な時系列モデルの表現が可能となります。可変パラメータモデル以外にも、観測されない潜在変数を織り込んだ次のモデルを考えることもあります。
\displaylines{
Y_t=α+βZ_t+U_t\\
Z_t=φZ_{t-1}+V_t
}
ダービン・ワトソン検定
以下のような時系列回帰モデルを考えます。
Y_t=β_1+β_2X_{2t}+...+β_KX_{Kt}+U_t\tag{3}
こちらの式は、線形回帰モデルであるため、最小二乗法で推定できます。ただし、最小二乗法で推定するためには、${U_t}$がホワイトノイズである必要があります。${U_t}$がホワイトノイズではなく、自己共分散が0ではない場合は、$U_t$の自己相関構造を考慮した最尤法や一般最小二乗法を用いた方が効率よく推定できます。
自己共分散が0でない場合、誤差項には自己相関もしくは、系列相関があると言われています。
誤差項に自己相関があるかどうかを仮説検定方法として、ダービン・ワトソン検定が行われることがあります。ダービン・ワトソン検定では、誤差項にAR(1)モデルを想定し、以下のように表します。
U_t=ρU_{t-1}+ε_t
ここで、$ε_t$はホワイトノイズです。$U_t$は、$ρ=0$の場合にホワイトノイズとなるので、ダービン・ワトソン検定では、帰無仮説と対立仮説は以下となります。
\displaylines{
帰無仮説 H_0:ρ=0\\
対立仮説 H_1:ρ>0
}
まず、式$(3)$を最小二乗法で推定し、得られた残差を$\hat{U_t}$とします。この残差を用いて、**ダービン・ワトソン比(DW比)**を検定統計量として用います。
DW=\frac{\sum_{t=2}^T(\hat{U_t}-\hat{U_{t-1}})^2}{\sum_{t=1}^T\hat{U_t}^2}
Tがある程度大きければ、以下のように近似できます。
DW=2(1-\hat{γ_1})
ただし、以下が成り立ちます。
\hat{γ_1}=\frac{\sum_{t=2}^T{\hat{U_t}\hat{U}_{t-1}}}{\sum_{t=1}^T\hat{U_t}^2}
$\hat{γ_1}$は1次の自己相関係数の推定量とみなすことができるので、帰無仮説の元では0に近く、対立仮説の元では、正の値を取ると予想され、以下の関係性が見えてきます。
\displaylines{
ρ=0⇆\hat{γ_1}≒0⇆DW≒2\\
ρ>0⇆\hat{γ_1}>0⇆DW<2
}
従って、DW比の値が2に近い場合に帰無仮説を受容し、2より十分小さな場合に帰無仮説を棄却すれば良いです。
Prophet
ProphetはFacebookによって開発された時系列予測モデルです。特に、日々または年次の観測値を持つ時系列データに対して、効果的な予測を行うことを目的として設計されています。このモデルは、柔軟性が高く、さまざまな時系列データに適応しやすいのが特徴です。Prophetのアルゴリズムには「トレンド」「季節性」「祝日効果」「エラー項」の４つの主要なコンポーネントが含まれています。
Prophetはその柔軟性とユーザーフレンドリーさから、商業的なビジネス予測から学術研究まで、幅広い分野で利用されています。特に、日常的なビジネスメトリクスの予測や、欠損データや大規模な時系列データセットを扱う場合に有効です。
状態空間モデルとProphetを使い分ける際に検討することとしては、「データの特性」と「モデルの目的」が挙げられます。
データの特性に関しては、Prophetは明確なトレンドと季節性を持つ時系列に適しています。一方、状態空間モデルはより複雑なダイナミクスや非定常性を持つデータに対して柔軟に対応できます。
モデルの目的に関しては、解釈しやすく、実装が簡単なモデルを求めている場合はProphetが適しています。一方、より詳細な制御やカスタマイズを必要とする複雑な時系列を扱う場合は状態空間モデルを検討します。
例題
