AR過程

直感 #

今日の最高気温は昨日の最高気温と 2〜3 ℃ しか違わず、昨日は一昨日に近い。けれど昨日の気温を「そのまま」引きずるのではなく、平年値に向かって少し戻る。この「昨日の値 × 1 未満の係数 + 揺らぎ」が AR(1) 過程で、係数 \(\phi\) が「昨日をどれだけ引きずるか」を表します。\(\phi\) が 1 だと戻る力がゼロになって株価のようにどこへでも行き、1.04 のようにわずかでも 1 を超えると雪だるま式に膨らんで 500 点のうちに桁が変わります。

式で書くと何を引きずっているか #

AR(p) 過程は、過去 \(p\) 期の自分自身の値に係数 \(\phi_1, \dots, \phi_p\) を掛けて足し、定数 \(c\) とノイズ \(\epsilon_t\) を加えたものです。

AR(1) なら \(y_t = c + \phi_1 y_{t-1} + \epsilon_t\) で、\(|\phi_1| < 1\) のとき平均は \(c / (1 - \phi_1)\) に落ち着きます。下のコードでは \(c = 1\)、\(\phi_1 = 0.1\) なので、系列は \(1 / 0.9 \approx 1.11\) の周りで揺れることになります。MA過程 との違いは、引きずる相手が過去の値 \(y_{t-1}\) か過去のノイズ \(\epsilon_{t-1}\) かという点です。

使うライブラリ #

1
2
3

import statsmodels.api as sm
import numpy as np
import matplotlib.pyplot as plt

AR過程のデータを作成する #

次の関数は上の式を for ループで書き下したものです。phis の長さが次数 \(p\) になり、c が定数項、sigma がノイズの標準偏差を決めます。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

def create_ARdata(phis=[0.1], N=500, init=1, c=1, sigma=0.3):
    """AR過程のデータを作成する"""
    print(f"==AR({len(phis)})過程の長さ{N}のデータを作成==")
    data = np.zeros(N)
    data[0] = init + np.random.normal(0, sigma)

    for t in range(2, N):
        res = c + np.random.normal(0, sigma)
        for j, phi_j in enumerate(phis):
            res += phi_j * data[t - j - 1]
        data[t] = res
    return data

係数が1より小さい場合（定常） #

\(\phi_1 = 0.1\) で 500 点を生成します。平均 \(1/0.9 \approx 1.11\) の周りを離れない、定常な系列になります。

1
2
3
4
5
6

plt.figure(figsize=(12, 6))
phis = [0.1]
ar1_1 = create_ARdata(phis=phis)
plt.plot(ar1_1)
plt.title(f"係数が{phis[0]}である場合", fontsize=15)
plt.show()

==AR(1)過程の長さ500のデータを作成==

係数が１である場合（ランダムウォーク） #

\(\phi_1 = 1\) にすると「昨日の値 + 1 + ノイズ」になり、平均に戻る力が消えます。

1
2
3
4
5
6

plt.figure(figsize=(12, 6))
phis = [1]
ar1_2 = create_ARdata(phis=phis)
plt.plot(ar1_2)
plt.title(f"係数が{phis[0]}である場合", fontsize=15)
plt.show()

==AR(1)過程の長さ500のデータを作成==

係数が１より大きい場合（発散） #

\(\phi_1 = 1.04\) と、1 をわずかに超えただけで系列は指数的に膨らみます。

1
2
3
4
5
6

plt.figure(figsize=(12, 6))
phis = [1.04]
ar1_2 = create_ARdata(phis=phis)
plt.plot(ar1_2)
plt.title(f"係数が{phis[0]}である場合", fontsize=15)
plt.show()

==AR(1)過程の長さ500のデータを作成==

何が起きているか（ステップで見る） #

コードを読む前に、この手法が図の上で何をしているかを追ってみましょう。「再生」を押すか、スライダーを動かすと 1 手順ずつ進みます。

見るポイント：φ が 1 に近いほど前の値を引きずり、波がゆっくりになる。負なら 1 期ごとに符号が反転しやすい。

AR(2) #

係数を 2 個にすると、2 期前の値も今日に影響します。まず \(\phi = (0.1, 0.3)\) の定常な例を 100 点だけ描きます。

1
2
3
4
5
6

plt.figure(figsize=(12, 6))
phis = [0.1, 0.3]
ar2_1 = create_ARdata(phis=phis, N=100)
plt.plot(ar2_1)
plt.title(f"係数が{phis}である場合", fontsize=15)
plt.show()

==AR(2)過程の長さ100のデータを作成==

次に \(\phi_2 = -1.11\) と、定常条件 \(|\phi_2| < 1\) を破る係数にしてみます。

1
2
3
4
5
6

plt.figure(figsize=(12, 6))
phis = [0.1, -1.11]
ar2_1 = create_ARdata(phis=phis)
plt.plot(ar2_1)
plt.title(f"係数が{phis}である場合", fontsize=15)
plt.show()

==AR(2)過程の長さ500のデータを作成==

モデルの推定 #

生成した ar1_1（真の係数 0.1）に AutoReg を当て、係数が復元できるかを確かめます。情報量規準 AIC・HQIC・BIC は次数選びの比較に使う値です。

1
2
3
4
5
6

from statsmodels.tsa.ar_model import AutoReg

res = AutoReg(ar1_1, lags=1).fit()

out = "AIC: {0:0.3f}, HQIC: {1:0.3f}, BIC: {2:0.3f}"
print(out.format(res.aic, res.hqic, res.bic))

AIC: 231.486, HQIC: 236.445, BIC: 244.124

1
2
3

print(res.params)
print(res.sigma2)
res.summary()

[1.03832755 0.07236388]
0.09199676371696269

使いどころと落とし穴 #

AR 過程を実データに当てる前に、まず図と ADF 検定で定常かどうかを確かめます。\(\phi = 1\) の例のように単位根があると AutoReg の係数は信頼できず、差分を取って定常にしてから当てる必要があります。次数 \(p\) は PACF がどのラグで途切れるかを目安にし、迷ったら lags を変えて AIC か BIC が最小になるものを選びます。AR は「過去の値を引きずる」成分しか表せないので、ショックが数期だけ尾を引く構造は MA過程、両方が混ざるなら ARMA を検討してください。