時系列分析で有名な「計量時系列分析」をpythonで実装しながら読み進めていきます。
汚いコードですが自主学習・記録が目的ですのでご了承ください。
様々な人がより分かりやすい記事を書いていますのでそちらをご覧ください。
そして、本格的に勉強したい方・細かい内容を知りたい方はぜひ参考書を手に取って勉強してください。
また、何か問題があればご連絡いただければと思います。

次の本を扱います。

参考書

前回

時系列分析の基礎概念

ARMA過程の性質

MA過程

移動平均(MA)過程は、ホワイトノイズの線形和で表される。1次MA過程(MA(1)過程)は、
$$
y_t=\mu+\varepsilon_t+\theta_1\varepsilon_{t-1},\ \varepsilon_t\sim W.N.(\sigma^2)
$$
で定義され、$y_t$がMA(1)に従うことは$y_t\sim MA(1)$と表記される。
$$
y_{t-1}=\mu+\varepsilon_{t-1}+\theta_1\varepsilon_{t-2}
$$
となるため、$y_t$と$y_{t-1}$が$\varepsilon_{t-1}$という共通項を持つので、$y_t$と$y_{t-1}$の間に相関が生じる。

MA(1)過程の期待値は

\begin{align}
E(y_t)&=E(\mu+\varepsilon_t+\theta_1\varepsilon_{t-1})\\
&=\mu
\end{align}

である。分散は、

\begin{align}
\gamma_0&=Var(y_t)\\
&=Var(\mu+\varepsilon_t+\theta_1\varepsilon_{t-1})\\
&=(1+\theta_1^2)\sigma^2
\end{align}

となる。

import numpy as np
import matplotlib.pyplot as plt

def make_ma1(mu, theta1, sigma, t=100):
    y = [0]
    e = [np.random.normal(scale=sigma)]
    for _ in range(t):
        e.append(np.random.normal(scale=sigma))
        y.append(mu+e[-1]+theta1*e[-2])
    return y

Mu = [0, 2, -2, 0, 2, -2]
Theta1 = [0.8, 0.5, 0.3, -0.3, -0.5, -0.8]
Sigma = [1, 0.5, 2, 1, 0.5, 2]
    
