自己回帰過程のサンプル

今回は移動平均過程の例をいくつか生成して紹介します。
$AR(1)$や$AR(2)$の例は色々なところで取り上げられているので、すべて$AR(3)$にしました。

AR特性方程式の解の絶対値が全て1より大きい定常な例を一つと、そうではない例が3つです。
定常ではない例は、AR特性方程式の解に絶対値が$1$より小さなものを含むものと、
$-1$を解に持つもの、$1$を解に持つものを選びました。
なお、撹乱項の分散は$1$に揃えています。
$$\varepsilon_t \sim iid N(0,1)$$

定常なAR(3)過程の例

まずは、次のAR(3)過程を取り扱います。
$$y_t=1+\frac{11}{15}y_{t-1}-\frac{1}{3}y_{t-2}+\frac{1}{15}y_{t-3}+\varepsilon_t$$
AR特性方程式の解は $3,1+2i,1-2i$の3個です。(意図的に虚数解を含むものを選びました。)
解の絶対値が全部1より大きいので、この過程は定常になるはずです。
pythonで時系列データを生成しプロットして見ます。


import sympy
import numpy as np
import matplotlib.pyplot as plt

phi_1 = 11/15
phi_2 = -1/3
phi_3 = 1/15

# 根の確認
z = sympy.symbols('z')
print(sympy.solve(1-phi_1*z-phi_2*z**2-phi_3*z**3))

# 過程の生成
y = np.zeros(200)
# 初期値
y[0] = np.random.randn()
y[1] = np.random.randn()
y[2] = np.random.randn()
for i in range(3, 200):
    y[i] = 1 + phi_1*y[i-1] + phi_2*y[i-2] + phi_3*y[i-3]  + np.random.randn()


fig = plt.figure(figsize=(12, 6))
title = '$y_t=1+\\frac{11}{15}y_{t-1}-\\frac{1}{3}y_{t-2}+\\frac{1}{15}y_{t-3}+\\varepsilon_t$'
ax = fig.add_subplot(1, 1, 1, title=title)
ax.plot(y)
plt.show()

出力がこちら。

いかにも定常な感じですね。

定常ではないAR(3)過程の例 1

次は定常ではないAR(3)過程の例としてこれを考えます。
$$y_t=1+\frac{1}{4}y_{t-1}+4y_{t-2}-y_{t-3}+\varepsilon_t$$
AR特性方程式の解は $-0.5,0.5,4$ であり、絶対値が1より小さいものを含みます。

コードは省略しますが、試しに出力したのがこちら。
撹乱項の値によって形が全然変わるので実行するたびに プラスに振れるのかマイナスに振れるのかも違います。

どう見ても定常ではありませんね。$t$が小さい時と大きい時の$y_t$の期待値が全く違います。

定常ではないAR(3)過程の例 2

もう一つ定常ではないAR(3)過程の例。
今度はAR特性方程式が$-1$を解に持ちます。
$$y_t=1-\frac{1}{6}y_{t-1}+\frac{2}{3}y_{t-2}-\frac{1}{6}y_{t-3}+\varepsilon_t$$

プロットしたのがこちら。

一見定常っぽくも見えますが、非常にギザギザしています。
実はこの過程、奇数番目と偶数番目の項の期待値が違います。
(撹乱項を取り払うと2つの値が交互に現れるので確認しやすいです)
このため、たしかに定常でないことがわかりました。

定常ではないAR(3)過程の例 3

次も定常ではない例ですが、今回はAR特性方程式が、$1$を解に持ちます。
$$y_t=1+\frac{11}{6}y_{t-1}-y_{t-2}+\frac{1}{6}y_{t-3}+\varepsilon_t$$
実はこの過程は定常ではないのですが、差分系列$\Delta y_t=y_t-y_{t-1}$は定常過程になります。
このような過程を単位根過程と言います。
単位根過程の詳しい話はそのうち紹介したいですがとりあえずグラフを出すと次のようになります。

綺麗なトレンドが出ましたね。

自己回帰過程の定義と定常になる条件、定常自己回帰過程の性質

例によって沖本本を参照し、自己回帰過程を紹介します。

自己回帰(AR)過程(autoregressive process)は、過程が自分の過去に回帰された形で表現される過程です。
p次AR過程(AR(p)過程)は、次で定義されます。
$$
y_t = c + \phi_1y_{t-1} + \phi_2y_{t-2} + \cdots + \phi_py_{t-p} + \varepsilon_t ,\ \
\varepsilon_t \sim W.N.(\sigma^2)
$$

