定常自己回帰過程AR(3)の自己相関係数を具体的に求めてみる

自己回帰過程の定義と定常になる条件、定常自己回帰過程の性質
という記事で、定常自己回帰過程の性質を紹介しました。
その中に、逐次的に自己相関係数を算出する次の公式があります。

$$
\rho_k = \phi_1\rho_{k-1}+\phi_2\rho_{k-2}+\cdots+\phi_p\rho_{k-p}, k\geq1
$$
これをユール・ウォーカー方程式と言います。
こいつを使って、具体的に自己相関過程の自己相関係数を求めてみようというのが今回の記事です。

$p=1$や$p=2$の例は本に載っているので、$p=3$でやってみます。
$$
\rho_k = \phi_1\rho_{k-1}+\phi_2\rho_{k-2}+\phi_p\rho_{k-3}, k\geq1
$$
また、例として扱う定常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$$

$\phi_1,\phi_2,\phi_3$の値がわかるので、あとは、$\rho_0,\rho_1,\rho_2$さえわかれば残りの$\rho_k$は順番に計算できます。

これは、次の4つの関係式を使って計算できます。最初の2つは、ユール・ウォーカー方程式に$k=1,2$を代入したもの、
残り2つは自己相関係数の性質です。

$$
\rho_1 = \phi_1\rho_0 + \phi_2\rho_{-1} + \phi_3\rho_{-2}\\
\rho_2 = \phi_1\rho_1 + \phi_2\rho_0 + \phi_3\rho_{-1}\\
\rho_0 = 1\\
\rho_{k} = \rho_{-k}
$$
これを整理すると、
$$
\rho_1 = \phi_1 + \phi_2\rho_1 + \phi_3\rho_2\\
\rho_2 = \phi_1\rho_1 + \phi_2 + \phi_3\rho_1\\
$$
となり、さらに
$$
(1-\phi_2)\rho_1 – \phi_3\rho_2= \phi_1\\
(\phi_1+\phi_3)\rho_1 -\rho_2 =-\phi_2\\
$$
と変形できます。
あとは$\phi_k$の値を代入して、
$$
\frac{4}{3}\rho_1 – \frac{1}{15}\rho_2= \frac{11}{15}\\
\frac{4}{5}\rho_1 -\rho_2 =\frac{1}{3}\\
$$
となり、これを解くと次の二つが得られます。
$$\rho_1=\frac{5}{9},\rho_2=\frac19$$

あとは順番に代入するだけなので、pythonにやってもらいましょう。


import numpy as np
phi_1 = 11/15
phi_2 = -1/3
phi_3 = 1/15
rho_values = np.zeros(6)
rho_values[0] = 1
rho_values[1] = 5/9
rho_values[2] = 1/9
for i in range(3, 6):
    rho_values[i] = phi_1*rho_values[i-1]\
                                + phi_2*rho_values[i-2]\
                                + phi_3*rho_values[i-3]
for rho in rho_values:
    print(rho)

# 以下出力
1.0
0.5555555555555556
0.1111111111111111
-0.037037037037037035
-0.027160493827160487
-0.0001646090534979357

これで計算できました。
自己相関が指数関数的に減衰している様子も確認できます。

pythonでアルファベットの大文字小文字変換

自然言語処理系の機械学習を行うときやテキストマイニングをやるとき、
前処理としてアルファベットの大文字小文字を統一することがよくあります。
(失われる情報もあるのでしない方が良いという主張も見たことがありますが、僕は大抵の場合小文字に統一します。)

これはpythonの文字列が持っている
str.lower()str.upper()
を使うことで実現できます。
pythonの組み込み型のドキュメントを見ると乗っています。


text = "Hello World!"
print(text.lower())  # hello world!
print(text.upper())  # HELLO WORLD!

これだけだと、記事にしなかったと思うのですが、ドキュメントを読んでいると他にも関数が準備されていることがわかりました。


text = "Hello World!"
# 大文字と小文字を入れ替える
print(text.swapcase())  # hELLO wORLD!

text = "HELLO world!"
# 各単語の1文字目を大文字に、残りを小文字に変換する
print(text.title())  # Hello World!

text = "HELLO world!"
# 最初の文字を大文字に、残りを小文字に変換する
print(text.capitalize())  # Hello world!