fig, ax = plt.subplots(3, 2, figsize=(10,10))
for i, (mu, theta1, sigma) in enumerate(zip(Mu, Theta1, Sigma)):
    y = make_ma1(mu, theta1, sigma, t=100)
    print('mu='+str(mu)+' theta='+str(theta1)+' sigma='+str(sigma))
    print('   var(theoritical): ',(1+theta1**2)*sigma**2)
    print('   var(observed): ', np.var(y))
    ax[i//2][i%2].plot(y)
    ax[i//2][i%2].set_title('$\mu=$'+str(mu)+' $\\theta=$'+str(theta1)+' $\sigma=$'+str(sigma))

mu=0 theta=0.8 sigma=1
   var(theoritical):  1.6400000000000001
   var(observed):  2.0290072630090026
mu=2 theta=0.5 sigma=0.5
   var(theoritical):  0.3125
   var(observed):  0.31979682722513064
mu=-2 theta=0.3 sigma=2
   var(theoritical):  4.36
   var(observed):  4.640198830929816
mu=0 theta=-0.3 sigma=1
   var(theoritical):  1.09
   var(observed):  0.9864495144307724
mu=2 theta=-0.5 sigma=0.5
   var(theoritical):  0.3125
   var(observed):  0.3265404932555239
mu=-2 theta=-0.8 sigma=2
   var(theoritical):  6.5600000000000005
   var(observed):  5.242312959410386

statsmodelsのARIMA関数でも同様にデータを生成することができる。

from statsmodels.tsa.arima.model import ARIMA

arima = ARIMA(endog=[0],order=(0, 0, 1))

fig, ax = plt.subplots(3, 2, figsize=(10,10))
for i, (mu, theta1, sigma) in enumerate(zip(Mu, Theta1, Sigma)):
    y = arima.simulate(params=(mu,theta1,sigma**2),nsimulations=100)
    print('mu='+str(mu)+' theta='+str(theta1)+' sigma='+str(sigma))
    print('   var(theoritical): ',(1+theta1**2)*sigma**2)
    print('   var(observed): ', np.var(y))
    ax[i//2][i%2].plot(y)
    ax[i//2][i%2].set_title('$\mu=$'+str(mu)+' $\\theta=$'+str(theta1)+' $\sigma=$'+str(sigma))

mu=0 theta=0.8 sigma=1
   var(theoritical):  1.6400000000000001
   var(observed):  1.647123665369126
mu=2 theta=0.5 sigma=0.5
   var(theoritical):  0.3125
   var(observed):  0.3435301990816707
mu=-2 theta=0.3 sigma=2
   var(theoritical):  4.36
   var(observed):  4.54071052662927
mu=0 theta=-0.3 sigma=1
   var(theoritical):  1.09
   var(observed):  0.9764823330598158
mu=2 theta=-0.5 sigma=0.5
   var(theoritical):  0.3125
   var(observed):  0.30295131940888587
mu=-2 theta=-0.8 sigma=2
   var(theoritical):  6.5600000000000005
   var(observed):  6.017607112367808

1次の自己共分散は

\begin{align}
\gamma_1&=Cov(y_t,y_{t-1})\\
&=Cov(\mu+\varepsilon_t+\theta_1\varepsilon_{t-1},\mu+\varepsilon_{t-1}+\theta_1\varepsilon_{t-2})\\
&=Cov(\varepsilon_t,\varepsilon_{t-1})+Cov(\varepsilon_t,\theta_1\varepsilon_{t-2})+Cov(\theta_1\varepsilon_{t-1},\varepsilon_{t-1})+Cov(\theta_1\varepsilon_{t-1},\theta_1\varepsilon_{t-2})\\
&=\theta_1Cov(\varepsilon_{t-1},\varepsilon_{t-1})\\
&=\theta_1\sigma^2
\end{align}

と求めることができる。これより、MA(1)過程の1次自己相関が

\rho_1=\frac{\gamma_1}{\gamma_0}=\frac{\theta_1}{1+\theta_1^2}

で与えられることがわかる。2次以降の自己共分散は、$k\geq2$として、

\begin{align}
\gamma_1&=Cov(y_t,y_{t-k})\\
&=Cov(\mu+\varepsilon_t+\theta_1\varepsilon_{t-1},\mu+\varepsilon_{t-k}+\theta_1\varepsilon_{t-k-1})\\
&=0
\end{align}

となることがわかる。

MA(1)過程の自己相関係数の計算を実装する。

def calc_autocor(y):
    """自己相関係数の計算
    """
    y_bar = np.mean(y)
    gamma_bar = [np.sum((y-y_bar)**2)]
    for i in range(1,T):
        gamma_bar.append(np.sum((y[:-i]-y_bar)*(y[i:]-y_bar)))
    rho_bar = gamma_bar / gamma_bar[0]
    return rho_bar

T = 100
mu=0
theta1 = 0.5
sigma=1

fig, ax = plt.subplots(2,2,figsize=(10,6))

Theta1 = [0.8, -0.8]
for i, theta1 in enumerate(Theta1):

    y = make_ma1(mu=mu, theta1=theta1, sigma=sigma, t=T)
    rho_bar = calc_autocor(y)

    print('ρ1: ', theta1/(1+theta1**2))

    ax[i][0].plot(y)
    ax[i][0].set_title('$\\theta_1=$'+str(theta1))
    
    ax[i][1].bar(np.arange(0,10), rho_bar[:10])
    ax[i][1].hlines(0, -1, 10, linestyle='--', color='black')
    ax[i][1].hlines([-1.96/np.sqrt(T),1.96/np.sqrt(T)], -1, 10, linestyle='--', color='red')
    ax[i][1].set_xticks(np.arange(0,10));

ρ1:  0.4878048780487805
ρ1:  -0.4878048780487805

一般的に$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)

で定義され、MA($q$)過程と表記される。

定理(MA(q)過程の性質)
MA($q$)過程は以下の性質を持つ
1.$E(y_t)=\mu$
2.$\gamma_0=Var(y_t)=(1+\theta_1^2+\theta_2^2+\cdots+\theta_q^2)\sigma^2$
3.

\gamma_k=\left\{\begin{array}{ll}
(\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{array}\right.

4.MA過程は常に定常である
5.

\rho_k=\left\{\begin{array}{ll}
\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{array}\right.

MA($q$)過程の自己相関係数の計算を実装する。

def calc_rho(theta, k):
    if k > len(theta):
        rho = 0.
    elif k == 0:
        rho = 1.
    else:
        rho = np.sum(theta[k-1:]*np.concatenate([np.ones(1), theta[:len(theta)-k]]))/(np.sum(theta**2)+1)
    return rho


def make_ma2(mu, Theta, sigma, t=100):
    y = [0]
    e = [np.random.normal(scale=sigma) for _ in range(len(Theta))]
    for _ in range(t):
        e.append(np.random.normal(scale=sigma))
        y.append(mu+e[-1]+Theta@np.array(e)[-1-len(Theta):-1][::-1])
    return y

T = 200
mu=0
Theta = np.array([0.8, 0.3])
sigma=1

Thetas = np.array([[0.8,0.3],
                   [0.3,0.8],
                   [0.8,-0.3],
                   [-0.8,-0.3]])

fig, ax = plt.subplots(4,2,figsize=(12,12))
for i, Theta in enumerate(Thetas):
    y = make_ma2(mu, Theta, sigma, T)

    rho_bar = calc_autocor(y)
    rho=np.array([calc_rho(Theta, k) for k in np.arange(20)])

    ax[i][0].plot(y)
    ax[i][0].set_title('$\\theta_1=$'+str(Theta[0])+', $\\theta_2=$'+str(Theta[1]))
    
    ax[i][1].bar(np.arange(0,10), rho[:10], alpha=.4, color='red')
    ax[i][1].bar(np.arange(0,10), rho_bar[:10], alpha=.4)
    
    ax[i][1].hlines(0, -1, 10, linestyle='--', color='black')
    ax[i][1].hlines([-1.96/np.sqrt(T),1.96/np.sqrt(T)], -1, 10, linestyle='--', color='red')
    ax[i][1].set_xticks(np.arange(0,10));

statsmodelsを使うことでより簡単に実装できる。

import statsmodels.api as sm

arima = ARIMA(endog=[0],order=(0, 0, 2))

fig, ax = plt.subplots(4,2,figsize=(12,12))
for i, (mu, Theta, sigma) in enumerate(zip(Mu, Thetas, Sigma)):
    y = arima.simulate(params=[mu]+list(Theta)+[sigma**2],nsimulations=T)

    rho=np.array([calc_rho(Theta, k) for k in np.arange(20)])

    ax[i][0].plot(y)
    ax[i][0].set_title('$\\theta_1=$'+str(Theta[0])+', $\\theta_2=$'+str(Theta[1]))
    
    ax[i][1].bar(np.arange(0,10), rho[:10], alpha=.4, color='red')
    sm.graphics.tsa.plot_acf(y, lags=np.arange(1,10), ax=ax[i][1], auto_ylims=True);

AR過程

自己回帰(AR)過程は、過程が自身の過去に回帰された形で表現される過程である。
1次AR過程(AR(1)過程)は、
$$
y_t=c+\phi_1y_{t-1}+\varepsilon_t,\ \varepsilon_t\sim W.N.(\sigma^2)
$$
で定義され、$y_t$がAR(1)に従うことは$y_t\sim AR(1)$と表記される。
AR(1)過程の場合は、$|\phi_1|<1$のときに過程は定常となる。一方$|\phi_1|>1$のモデルは、過程が指数的に上昇しており、このような過程は爆発的と呼ばれる。$|\phi_1|=1$のモデルは単位根過程と呼ばれる。
AR(1)過程の期待値は
$$
\mu=\frac{c}{1-\phi_1}
$$
で得られる。また、分散は
$$
\gamma_0=\frac{\sigma^2}{1-\phi_1^2}
$$
で与えられる。
$k$次自己共分散は、
$$
\gamma_k=\phi_1\gamma_{k-1}
$$
が得られる。両辺を$\gamma_0$で割ると、
$$
\rho_k=\phi_1\rho_{k-1}
$$
というユール・ウォーカー方程式が得られる。
一般的に自己相関を求めることも容易であり
$$
\rho_k=\phi_1^2\rho_{k-2}=\phi_1^3\rho_{k-3}=\cdots=\phi_1^k\rho_0=\phi_1^k
$$
となる。
この結果と$|\phi_k|<1$より、AR(1)過程の自己相関の絶対値は指数的に減衰していくことがわかる。

def make_ar1(c, phi1, sigma=1, t=200):
    y = [c]
    for _ in range(t):
        e = np.random.normal(scale=sigma)
        y.append(c+phi1*y[-1]+e)
    return y

C = [2, -2, 0, -2, 2, 2]
Phi1 = [0.8, 0.3, -0.3, -0.8, 1, 1.1]
Sigma = [1, 0.5, 2, 1, 1, 1]

fig, ax = plt.subplots(3, 2, figsize=(10,10))
for i, (c, phi1, sigma) in enumerate(zip(C, Phi1, Sigma)):
    y = make_ar1(c, phi1, sigma, t=200)
    if np.abs(phi1) < 1:
        print('c='+str(mu)+' phi='+str(theta1)+' sigma='+str(sigma))
        print('   mu(theoritical): ',c/(1-phi1))
        print('   mu(observed): ', np.mean(y[50:]))
        print('   var(theoritical): ',sigma**2/(1-phi1**2))
        print('   var(observed): ', np.var(y[50:]))
    ax[i//2][i%2].plot(y)
    ax[i//2][i%2].set_title('$c=$'+str(mu)+' $\\phi=$'+str(theta1)+' $\sigma=$'+str(sigma))

c=0 phi=-0.8 sigma=1
   mu(theoritical):  10.000000000000002
   mu(observed):  10.423554481938877
   var(theoritical):  2.7777777777777786
   var(observed):  2.0199785950212443
c=0 phi=-0.8 sigma=0.5
   mu(theoritical):  -2.857142857142857
   mu(observed):  -2.8970011449602753
   var(theoritical):  0.2747252747252747
   var(observed):  0.27867309435515886
c=0 phi=-0.8 sigma=2
   mu(theoritical):  0.0
   mu(observed):  -0.04118207753762859
   var(theoritical):  4.395604395604395
   var(observed):  4.4063065350110735
c=0 phi=-0.8 sigma=1
   mu(theoritical):  -1.1111111111111112
   mu(observed):  -1.1117681465144518
   var(theoritical):  2.7777777777777786
   var(observed):  3.083649828299484

AR(1)過程について自己相関係数を計算する。

T = 200
c=2
sigma=1

fig, ax = plt.subplots(2,2,figsize=(10,6))

Phi1 = [0.8, -0.8]
for i, phi1 in enumerate(Phi1):

    y = make_ar1(c, phi1, sigma, t=T+10)[10:]
    
    rho_bar = calc_autocor(y)
    rho = phi1**np.arange(20)

    ax[i][0].plot(y)
    ax[i][0].set_title('$\\phi_1=$'+str(phi1))
    
    ax[i][1].bar(np.arange(0,20), rho[:20], alpha=.4, color='red')
    ax[i][1].bar(np.arange(0,20), rho_bar[:20], alpha=.4)
    
    ax[i][1].hlines(0, -1, 20, linestyle='--', color='black')
    ax[i][1].hlines([-1.96/np.sqrt(T),1.96/np.sqrt(T)], -1, 20, linestyle='--', color='red')
    ax[i][1].set_xticks(np.arange(0,20));

ここでもstatsmodelを使用して簡略化したものを示す。

arima = ARIMA(endog=[0],order=(1, 0, 0))

fig, ax = plt.subplots(2,2,figsize=(10,6))
for i, phi1 in enumerate(Phi1):
    y = arima.simulate(params=[c,phi1,sigma**2],nsimulations=T)

    rho = phi1**np.arange(20)

    ax[i][0].plot(y)
    ax[i][0].set_title('$\\phi_1=$'+str(phi1))
    
    ax[i][1].bar(np.arange(0,10), rho[:10], alpha=.4, color='red')
    sm.graphics.tsa.plot_acf(y, lags=np.arange(1,10), ax=ax[i][1], auto_ylims=True);

$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)
$$
で定義される。

定理(定常AR(p)過程の性質)
定常AR($p$)過程は以下の性質を持つ

$\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}$

自己共分散と自己相関は$y_t$が従うAR過程と同一の係数を持つ以下の$p$次差分方程式に従う
$$
\begin{align}
\gamma_k &= \phi_1\gamma_{k-1}+\phi_2\gamma_{k-2}+\cdots+\phi_p\gamma_{k-p}\ k\geq 1\
\rho_k &= \phi_1\rho_{k-1}+\phi_2\rho_{k-2}+\cdots+\phi_p\rho_{k-p}\ k\geq 1
\end{align}
$$
第2式はユール・ウォーカー方程式と呼ばれる。

AR過程の自己相関は指数的に減衰する

AR(2)過程のユール・ウォーカー方程式は

$$
\rho_k=\phi_1\rho_{k-1}+\phi_2\rho_{k-2}
$$
となる。
$k=1$を代入して$\rho_0=1$と$\rho_1=\rho_{-1}$に注意すると、
$$
\rho_1=\phi_0+\phi_2\rho_{-1}=\phi_1+\phi_2\rho_1
$$
が得られる。したがって、
$$
\rho_1=\frac{\phi_1}{1-\phi_2}
$$
次に、$k=2$を代入すると、
$$
\rho_2=\frac{\phi_1^2+\phi_2-\phi_2^2}{1-\phi_2}
$$
となる。

AR($p$)についての自己相関係数を計算する。

def make_ar2(c, Phi, sigma=1, t=200):
    y = [c, c]
    for _ in range(t):
        e = np.random.normal(scale=sigma)
        y.append(c+Phi[0]*y[-1]+Phi[1]*y[-2]+e)
    return y

def calc_rho2(Phi, t=200):
    rho = [1,Phi[0]/(1-Phi[1])]
    for _ in range(t):
        rho.append(Phi[0]*rho[-1]+Phi[1]*rho[-2])
    return rho

T = 200
c=2
sigma=1

fig, ax = plt.subplots(4,2,figsize=(12,12))

Phis = [[0.5,0.35], [0.1,0.5],[0.5,-0.8], [0.9,-0.8]]
for i, Phi in enumerate(Phis):

    y = make_ar2(c, Phi, sigma, t=T+10)[10:]
    
    rho_bar = calc_autocor(y)
    rho = calc_rho2(Phi, 20)

    ax[i][0].plot(y)
    ax[i][0].set_title('$\\phi_1=$'+str(Phi[0])+', $\\phi_2=$'+str(Phi[1]))
    
    ax[i][1].bar(np.arange(0,20), rho[:20], alpha=.4, color='red')
    ax[i][1].bar(np.arange(0,20), rho_bar[:20], alpha=.4)
    
    ax[i][1].hlines(0, -1, 20, linestyle='--', color='black')
    ax[i][1].hlines([-1.96/np.sqrt(T),1.96/np.sqrt(T)], -1, 20, linestyle='--', color='red')
    ax[i][1].set_xticks(np.arange(0,20));

次に、SARIMAX関数を使用して出力したデータに対して自己相関係数を求める。

Phis = [[0.5,0.35], [0.1,0.5],[0.5,-0.8], [0.9,-0.8]]
arima = sm.tsa.SARIMAX([0], order=(2,0,0))

fig, ax = plt.subplots(4,2,figsize=(12,12))
for i, Phi in enumerate(Phis):

    y = arima.simulate(params=list(Phi)+[c,sigma**2],nsimulations=T)
    
    rho_bar = calc_autocor(y)
    rho = calc_rho2(Phi, 10)

    ax[i][0].plot(y)
    ax[i][0].set_title('$\\phi_1=$'+str(Phi[0])+', $\\phi_2=$'+str(Phi[1]))
    
    ax[i][1].bar(np.arange(0,10), rho[:10], alpha=.4, color='red')
    sm.graphics.tsa.plot_acf(y, lags=np.arange(1,10), ax=ax[i][1], auto_ylims=True);

ARMA過程

自己回帰移動平均(ARMA)過程は自己回帰項と移動平均項を両方含んだ項である。
$(p,q)$次ARMA過程(ARMA($p,q$)過程)は
$$
y_t=c+\phi_1y_{t-1}\cdots+\phi_py_{t-p}+\varepsilon_t+\theta_1\varepsilon_{t-1}\cdots+\theta_q\varepsilon_{t-q}\ \varepsilon_t\sim W.N.(\sigma^2)
$$
で定義される。

定義(ARMA(p,q)過程の性質)
定常ARMA($p,q$)過程は以下の性質を持つ。

$\mu=E(y_t)=\frac{c}{1-\phi_1-\phi_2-\cdots-\phi_p}$

$q+1$次以降の自己共分散と自己相関は$y_t$が従うARMA過程のAR部分と同一の係数をもつ$q$次差分方程式(ユール・ウォーカー方程式)に従う。
$$
\begin{align}
\gamma_k &= \phi_1\gamma_{k-1}+\phi_2\gamma_{k-2}+\cdots+\phi_p\gamma_{k-p}\ k\geq q+1\
\rho_k &= \phi_1\rho_{k-1}+\phi_2\rho_{k-2}+\cdots+\phi_p\rho_{k-p}\ k\geq q+1
\end{align}
$$

ARMA過程の自己相関は指数的に減衰する

$q$次までの自己共分散と自己相関はMA項の影響があるため、一般的に表現するのは難しい。

def make_arma(c, Phi, Theta, sigma=1, t=200):
    y = [c]*(len(Phi))
    e = [np.random.normal(scale=sigma) for _ in range(len(Theta))]
    for _ in range(t+10):
        et = np.random.normal(scale=sigma)
        y.append(c+Phi@np.array(y)[-len(Phi):][::-1]+Theta@np.array(e)[-len(Theta):][::-1]+et)
        e.append(et)
    return y[len(Phi)+10:]

c = 1
Phi=np.array([0.3,0.5])
Theta=np.array([0.3,0.5])

y = make_arma(c=c, Phi=Phi, Theta=Theta, sigma=.2, t=300)

print('mu(theoritical): ',c/(1-np.sum(Phi)))
print('mu(observed): ', np.mean(y[50:]))

plt.plot(y)

mu(theoritical):  5.000000000000001
mu(observed):  5.045790843274028

ここでは、arma_generate_sample関数を使用して同様の仮定を出力する。
ARIMA関数でもSARIMAX関数でも同様に出力できる。

from statsmodels.tsa.arima_process import arma_generate_sample

y = arma_generate_sample(ar=np.r_[1,-Phi],
                         ma=np.r_[1,Theta],
                         scale=.2,
                         nsample=300)+c/(1-np.sum(Phi))

print('mu(theoritical): ',c/(1-np.sum(Phi)))
print('mu(observed): ', np.mean(y[50:]))

plt.plot(y)

mu(theoritical):  5.000000000000001
mu(observed):  4.944490318078968

ARMA過程の定常性と反転可能性

AR過程の定常性

AR($p$)過程は、
$$
1-\phi_1z-\cdots-\phi_pz^p=0
$$
の方程式のすべての解の絶対値が1より大きいとき、AR過程は定常となる。
この方程式はAR特性方程式と呼ばれ、特性方程式の左辺の多項式はAR多項式と呼ばれる。

AR(1)過程
$$
y_t=c+\phi_1y_{t-1}+\varepsilon_t\ \varepsilon_t\sim W.N.(\sigma^2)
$$
の定常条件は、特性方程式$1-\phi_1z=0$の解が$z=\phi_1^{-1}$で与えられることから、$|\phi_1|<1$である。

AR(2)過程
$$
y_t=c+\phi_1y_{t-1}+\phi_2y_{t-2}+\varepsilon_t\ \varepsilon_t\sim W.N.(\sigma^2)
$$
の定常条件は、

\begin{align}
\phi_2&>-1\\
\phi_2&<1+\phi_1\\
\phi_2&<1-\phi_1
\end{align}

で与えられることが知られている。

定常AR過程はMA($\infty$)過程で書き直すことができるという性質もある。
AR(1)過程
$$
y_t=c+\phi_1y_{t-1}+\varepsilon_t\ \varepsilon_t\sim W.N.(\sigma^2)
$$
を繰り返し代入していくことで、

\begin{align}
y_t&=\phi_1y_{t-1}+\varepsilon_t\\
&=\phi_1(\phi_1y_{t-2}+\varepsilon_{t-1})+\varepsilon_{t-1}\\
&\vdots\\
&=\phi_1^my_{t-m}+\sum_{k=0}^{m-1}\phi_1^k\varepsilon_{t-k}
\end{align}

となることがわかる。
ここで、$|\phi_1|<1$ならば、$m→\infty$のとき$\phi_1^my_{t-m}→0$が成立するので、
$$
y_t=\sum_{k=0}^\infty\phi_1^k\varepsilon_{t-k}
$$
と書き直すことができる。

これは次のように固有方程式を解くことで確認できる。

Phis = [[0.5,0.35], [1.01, .01]]

fig, ax = plt.subplots(1,2,figsize=(12,4))
for i, Phi in enumerate(Phis):
    A = np.concatenate([np.array(Phi).reshape(1,-1), np.eye(len(Phi))], axis=0)[:-1]
    eig_values, eig_vectors = np.linalg.eig(A)

    print('φ1='+str(Phi[0])+', φ2='+str(Phi[1]))
    print('eigen values: ', eig_values)
    print('judge: ', np.all(eig_values<1))
    
    y = make_ar2(c, Phi, sigma, t=T+10)[10:]
    ax[i].plot(y)
    ax[i].set_title('$\\phi_1=$'+str(Phi[0])+', $\\phi_2=$'+str(Phi[1]))

φ1=0.5, φ2=0.35
eigen values:  [ 0.89226163 -0.39226163]
judge:  True
φ1=1.01, φ2=0.01
eigen values:  [ 1.01980579 -0.00980579]
judge:  False

MA過程の反転可能性

MA過程は常に定常であるが、同一の期待値と自己相関構造を持つ異なるMA過程が存在するという問題がある。
どのMA過程を用いるかの判断基準の1つとして、MA過程の反転可能性がある。

定義(MA過程の反転可能性)
MA過程がAR($\infty$)に書き直せるとき、MA過程は反転可能といわれる

MA過程が反転可能のとき、撹乱項$\varepsilon_t$は過去の$y_t$の関数として表現でき、さらに過去の$y$を用いて$y_t$を予測したときの予測誤差と解釈できる。
このため、反転可能表現に伴う$\varepsilon_t$を$y_t$の本源的な撹乱項と呼ぶこともある。
同一の期待値と自己相関構造をもつ複数のMA過程のうち、反転可能な過程を選択する方が望ましい。
MA($q$)過程の反転可能条件は、AR過程の定常条件と同様のものであり、

$$
1+\theta_1z+\theta_2z^2+\cdots+\theta_pz^p=0
$$
というMA特性方程式を用いて調べることができる。
すべての解の絶対値が1より大きいとき、MA過程は反転可能となる。
例えば、MA(1)過程
$$
y_t=\varepsilon_t+\theta_1\varepsilon_{t-1},\ \varepsilon_t\sim W.N.(\sigma^2)
$$
の反転可能条件は$|\theta_1|<1$となる。

$\varepsilon_t$をAR過程で書き直すと、

\begin{align}
\varepsilon_t&=-\theta\varepsilon_{t-1}+y_t\\
&=(-\theta)^m\varepsilon_{t-m}+\sum_{k=0}^{m-1}(-\theta)^ky_{t-k}\\
&→\sum_{k=0}^{m-1}(-\theta)^ky_{t-k}
\end{align}

したがって、MA(1)過程は反転可能なとき、
$$
y_t=-\sum_{k=0}^{m-1}(-\theta)^ky_{t-k}+\varepsilon_t
$$
というAR($\infty$)過程で書き直すことができる。

ARMA過程の定常・反転可能性

MA過程は常に定常であるので、AR過程の部分が定常であればARMA過程は定常となる。
AR特性方程式のすべての解の絶対値が1より大きければ、ARMA過程は定常となる。
反転可能となるには、MA過程の部分が反転可能となればよい。
MA過程の部分のMA特性方程式のすべての解の絶対値が1より大きければ、ARMA過程は反転可能となる。

ARMAモデルの推定

最小2乗法

最も基本的な方法は最小2乗法(OLS) である。
以下では、
$$
y_t=c+\phi y_{t-1}+\varepsilon_t\ \varepsilon_t\sim iid(0,\sigma^2)
$$

というAR(1)モデルを用いてOLSの説明を行う。
OLSでは推定するモデルを回帰モデルと呼び、左辺の変数$y_t$のことを被説明変数と呼ぶ。$\varepsilon_t$のことを誤差項と呼び、誤差項以外の右辺の変数のことを説明変数と呼ぶ。
回帰モデルにおける$c$と$\phi$の任意の推定量を$\tilde{c}$と$\tilde{\phi}$とすると、OLSでは、モデルが説明できない部分
$$
e_t=y_t-\tilde{c}-\tilde{\phi}y_{t-1}
$$
を残差と呼び、残差$e_t$を平均的に0に近くすることを考える。OLSでは残差平方和(SSR)
$$
SSR=\sum_{t=1}^Te_t^2=\sum_{t=1}^T(y_t-\tilde{c}-\tilde{\phi}y_{t-1})^2
$$
が最小になるように$\tilde{c}$と$\tilde{\phi}$を決める。
正規方程式を解くことでOLS推定量が、

\begin{align}
\hat{c}&=\bar{y}_t-\hat{\phi}\bar{y}_{t-1}\\
\hat{\phi}&=\frac{\sum_{t=1}^T(y_t-\bar{y}_t)(y_{t-1}-\bar{y}_{t-1})}{\sum_{t=1}^T(y_{t-1}-\bar{y}_{t-1})^2}
\end{align}

で与えられることがわかる。ここで、$\bar{y}t=\frac{1}{T}\sum{t=1}^Ty_t$である。
OLS推定量の計算には初期値$y_0$が必要となる。一般的に、AR($p$)モデルの推定には$p$個の初期値が必要となる。
また、撹乱項の分散$\sigma^2$のOLS推定量$s^2$は以下のOLS残差に基づいた不変標本分散で与えられる。
$$
s^2=\frac{1}{T-2}\sum_{t=1}^T(y_t-\hat{c}-\hat{\phi}y_{t-1})^2
$$

定理(OLS推定量の性質)
ARモデルにおけるOLS推定量は以下の性質を持つ

OLS推定量は一致推定量である

OLS推定量を基準化したものは漸近的に正規分布に従う

$\varepsilon_t\sim iid(0,\sigma^2)$のとき、OLS推定量は一致推定量の中で最小の分散共分散行列を持つ

OLSは不変推定量であり、すべての線形不変推定量の中で、最小の分散共分散行列をもつ最良線形不変推定量(BLUE) となることが知られている。

AR(1)過程についての推定結果を示す。

Phi_hat = []
C_hat = []
S2 = []

for _ in range(200):
    y = make_ar1(c=4, phi1=0.2, sigma=1.5, t=200)
    yt = y[1:]
    yt_1 = y[:-1]

    yt_bar = np.mean(yt)
    yt_1_bar = np.mean(yt_1)

    phi_hat = np.sum((yt-yt_bar)*(yt_1-yt_1_bar)) / np.sum((yt_1-yt_1_bar)**2)
    c_hat = yt_bar - phi_hat*yt_1_bar
    s2 = 1/(len(yt)-2)*np.sum((yt-c_hat-phi_hat*yt_1_bar)**2)
    
    Phi_hat.append(phi_hat)
    C_hat.append(c_hat)
    S2.append(s2)

fig, ax = plt.subplots(1,3,figsize=(12,3))

ax[0].hist(Phi_hat, bins=15);
ax[1].hist(C_hat, bins=15);
ax[2].hist(S2, bins=15);

AR(2)についても推定を行う。
sklearnのLinearRegressioを利用した。

from sklearn.linear_model import LinearRegression

Phi1_hat = []
Phi2_hat = []
C_hat = []
S2 = []

for _ in range(200):
    y = make_ar2(c=4, Phi=np.array([0.1,-0.8]), sigma=1.5, t=200)
    yt = np.array(y[2:])
    yt_1 = np.array([y[i:i+2] for i in range(len(y)-2)])
    
    model = LinearRegression().fit(yt_1,yt)
    phi1_hat = model.coef_[0]
    phi2_hat = model.coef_[1]
    
    Phi1_hat.append(model.coef_[0])
    Phi2_hat.append(model.coef_[1])
    C_hat.append(model.intercept_)
    S2.append(np.sum((yt-model.predict(yt_1))**2)/(len(yt)-2))

fig, ax = plt.subplots(1,4,figsize=(14,3))
ax[0].hist(Phi1_hat, bins=15);
ax[1].hist(Phi2_hat, bins=15);
ax[2].hist(C_hat, bins=15);
ax[3].hist(S2, bins=15);

次のようにstatsmodelsの関数を利用できる。

from statsmodels.tsa import ar_model

model = ar_model.AutoReg(y, lags=2).fit()
print(model.summary())

                            AutoReg Model Results                             
==============================================================================
Dep. Variable:                      y   No. Observations:                  202
Model:                     AutoReg(2)   Log Likelihood                -364.608
Method:               Conditional MLE   S.D. of innovations              1.498
Date:                Wed, 11 Jan 2023   AIC                            737.216
Time:                        09:38:05   BIC                            750.410
Sample:                             2   HQIC                           742.555
                                  202                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const          3.7588      0.184     20.441      0.000       3.398       4.119
y.L1           0.1359      0.047      2.878      0.004       0.043       0.229
y.L2          -0.7417      0.047    -15.723      0.000      -0.834      -0.649
                                    Roots                                    
=============================================================================
                  Real          Imaginary           Modulus         Frequency
-----------------------------------------------------------------------------
AR.1            0.0916           -1.1575j            1.1611           -0.2374
AR.2            0.0916           +1.1575j            1.1611            0.2374
-----------------------------------------------------------------------------

ARMA過程については例えばstatsmodelsのSARIMAXを利用する。

import statsmodels.api as sm

T = 300
c = 3
sigma = .2
Phi = np.array([0.4,0.5])
Theta = np.array([0.3,-0.5,0.2,0.4])

y = make_arma(c, Phi, Theta, sigma, T)

model = sm.tsa.statespace.SARIMAX(y, trend='c', order=(2,0,4)).fit()
print(model.summary())

                               SARIMAX Results                                
==============================================================================
Dep. Variable:                      y   No. Observations:                  300
Model:               SARIMAX(2, 0, 4)   Log Likelihood                  -2.339
Date:                Wed, 11 Jan 2023   AIC                             20.677
Time:                        09:38:06   BIC                             50.308
Sample:                             0   HQIC                            32.535
                                - 300                                         
Covariance Type:                  opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
intercept      0.1327      0.090      1.471      0.141      -0.044       0.310
ar.L1          0.7969      0.118      6.743      0.000       0.565       1.028
ar.L2          0.1978      0.118      1.682      0.093      -0.033       0.428
ma.L1          0.2325      0.113      2.051      0.040       0.010       0.455
ma.L2         -0.3062      0.041     -7.431      0.000      -0.387      -0.225
ma.L3          0.3826      0.055      6.971      0.000       0.275       0.490
ma.L4          0.5232      0.064      8.185      0.000       0.398       0.649
sigma2         0.0583      0.005     12.439      0.000       0.049       0.068
===================================================================================
Ljung-Box (L1) (Q):                   1.50   Jarque-Bera (JB):                 4.72
Prob(Q):                              0.22   Prob(JB):                         0.09
Heteroskedasticity (H):               0.82   Skew:                             0.25
Prob(H) (two-sided):                  0.31   Kurtosis:                         3.34
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).


C:\Users\00220900095\.virtualenvs\01.syouseiro-HheMhbJR\lib\site-packages\statsmodels\base\model.py:604: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
  warnings.warn("Maximum Likelihood optimization failed to "

最尤法

最尤法をAR(1)モデルに応用することを考える。$\varepsilon_t\sim iidN(0,\sigma^2)$と仮定して議論する。
$\boldsymbol{\theta}$を推定したいパラメータとすると、AR(1)モデルの場合は、$\boldsymbol{\theta}=(c,\phi,\sigma^2)'$となる。
確率密度をパラメータ$\boldsymbol{\theta}$の関数として見たものを尤度と呼び、この尤度を最大化するものが最尤推定量となる。
ベイズの法則より、$y_1$と$y_2$の同時密度に関して、

$$
\begin{align}
f_{Y_2,Y_1|Y_0}(y_2,y_1|y_0;\boldsymbol{\theta})&=f_{Y_2|Y_1,Y_0}(y_2|y_1,y_0;\boldsymbol{\theta})・f_{Y_1|Y_0}(y_1|y_0;\boldsymbol{\theta})\
&=f_{Y_2|Y_1}(y_2|y_1;\boldsymbol{\theta})・f_{Y_1|Y_0}(y_1|y_0;\boldsymbol{\theta})
\end{align}
$$
が成立する。最後の等号はAR(1)過程の場合、直近1期の値にしか条件付分布が依存しないことによる。
同様に考えると、
$$
f_{Y_T,Y_{T-1},\cdots,Y_1|Y_0}(y_T,y_{T-1},\cdots,y_1|y_0;\boldsymbol{\theta})=\prod_{t=1}^Tf_{Y_t|Y_{t-1}}(y_t|y_{t-1};\boldsymbol{\theta})
$$
が成立することがわかる。したがって、対数尤度は
$$
\mathcal{L}(\boldsymbol{\theta})=\sum_{t=1}^T\log f_{Y_t|Y_{t-1}}(y_t|y_{t-1};\boldsymbol{\theta})
$$
と書ける。ここで、$\varepsilon_t\sim iidN(0,\sigma^2)$のとき、$y_t|y_{t-1}\sim N(c+\phi y_{t-1},\sigma^2)$であるので、
$$
f_{Y_t|Y_{t-1}}(y_t|y_{t-1};\boldsymbol{\theta})=\frac{1}{\sqrt{2\pi\sigma^2}}\exp\bigg[\frac{-(y_t-c-\phi y_{t-1})^2}{2\sigma^2} \bigg]
$$
となる。これより対数尤度が、
$$
\mathcal{L}(\boldsymbol{\theta})=-\frac{T}{2}\log (2\pi)-\frac{T}{2}\log(\sigma^2)-\sum_{t=1}^T\bigg[\frac{-(y_t-c-\phi y_{t-1})^2}{2\sigma^2} \bigg]
$$
これを最大にするような$c$と$\phi$は$(y_t-c-\phi y_{t-1})^2$を最小にするような$c$と$\phi$で与えられるが、これはOLS推定量そのものである。

異常の議論を一般化すると、$t$期までに利用可能な情報の集合、すなわち時点$t-1$の情報集合を$\Sigma_{t-1}$で表すとすると、対数尤度は一般的に
$$
\mathcal{L}(\boldsymbol{\theta})=\sum_{t=1}^T\log f_{Y_t|\Omega_{t-1}}(y_t|\Omega_{t-1};\boldsymbol{\theta})
$$

という形で書くことができるので、この式を最大化するような$\boldsymbol{\theta}$が最尤推定量となる。

ARMAモデルの選択

モデル候補の選択

モデル候補の選択の際に重要なツールとなるのが、標本自己相関係数である。
MA($q$)過程の自己相関$q+1$次以降0になる性質があるが、AR過程とARMA過程の識別は難しい。
ここで、k次偏自己相関とは、単純に$y_t$と$y_{t-k}$から$y_{t-1},\cdots,y_{t-k+1}$の影響を取り除いたものの間の相関を考えたものである。
偏自己相関を$k$の関数として見たものは偏自己相関関数と呼ばれる。
AR($p$)モデルの場合、$y_t$は$y_{t-1},\cdots,y_{t-p}$の線形関数でかけるので、$p+1$次以降の偏自己相関は0となる。
MA過程はAR($\infty$)過程で書き直すことができるので、$y_t$と$y_{t-k}$が相関を持つような$k$は限りなく存在することになる。

モデル	自己相関	偏自己相関
AR(p)モデル	減衰していく	p+1次以降0
MA(q)モデル	p+1次以降0	減衰していく
ARMA(p,q)モデル	減衰していく	減衰していく

注意点として、標本自己相関と標本偏自己相関は、ともに推定誤差を含むことがある。
自己相関や偏自己相関が0であるかどうかの判断は慎重に行わなければならない。
また、減衰していることと切断されていることの差も微妙であり、その判断は主観的にならざるを得ないところがある。
さらに、そもそも標本自己相関と標本偏自己相関だけでは、ARMA過程の次数を決めることはできない。

各過程の自己相関係数、偏自己相関係数を確認する。

def calc_pautocor(y):
    """偏自己相関係数
    """
    rho = []
    for k in range(1,20):
        _y = np.array(y[k:])
        _x = np.array([y[i:i+k] for i in range(len(y)-k)])

        model = LinearRegression().fit(_x,_y)
        rho.append(model.coef_[0])
    return rho

T = 300
c = 1
mu = 1
sigma = 1
Phi = np.array([0.4,0.2])
Theta = np.array([0.7,0.5,0.2,0.8])

y1 = make_ma2(mu, Theta, sigma, T)
y2 = make_ar2(c, Phi, sigma, T)
y3 = make_arma(c, Phi, Theta, sigma, T)

fig, ax = plt.subplots(3,3,figsize=(12,8))
for i, y in enumerate([y1,y2,y3]):
    rho_bar = calc_autocor(y)
    prho = calc_pautocor(y)

    ax[i][0].plot(y)
    ax[i][1].bar(np.arange(0,20), rho_bar[:20], alpha=.4, color='blue')
    ax[i][2].bar(np.arange(1,20), prho[:20], alpha=.4, color='red')
    
    ax[i][1].hlines(0, -1, 20, linestyle='--', color='black')
    ax[i][1].hlines([-1.96/np.sqrt(T),1.96/np.sqrt(T)], -1, 20, linestyle='--', color='red')
    ax[i][2].hlines(0, -1, 20, linestyle='--', color='black')
    ax[i][2].hlines([-1.96/np.sqrt(T),1.96/np.sqrt(T)], -1, 20, linestyle='--', color='red')

statsmodelの関数を利用してより簡単に確認できる。

fig, ax = plt.subplots(3,3,figsize=(12,8))
for i, y in enumerate([y1,y2,y3]):
    ax[i][0].plot(y)
    sm.graphics.tsa.plot_acf(y, lags=np.arange(1,20), ax=ax[i][1], auto_ylims=True); 
    sm.graphics.tsa.plot_pacf(y, lags=np.arange(1,20), ax=ax[i][2], auto_ylims=True, method='ywm');

情報量基準

情報量基準は最尤法の推定結果を基に最適なモデルを選択する客観的な基準であり、一般的に
$$
IC=-2\mathcal{L}(\hat{\boldsymbol{\theta}})+p(T)k
$$
で定義される。ここでは応用で最も頻繁に用いられる2つの情報量基準を紹介する。
1つ目は赤池情報量基準(AIC) と呼ばれるものであり、
$$
AIC=-2\mathcal{L}(\hat{\boldsymbol{\theta}})+2k
$$で定義される。
もう1つはSchwarz情報量基準(SIC) もしくはベイズ情報量基準(BIC) と呼ばれるものであり、
$$
SIC=-2\mathcal{L}(\hat{\boldsymbol{\theta}})+\log(T)k
$$
で定義される。
真のモデルをAR($p$)過程とし、$p$が有限であるとすると、AICが選択する$\hat{p}$に関して
$$
\lim_{T→\infty} P(\hat{p}<p)=0\
\lim_{T→\infty} P(\hat{p}>p)>0
$$
が成立することが知られている。
つまり、標本数$T$が大きくなるとき、AICはモデル次数を過少に評価することはないが、過大なモデルを選択する確率が0とはならない。
AICはこの意味で一致性を持たない。それに対して、SICが選択する$\tilde{p}$に関しては
$$
\lim_{T→\infty} P(\tilde{p}=p)=1
$$
が成立し、SICは一致性を持つ。
AICとSICが異なるモデルを最適とした場合、どちらのモデルを採用するかは難しい問題である。

モデルの診断

真のモデルにおいては$\varepsilon_t$はホワイトノイズであるので、自己相関をもたない。
したがって、モデルが適切であれば、モデルから推定された$\hat{\varepsilon}_t$もほとんど自己相関を持たないはずである。
$\hat{\varepsilon}_t$の自己相関を検定することで、モデルの診断を行うことができる。
ARMA($p,q$)モデルを推定した場合、$m$次までの自己相関を用いたかばん検定の漸近分布は$\chi^2(m-p-q)$となる。
$Q(m)$の方が大きければ、帰無仮説を棄却する。帰無仮説が棄却された場合、モデルはあまり良いものとはいえないので、何らかの修正を考える必要がある。

AR過程についてAICを求めて診断を行う。

def calc_logL(y, k):
    """対数尤度
    """
    T = len(y)-k
    _y = np.array(y[k:])
    _x = np.array([y[i:i+k] for i in range(len(y)-k)])
    model = LinearRegression().fit(_x,_y)
    s2 = np.sum((_y-model.predict(_x))**2) / T
    return -T/2*np.log(2*np.pi*s2)-np.sum((_y-model.predict(_x))**2)/(2*s2)

def calc_AIC(y):
    """赤池情報量基準(AIC)
    """
    return np.array([-2*calc_logL(y, k)+2*(k+1) for k in range(1,10)])

T = 200
c = 1
sigma = .1
Phi = np.array([0.2,0.5])

y = make_ar2(c, Phi, sigma, T)

plt.plot(np.arange(1,10), calc_AIC(y), '-o')

statsmoedlで用意されている関数での診断は次のようになる。

selectmodel = ar_model.ar_select_order(y, 10)
display(selectmodel.bic)
print(selectmodel.ar_lags)

best_model = selectmodel.model.fit()
print(best_model.summary())

{(1, 2): -323.80776145370356,
 (1, 2, 3): -321.90955694436394,
 (1, 2, 3, 4): -316.7728912151294,
 (1, 2, 3, 4, 5): -311.9705892663207,
 (1, 2, 3, 4, 5, 6): -306.7411171455138,
 (1, 2, 3, 4, 5, 6, 7): -304.1614486813692,
 (1, 2, 3, 4, 5, 6, 7, 8): -299.2600822891575,
 (1, 2, 3, 4, 5, 6, 7, 8, 9): -294.0254086724308,
 (1, 2, 3, 4, 5, 6, 7, 8, 9, 10): -288.90507451328494,
 (1,): -277.49759575246037,
 0: -253.38001481990034}


[1, 2]
                            AutoReg Model Results                             
==============================================================================
Dep. Variable:                      y   No. Observations:                  202
Model:                     AutoReg(2)   Log Likelihood                 176.502
Method:               Conditional MLE   S.D. of innovations              0.100
Date:                Wed, 11 Jan 2023   AIC                           -345.004
Time:                        09:38:08   BIC                           -331.811
Sample:                             2   HQIC                          -339.665
                                  202                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const          1.0724      0.083     12.947      0.000       0.910       1.235
y.L1           0.1440      0.051      2.827      0.005       0.044       0.244
y.L2           0.5330      0.044     12.025      0.000       0.446       0.620
                                    Roots                                    
=============================================================================
                  Real          Imaginary           Modulus         Frequency
-----------------------------------------------------------------------------
AR.1            1.2413           +0.0000j            1.2413            0.0000
AR.2           -1.5114           +0.0000j            1.5114            0.5000
-----------------------------------------------------------------------------

ARMA過程についてはarma_order_select_ic関数を利用する。

c = 1
Phi=np.array([0.3,0.5, 0.1])
Theta=np.array([0.3,0.5])

y = make_arma(c=c, Phi=Phi, Theta=Theta, sigma=.2, t=300)

res_selection = sm.tsa.arma_order_select_ic(y, ic='aic', trend='n')
res_selection

{'aic':              0            1            2
 0  2214.461593  1830.335429  1479.019134
 1    64.876534    10.896186   -82.081280
 2   -47.795628   -67.470644  -105.906962
 3   -90.940931   -85.412068  -103.662765
 4  -109.595758   -90.944288  -106.563590,
 'aic_min_order': (4, 0)}

arima_model = sm.tsa.statespace.SARIMAX(y,
                             order=(3,0,2), 
                             enforce_stationarity = False, 
                             enforce_invertibility = False).fit()
print(arima_model.summary())

                               SARIMAX Results                                
==============================================================================
Dep. Variable:                      y   No. Observations:                  300
Model:               SARIMAX(3, 0, 2)   Log Likelihood                  69.723
Date:                Wed, 11 Jan 2023   AIC                           -127.446
Time:                        09:38:11   BIC                           -105.283
Sample:                             0   HQIC                          -118.573
                                - 300                                         
Covariance Type:                  opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ar.L1          0.9183      0.102      9.023      0.000       0.719       1.118
ar.L2          0.5279      0.090      5.889      0.000       0.352       0.704
ar.L3         -0.4460      0.088     -5.046      0.000      -0.619      -0.273
ma.L1         -0.3629      0.113     -3.205      0.001      -0.585      -0.141
ma.L2          0.1655      0.091      1.824      0.068      -0.012       0.343
sigma2         0.0366      0.003     10.670      0.000       0.030       0.043
===================================================================================
Ljung-Box (L1) (Q):                   0.02   Jarque-Bera (JB):                 4.60
Prob(Q):                              0.88   Prob(JB):                         0.10
Heteroskedasticity (H):               1.07   Skew:                            -0.25
Prob(H) (two-sided):                  0.73   Kurtosis:                         2.64
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).


C:\Users\00220900095\.virtualenvs\01.syouseiro-HheMhbJR\lib\site-packages\statsmodels\base\model.py:604: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
  warnings.warn("Maximum Likelihood optimization failed to "

演習のデータを見ていく。

import pandas as pd

def plot_y_acf_pacf(y, k=20):
    fig, ax = plt.subplots(1,3,figsize=(15,3))
    ax[0].plot(y)
    sm.graphics.tsa.plot_acf(y, lags=np.arange(1,k), ax=ax[1], auto_ylims=True); 
    sm.graphics.tsa.plot_pacf(y, lags=np.arange(1,k), ax=ax[2], auto_ylims=True, method='ywm');  

# 各種データの読み込み
df = pd.read_excel('economicdata.xls')
y = df['indprod'].values

# 対数差分系列の計算
y_rate = (np.log(y)[1:] - np.log(y)[:-1])*100
T = len(y_rate)

plot_y_acf_pacf(y_rate, k=20)

from statsmodels.stats.diagnostic import acorr_ljungbox

def check_criterion_and_plot_resid(y, PQ):
    print('(P,Q)   AIC    BIC')
    fig, ax = plt.subplots(6,3,figsize=(15,18))
    for i, (p,q) in enumerate(PQ):
        arima_model = sm.tsa.statespace.SARIMAX(y, 
                                 order=(p,0,q),
                                 enforce_stationarity = False, 
                                 enforce_invertibility = False).fit()
        print('('+str(p)+','+str(q)+')   '+str(np.round(arima_model.aic/T,2))+'   '+str(np.round(arima_model.bic/T,2)))
    
        # かばん検定の統計量の計算
        result = acorr_ljungbox(arima_model.resid, lags = 20)
        Q = result['lb_pvalue']
    
        sm.graphics.tsa.plot_acf(arima_model.resid, lags=np.arange(1,20), ax=ax[i][0], auto_ylims=True); 
        sm.graphics.tsa.plot_pacf(arima_model.resid, lags=np.arange(1,20), ax=ax[i][1], auto_ylims=True, method='ywm');    
        ax[i][2].plot(Q[:15])
        ax[i][2].hlines([0.05], 0,15, color='red', linestyle='--')
    
        ax[i][0].set_title('('+str(p)+','+str(q)+')', loc='left',fontdict={'fontsize':10})

PQ = np.array([[4,0],
               [0,3],
               [1,1],
               [2,1],
               [1,2],
               [2,2]])

check_criterion_and_plot_resid(y_rate, PQ)

(P,Q)   AIC    BIC
(4,0)   3.21   3.26
(0,3)   3.23   3.27
(1,1)   3.35   3.39
(2,1)   3.28   3.32
(1,2)   3.22   3.26
(2,2)   3.22   3.27

演習の2つ目のデータを見ていく。

df = pd.read_excel('arma.xls')

y1 = df['y1'].values
y2 = df['y2'].values
y3 = df['y3'].values

for y in [y1, y2, y3]:
    plot_y_acf_pacf(y, k=20)

PQ = np.array([[3,0],
               [0,6],
               [1,2],
               [1,3],
               [2,2],
               [2,3]])

check_criterion_and_plot_resid(y1, PQ)

(P,Q)   AIC    BIC
(3,0)   2.44   2.48
(0,6)   2.42   2.49
(1,2)   2.44   2.48
(1,3)   2.44   2.49
(2,2)   2.44   2.49
(2,3)   2.43   2.5

PQ = np.array([[2,0],
               [0,2],
               [1,1],
               [1,2],
               [2,1],
               [2,2]])

check_criterion_and_plot_resid(y2, PQ)

(P,Q)   AIC    BIC
(2,0)   2.41   2.44
(0,2)   2.42   2.45
(1,1)   2.39   2.43
(1,2)   2.39   2.43
(2,1)   2.4   2.44
(2,2)   2.4   2.45

PQ = np.array([[3,0],
               [0,2],
               [1,2],
               [1,3],
               [2,2],
               [2,3]])

check_criterion_and_plot_resid(y3, PQ)

(P,Q)   AIC    BIC
(3,0)   2.42   2.47
(0,2)   2.43   2.46
(1,2)   2.42   2.46
(1,3)   2.41   2.46
(2,2)   2.41   2.46
(2,3)   2.4   2.46

本参考書ではARIMAやSARIMAに関する記述がないためここで補足する。

補足：ARIMAモデル

自己回帰和分移動平均モデル(ARIMA) はARMAモデルを非定常過程に対応させたものである。
和分過程$I(d)$とは、$d-1$階差分を取った系列は非定常であるが、$d$階差分を取った系列が定常過程に従う過程である。

# 各種データの読み込み
df = pd.read_excel('economicdata.xls')
y = df['indprod'].values
T = len(y)

plot_y_acf_pacf(y, k=20)

yを加工せずにARMAモデルを作成し、残差の自己相関係数を確認する。
次数は前の演習に基づいて決めた。

arima_model = sm.tsa.statespace.SARIMAX(y, 
                                 order=(1,0,2),
                                 enforce_stationarity = False, 
                                 enforce_invertibility = False).fit()
plot_y_acf_pacf(arima_model.resid[2:], k=10)

1次の階差を取るARIMAモデルを作成し、残差の自己相関係数を確認する。
残差の自己相関係数はゼロに近づくことが確認できた。

arima_model = sm.tsa.statespace.SARIMAX(np.log(y), 
                                 order=(1,1,2),
                                 enforce_stationarity = False, 
                                 enforce_invertibility = False).fit()
plot_y_acf_pacf(arima_model.resid[2:], k=10)

補足：SARIMAモデル

季節変動自己回帰和分移動平均モデル(SARIMA) は階差を取ることで季節性を除去する操作をARMAモデルに加え、さらにARIMAモデルを適用したモデルである。

df = pd.read_csv('AirPassengers.csv')
y = df['#Passengers'].values
plot_y_acf_pacf(y, k=20)

まずはARIMAモデルで残差の確認をする。
12か月目のところで大きな相関が見られ、周期的な変化を除去できていないことがわかる。

arima_model = sm.tsa.statespace.SARIMAX(y, 
                                 order=(1,1,2),
                                 enforce_stationarity = False, 
                                 enforce_invertibility = False).fit()
plot_y_acf_pacf(arima_model.resid[2:], k=20)

SARIMAでモデルを作成すると、残差は有意にゼロとなり季節成分を除去できたことがわかる。
パラメータの調整などはここでは行わない。

sarima_model = sm.tsa.statespace.SARIMAX(y, 
                                 order=(1,1,2),
                                 seasonal_order=(1,1,1,12),
                                 enforce_stationarity = False, 
                                 enforce_invertibility = False).fit()
plot_y_acf_pacf(sarima_model.resid[2:], k=20)

以上となります。

次回

予測