要するに$y_t$を過去p期間の値に回帰したモデルです。
AR(p)過程は常に定常になるとは限りません。

AR(p)過程が定常になるかどうかは次の方法で確認することができます。
まず、次の変数zの方程式を考えます。 ($\phi_t$は上の定義で登場した係数)
$$
1-\phi_1z-\phi_2z^2-\cdots-\phi_pz^p = 0
$$
これはAR特性方程式とよばれます。また、左辺の多項式はAR多項式と呼ばれるそうです。
この方程式の全ての解の絶対値が1より大きい時、AR過程は定常になります。

AR過程が定常な場合、下記の性質が成り立ちます。
(改めてですが、AR過程は定常でないこともあるので注意です。)

$$
\begin{eqnarray}
\mu & = E(y_t) = \frac{c}{1-\phi_1-\phi_2-\cdots-\phi_p}\\
\gamma_0 & = Var(y_t) = \frac{\sigma^2}{1-\phi_1\rho_1-\phi_2\rho_2-\cdots-\phi_p\rho_p}\\
\gamma_k & = \phi_1\gamma_{k-1}+\phi_2\gamma_{k-2}+\cdots+\phi_p\gamma_{k-p}, k\geq1\\
\rho_k & = \phi_1\rho_{k-1}+\phi_2\rho_{k-2}+\cdots+\phi_p\rho_{k-p}, k\geq1
\end{eqnarray}
$$
そして、AR過程の自己相関は指数的に減衰します。

移動平均過程の定義と性質

今回も出典は沖本先生の経済・ファイナンスデータの計量時系列分析から。

ホワイトノイズの線形和で表される過程を移動平均過程(moving average process)と呼びます。
$q$を1以上の整数ととすると、q次移動平均過程$MA(q)$は次で定義されます。
$$y_t=\mu+\varepsilon_t+\theta_1\varepsilon_{t-1}+\theta_2\varepsilon_{t-2}+\cdots+\theta_q\varepsilon_{t-q},\ \ \varepsilon_t\sim W.N.(\sigma^2)$$

そして、$y_t$が$MA(q)$過程に従うことを $y_t\sim MA(q)$と書きます。
$MA(q)$過程は常に定常であり、さらに次の性質を持ちます。