正直、使う場面を思いつかないのですが面白いですね。
英語ネイティブな人たちが使うのでしょうか。

pythonで文字と文字コードの相互変換

pythonで文字を文字コードに変えたり、文字コード(整数)をそれが表す文字に変換したりする方法のメモです。
これらはそれぞれ、 ord と chr という組み込み関数で実現できます。

ドキュメントはこちら。
組み込み関数

これらをアスキーコードへの変換やアスキーコードから文字への変換だと説明しているサイトもあるようですが、
ドキュメントに書かれている通りUnicodeに対応しています。
(もしかしたらpython2系の時代はアスキー文字だけだったのかな?)

組み込み関数なので何もインポートせずに利用可能です。
ただし、ordは文字のみを受けつけ、文字列を渡すとエラーになるので注意してください。

単純なサンプル。


print(ord("a"))  # 97
print(ord("あ"))   # 12354
print(chr(97))   # a
print(chr(12354))   # あ

numpyの線形代数モジュールで連立一次方程式を解く

numpyで作成できる配列によく似たオブジェクトのarrayですが、配列同様にネストして多次元のarrayを作成できます。
(日頃から意識せず普通に使ってますね。)
特に2次元のarrayについては、行列として扱うことが可能です。

そしてnumpyには、numpy.linalgという線形代数のモジュールがあり、
行列式や固有値、逆行列の計算などができます。

ドキュメントはこちら。
Linear algebra (numpy.linalg)

今回は基本的な例として次の連立一次方程式を解いてみましょう。
$$
\begin{eqnarray}
&x_0-4&x_1+2&x_2&=7\\
9&x_0-5&x_1+2&x_2&=-2\\
3&x_0-10&x_1+5&x_2&=3
\end{eqnarray}
$$

これは
$$
A = \left[ \begin{matrix}
1&-4&2\\
9&-5&2\\
3&-10&5
\end{matrix}\right]
$$
$$
b = \left[ \begin{matrix}
7\\
-2\\
-3
\end{matrix}\right]
$$
とおいて、 Aに逆行列が存在すれば、(つまりAの行列式が0でなければ)$A^{-1}b$を計算することで解けます。

それではやってみましょう。
行列式はlinalg.det(A)
逆行列はlinalg.inv(A)で算出できます。
行列の式はnp.dot(A,B)A.dot(B)か、A@Bなどの書き方があります。


# 行列の定義
A = np.array(
        [
            [1,  -4, 2],
            [9,  -5, 2],
            [3, -10, 5]
        ]
    )
b = np.array([7, -2, 3]).T

# Aの行列式の確認 (実装の都合でわずかに誤差が出ます。)
print("dat(A)=", np.linalg.det(A))

# 方程式の解
print("[x_0, x_1, x_2]=", np.linalg.inv(A)@b)

# 以下出力
dat(A)= 0.9999999999999893
[x_0, x_1, x_2]= [ -29. -223. -428.]

これで、連立一次方程式が解けました。
Aの行列式が微妙に1にならないのは実装されていているアルゴリズムの都合のようです。

Sympyで方程式を解く

一つ前の自己回帰過程のサンプルを算出する記事のコードの中で、
こっそりとSympyというのを使っていました。

これはpythonで数式処理を行うためのライブラリです。
簡単な使い方などは公式のチュートリアルのイントロを見ていただくとわかります。
SymPy Tutorial » Introduction

これで記事を終えてしまうとあんまりなので、今回は方程式を解く方法紹介します。
シンプルに
$$x^2-5x+6=0$$
を$x$について解いて $x=2,3$という解を得るコードです。


import sympy
# 記号として使いたい文字をsymbols関数で定義
x = sympy.symbols('x')
# 多項式の定義
f = x**2 - 5*x + 6
print(f)  # 出力 : x**2 - 5*x + 6
# solveに渡すことで、 f = 0 という方程式を xについて解いてくれます。
sympy.solve(f, x)  # 出力 [2, 3]

多項式の展開や因数分解、微積分など、結構色々なことができるライブラリなので便利です。
このブログでも順次各機能を紹介したいと思います。

自己回帰過程のサンプル

今回は移動平均過程の例をいくつか生成して紹介します。
$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点の周期性を持っていること意味するものではありません。