$$
\begin{align}
E(y_t) &= \mu\\
\gamma_0 &= Var(y_t) = (1+\theta_1^2++\theta_2^2+\cdots+\theta_q^2)\sigma^2\\
\gamma_k & = \left\{
\begin{matrix}
\ (\theta_k+\theta_1\theta_{k+1}+\cdots+\theta_{q-k}\theta_{q})\sigma^2, & (1\leq k \leq q)\\
0,& (k\geq q+1)
\end{matrix}
\right.\\
\rho_k & = \left\{
\begin{matrix}
\frac{\theta_k+\theta_1\theta_{k+1}+\cdots+\theta_{q-k}\theta_{q}}
{1+\theta_1^2++\theta_2^2+\cdots+\theta_q^2} , & (1\leq k \leq q)\\
0,& (k\geq q+1)
\end{matrix}
\right.
\end{align}
$$
式から分かる通り、MA(q)過程の q+1次以降の自己相関は常に0になります。
そのため、長期間にわたる自己相関を移動平均過程でモデル化するには多くのパラメーターが必要になります。

カバン検定を自分で実装してみる

前の記事で、statsmodels に実装されている acorr_ljungbox 関数を使って、カバン検定を行ってみました。
statsmodelsでかばん検定 (自己相関の検定)

今回は学習のため、numpyで実装しました。
実装にあったって数式から説明します。
検定の帰無仮説や対立仮説については前の記事を参照してください。

$\{y_t\}_{0}^{T-1}$を時系列データとします。 (沖本本の定義とインデックスが一つずれているので注意。)

標本平均$\bar{y}$, 標本自己共分散$\hat{\gamma}_k$, 標本自己相関係数$\hat{\rho}_k$を次のように定義します。
\begin{eqnarray}
\bar{y} & = & \frac{1}{T}\sum_{t=0}^{T-1}y_t\\
\hat{\gamma}_k & = & \frac{1}{T}\sum_{t=k}^{T-1}(y_t-\bar{y})(y_{t-k}-\bar{y}) , \hspace{1em} k = 0, 1, 2, \cdots\\
\hat{\rho}_k & = & \frac{\hat{\gamma}_k}{\hat{\gamma}_0},\hspace{1em} k = 1, 2, 3, \cdots
\end{eqnarray}

この時、Ljung と Box が考案した次の統計量$Q(m)$は漸近的に自由度$m$のカイ2乗分布に従います。
$$Q(m) = T(T+2)\sum_{k=1}^{m}\frac{\hat{\rho}_k^2}{T-k}$$

今回も7点周期のデータを準備します。
m = 7 で検定を行いますので、帰無仮説が棄却されれば、
LAG1〜7のいずれかの自己相関係数が0でないことが主張されます。
(p値は5%を基準に判断します。)


import numpy as np
import pandas as pd
from scipy.stats import chi2
# 答え合わせ用
from statsmodels.stats.diagnostic import acorr_ljungbox

# 7点ごとに周期性のあるデータを準備
series = pd.Series([1, 1, 1, 1, 1, 1, 5]*10)
# 乱数追加
series += np.random.randn(70)

# データの個数
T = len(series)
# 検定対象のm
m = 7
# 標本平均
y_ = series.mean()
# 標本自己共分散の列
gamma_list = np.array([
                np.sum([
                    (series[t]-y_)*(series[t-k]-y_) for t in range(k, T)
                ])/T for k in range(m+1)
            ])
# 標本自己相関係数の列
ro_list = gamma_list/gamma_list[0]

# Q(m)の計算
Q = T*(T+2)*np.sum([ro_list[k]**2 / (T-k) for k in range(1, m+1)])
print("lb値:", Q)
print("p値:", chi2.sf(Q, m))

# ライブラリを使った計算結果
lbvalues, pvalues = acorr_ljungbox(series, lags=7)
print(lbvalues[6])
print(pvalues[6])

# 出力結果
lb値: 43.78077348272421
p値: 2.3564368405873674e-07
43.780773482724186
2.3564368405873896e-07

p値が非常に小さな値であることがわかりました。
また、ライブラリを使った結果と四捨五入による誤差がありますが、それ以外は同じ値なのでミスも無さそうです。

この計算の注意点としては、
標本自己相関係数を使うところでしょうか。
pandasのautocorrは自己相関係数を出力してくれますが、
これが標本自己相関係数とは微妙に定義が異なり、autocorrを使うと少し値がずれてしまいます。

statsmodelsでかばん検定 (自己相関の検定)

時系列データを分析をするとき、そのデータが自己相関を持つかどうかはとても重要です。
データが自己相関を持っていたらその構造を記述できるモデルを構築して、予測等に使えるからです。
逆に自己相関を持っていないと、時系列分析でできることはかなり限られます。
過去のデータが将来のデータと関係ないわけですから当然ですね。
(どちらも沖本本の 1.4 自己相関の検定から)

ということで今回行うのは自己相関の検定です。

自己相関が全てゼロという帰無仮説、つまり
$H_0:\rho_1=\rho_2=\cdots=\rho_m=0$を、
$H_1:$少なくとも1つの$k\in[1,m]$において、$\rho_k\neq0$
という対立仮説に対して検定します。

この検定はかばん検定(portmanteau test)と呼ばれているそうです。
検定量は色々考案されているそうですが、 Ljung and Box が考案されたものがメジャーとのこと。

具体的な数式や、numpyでの計算例は次回に譲るとして、とりあえずpythonのライブラリでやってみましょう。
statsmodels に acorr_ljungbox という関数が用意されています。

statsmodels.stats.diagnostic.acorr_ljungbox
(完全に余談ですが、statsmodelsのドキュメントで portmanteau という関数を探していたので、これを見つけるのに結構苦労しました)

あからさまに7点周期を持つデータを準備し、1から10までのmに対して検定を実施したコードがこちらです、
acorr_ljungbox は lb値と、p値をそれぞれのmに対して返します。


import pandas as pd
import numpy as np
from statsmodels.stats.diagnostic import acorr_ljungbox

# 7点ごとに周期性のあるデータを準備
series = pd.Series([1, 1, 1, 1, 1, 1, 5]*10)
# 乱数追加
series += np.random.randn(70)

lbvalues, pvalues = acorr_ljungbox(series, lags=10)
lag = 1
for lb, p in zip(lbvalues, pvalues):
    print(lag, lb, p)
    lag += 1

# 出力
1 4.095089120802025 0.04300796255246615
2 4.223295881654548 0.1210383379718042
3 6.047027807671336 0.10934450247698609
4 6.312955422660912 0.1769638685520115
5 6.457291061922424 0.26422887039323306
6 9.36446186985462 0.15409458108816818
7 40.47102807634717 1.0226763047711032e-06
8 45.2406378234939 3.312995830103468e-07
9 45.24892593869829 8.298135787344997e-07
10 46.35530787049628 1.2365386593885822e-06

lag が7未満の時は、p値が0.05を超えているので、帰無仮説$H_0$を棄却できず、
7以上の時は棄却できていることがわかります。

念の為ですが、lag が 8,9,10の時に棄却できているのは、
どれもデータが7点周期を持っていることが理由であり、
8,9,10点の周期性を持っていること意味するものではありません。

定常過程の例としてのiid過程とホワイトノイズ

前回の記事で定義した定常過程の例を紹介します。
(今回も沖本本からの引用です。)

iid系列

iid系列は、強定常過程の例です。

定義:
各時点のデータが互いに独立でかつ同一の分布に従う系列はiid系列と呼ばれる。

iidとは、 independent and identically distributedの略です。
Wikipedia : 同時独立分布

$y_t$が期待値$\mu$,分散$\sigma^2$のiid系列であることを $y_t\sim iid(\mu,\sigma^2)$と書きます。

ホワイトノイズ

独立性や同一分布性は非常に強い仮定なのですが、
実用上、モデルの錯乱項などとして使うのであればもう少し弱い仮定で十分です。
そこで次のホワイトノイズというものが使われます。
ホワイトノイズは弱定常過程です。

定義:
すべての時点$t$において
$$
\begin{align}
E(\varepsilon _t) &=0\\
\gamma _k &= E(\varepsilon _t \varepsilon _{t-k}) =\left\{
\begin{array}{ll}
\sigma^2, & (k=0)\\
0, & (k\neq 0)
\end{array}\right.
\end{align}
$$
が成立するとき, $\varepsilon _t$はホワイトノイズ(white noise)と呼ばれる.

$\varepsilon _t$が分散$\sigma^2$のホワイトノイズであることを $\varepsilon _t\sim W.N.(\sigma^2)$と書きます。
次のようにホワイトノイズを使って、基礎的な弱定常過程の例を作れます。
$$
y_t = \mu + \varepsilon _t,\ \ \ \varepsilon _t\sim W.N.(\sigma^2)
$$

自分だけかもしれないのですが、ホワイトノイズと聞くと無意識に$(\varepsilon _t,\cdots,\varepsilon _{t+k})$の同時分布が期待値0の多変量正規分布であるかのようにイメージしてしまいます。
しかし定義の通り、ホワイトノイズにはそこまでの強い仮定はないので注意です。

ちなみに $(y_t,\cdots,y_{t+k})$の同時分布が多変量正規分布の場合、それは正規過程と呼ばれ、
正規過程に従うホワイトノイズは正規ホワイトノイズと呼ばれます。(沖本本のp9,p12)
正規ホワイトノイズはiid過程なので、強定常過程になります。

確率過程の弱定常性と強定常性

時系列データの分析において重要な定常性の定義のメモです。
用語は沖本先生の経済・ファイナンスデータの計量時系列分析の定義に従っています。

確率過程の定常性とは、その過程の同時分布や基本統計量の時間不変性に関するものです。
何を不変とするかによって、弱定常性 (weak stationarity)と強定常性(strict stationarity)に分類されます。

以下、$y_t$を確率過程とします。

弱定常性

弱定常性は、過程の期待値と自己共分散が一定であることです。

定義:
任意の$t$と$k$に対して、
$$
\begin{align}
E(y_t) & = \mu\\
Cov(y_t, y_{t-k}) & = E[(y_t-\mu)(y_{t-k}-\mu)] = \gamma_k
\end{align}
$$
が成立する場合、過程は弱定常(weak stationary)と言われます。

$y_t$の期待値は任意の$t$に対して一定で、$y_t$と$y_{t-k}$の共分散は$k$の値によってのみ決まります。
また、この二つが時点に依存しないので、自己相関も時点に依存しません。

経済・ファイナンスの分野では単に定常性というと、弱定常性をさすことが多いそうです。(沖本本 P10)

強定常性

強定常性は、同時分布が不変であることを要求します。

定義:
任意の$t$と$k$に対して、$(y_t,y_{t+1},\cdots,y_{t+k})$の同時分布が同一となる場合、
過程は強定常(strict stationary)と言われます。

弱定常は期待値と共分散だけに対する仮定だったので、強定常のほうが名前の通り強い仮定になります。

ただし、強定常ならば絶対に弱定常というわけではなく、
“過程の分散が有限であるならば、”強定常過程は弱定常過程である、
という風に一つ条件が必要です。

pandasで時系列データの自己相関係数算出

時系列データとそれをn時間分シフトさせたデータ間の相関係数を自己相関係数と呼びます。
データに周期性があるかどうかを確認するときなどに調査します。

単に配列をスライドさせてnumpyか何かで相関係数を出してもいいのですが、
pandasに専用の関数が定義されているのでそれを使って求めることもできます。

pandas.Series.autocorr

試しにやってみましょう。


import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

# 7点ごとに周期性のあるデータを準備
series = pd.Series([5, 4, 5, 6, 5, 3, 1]*10)
# 乱数追加
series += np.random.randn(70) * 0.5

# lag=0から29までの自己相関係数
auto_series = [series.autocorr(lag=i) for i in range(30)]

# 可視化
fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(2, 1, 1, title="周期性のあるデータ")
ax.plot(series.index, series)
ax = fig.add_subplot(2, 1, 2, title="自己相関係数")
ax.bar(range(30), auto_series)
plt.show()

出力されるグラフがこちら。

7点間隔の周期が綺麗に現れているのがわかります。

LaTeXで行列の成分を省略したときのドットを書く

大学院生のときは頻繁に書いていたのに、すっかり忘れてしまっていたので$\LaTeX$の復習です。

まず、要素を省略せずに3*3行列を書く例。

$$
\textbf{A} =
\left(
\begin{array}{ccc}
a_{11} & a_{12} & a_{13} \\
a_{21} & a_{22} & a_{23} \\
a_{31} & a_{32} & a_{33}
\end{array}\right)
$$

結果。
$$
\textbf{A} =
\left(
\begin{array}{ccc}
a_{11} & a_{12} & a_{13} \\
a_{21} & a_{22} & a_{23} \\
a_{31} & a_{32} & a_{33}
\end{array}\right)
$$

次に、ドットを使って書く例。
横並びは\cdots, 縦並びは\vdots, 斜めは\ddots を使います。


$$
\textbf{W} =
\left(
\begin{array}{ccc}
w_{11} & \cdots & w_{1n} \\
\vdots & \ddots & \vdots \\
w_{n1} & \cdots & w_{nn}
\end{array}\right)
$$

結果。
$$
\textbf{W} =
\left(
\begin{array}{ccc}
w_{11} & \cdots & w_{1n} \\
\vdots & \ddots & \vdots \\
w_{n1} & \cdots & w_{nn}
\end{array}\right)
$$

scipyで手軽にカイ二乗検定

業務でカイ二乗検定を行う場面は、割と多くあります。
自分の場合は、ABテストの評価で使うことが多いです。

Excelでも自分でコードを書いてでもさっとできるのですが、
scipyに専用の関数が用意されているのでその紹介です。

scipy.stats.chi2_contingency

これを使うと、引数にデータのクロス表を放り込むだけで、
カイ二乗値、p値、自由度、期待度数表を一発で算出してくれます。
3つの値と、期待度数がarray形式ではいったタプルが戻るので、引数を4つ用意して受け取ると便利です。


import pandas as pd
from scipy.stats import chi2_contingency

# 検定対象のデータ(サンプルなので値はダミー)
df = pd.DataFrame([[300, 100], [800, 200]])

chi2, p, dof, expected = chi2_contingency(df, correction=False)

print("カイ二乗値:", chi2)
print("p値:", p)
print("自由度:", dof)
print("期待度数:", expected)

# 以下出力
カイ二乗値: 4.242424242424244
p値: 0.03942583262944296
自由度: 1
期待度数: [[314.28571429  85.71428571]
 [785.71428571 214.28571429]]

このようにめんどな手間がなく一瞬で検定ができました。
ちなみに、
correction=False
はイエーツの修正を”行わない”設定です。
デフォルトがTrueなので注意してください。