pythonで「計量時系列分析」を読む2 ARMA過程

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

次の本を扱います。

### 参考書

https://www.amazon.co.jp/%E7%B5%8C%E6%B8%88%E3%83%BB%E3%83%95%E3%82%A1%E3%82%A4%E3%83%8A%E3%83%B3%E3%82%B9%E3%83%87%E3%83%BC%E3%82%BF%E3%81%AE%E8%A8%88%E9%87%8F%E6%99%82%E7%B3%BB%E5%88%97%E5%88%86%E6%9E%90-%E7%B5%B1%E8%A8%88%E3%83%A9%E3%82%A4%E3%83%96%E3%83%A9%E3%83%AA%E3%83%BC-%E6%B2%96%E6%9C%AC-%E7%AB%9C%E7%BE%A9/dp/4254127928

https://www.amazon.co.jp/%E6%99%82%E7%B3%BB%E5%88%97%E8%A7%A3%E6%9E%90-%E8%87%AA%E5%B7%B1%E5%9B%9E%E5%B8%B0%E5%9E%8B%E3%83%A2%E3%83%87%E3%83%AB%E3%83%BB%E7%8A%B6%E6%85%8B%E7%A9%BA%E9%96%93%E3%83%A2%E3%83%87%E3%83%AB%E3%83%BB%E7%95%B0%E5%B8%B8%E6%A4%9C%E7%9F%A5-Advanced-%E5%B3%B6%E7%94%B0-%E7%9B%B4%E5%B8%8C/dp/4320125010/ref=sr_1_2?__mk_ja_JP=%E3%82%AB%E3%82%BF%E3%82%AB%E3%83%8A&crid=2B3QRAXCAF0W8&keywords=%E6%99%82%E7%B3%BB%E5%88%97%E5%88%86%E6%9E%90+python&qid=1672965711&sprefix=%E6%99%82%E7%B3%BB%E5%88%97%E5%88%86%E6%9E%90+python%2Caps%2C180&sr=8-2

### 前回
時系列分析の基礎概念

https://qiita.com/tanaka_benkyo/items/2bcbb57be8ad38880a6c


### 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)過程の期待値は

```math
\begin{align}
E(y_t)&=E(\mu+\varepsilon_t+\theta_1\varepsilon_{t-1})\\
&=\mu
\end{align}
```
である。分散は、

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



```python
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
    

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/0c0c3d59-91be-b271-8df7-7fab965f94a5.png)

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


```python
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
    
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/3fa924f6-eeb9-ad36-e457-5c2880449767.png)


1次の自己共分散は

```math
\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次自己相関が

```math
\rho_1=\frac{\gamma_1}{\gamma_0}=\frac{\theta_1}{1+\theta_1^2}
```
で与えられることがわかる。2次以降の自己共分散は、$k\geq2$として、

```math
\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)過程の自己相関係数の計算を実装する。


```python
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
    

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/f19e40b9-b3ed-cc45-d613-3a81fbfa0319.png)



一般的に$q$次移動平均過程は、

```math
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.

```math
\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.
```math
\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$)過程の自己相関係数の計算を実装する。


```python
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));
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/4ef6e2d7-2149-b176-168b-5753b041f57d.png)


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


```python
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); 
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/d11004c5-4bdd-2d0d-fa34-6779473011d4.png)


#### 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)過程の自己相関の絶対値は指数的に減衰していくことがわかる。


```python
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
    

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/ac27c9b2-c21e-08b7-757b-4185bef95112.png)



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


```python
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));
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/d8985ce5-905d-db51-96b6-2e2a3a71fe11.png)


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


```python
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); 
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/931eb47e-98a8-8683-ecfb-17fdfe828f37.png)


$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$)過程は以下の性質を持つ
>1. $\mu=E(y_t)=\frac{c}{1-\phi_1-\phi_2-\cdots-\phi_p}$
>2. $\gamma_0=Var(y_t)=\frac{\sigma^2}{1-\phi_1\rho_1-\phi_2\rho_2-\cdots-\phi_p\rho_p}$
>3. 自己共分散と自己相関は$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式はユール・ウォーカー方程式と呼ばれる。
>4. 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$)についての自己相関係数を計算する。


```python
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));
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/ae838ada-8e47-6df1-641d-bb762febf927.png)



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


```python
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); 
```


![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/8c7f160f-06be-1519-17f3-660a3edb21c1.png)

    


#### 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$)過程は以下の性質を持つ。
>1. $\mu=E(y_t)=\frac{c}{1-\phi_1-\phi_2-\cdots-\phi_p}$
>2. $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}
$$
>3. ARMA過程の自己相関は指数的に減衰する

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


```python
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
    

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/a55ddb25-8196-dc2d-f7fb-96cb97397be3.png)


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


```python
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
    


![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/a27b5b49-2ced-1f10-e350-1dacd79905a9.png)


### 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)
$$
の定常条件は、

```math
\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)
$$
を繰り返し代入していくことで、

```math
\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}
$$
と書き直すことができる。

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


```python
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
    

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/cac5ef0b-f98c-eece-1a46-0eb8f49193d1.png)

    


#### 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過程で書き直すと、

```math
\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推定量が、

```math
\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推定量は以下の性質を持つ
>1. OLS推定量は**一致推定量**である
>2. OLS推定量を基準化したものは漸近的に正規分布に従う
>3. $\varepsilon_t\sim iid(0,\sigma^2)$のとき、OLS推定量は一致推定量の中で最小の分散共分散行列を持つ

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

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


```python
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);
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/b6a8ec70-431c-172a-79fa-a2e88946182e.png)



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


```python
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);
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/75708821-fa53-63c8-faaf-4f306684d79c.png)



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


```python
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を利用する。


```python
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過程の次数を決めることはできない。

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


```python
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')
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/c6d81955-76dd-95ff-ef86-f509a791573e.png)



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


```python
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'); 
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/b50fab7e-5bf5-934e-aa45-645569d39cee.png)



#### 情報量基準
**情報量基準**は最尤法の推定結果を基に最適なモデルを選択する客観的な基準であり、一般的に
$$
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を求めて診断を行う。


```python
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')
```



![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/f663f034-ee8d-a548-d50e-fe1007215f8e.png)



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


```python
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関数を利用する。


```python
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)}




```python
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 "
    

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


```python
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)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/d189a79f-93c8-2efb-8422-43d2c5cf3550.png)



```python
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
    

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/91a1cf58-f9dd-ae36-28ea-b01861320c30.png)

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/0396108e-237e-d37f-4140-23a5af5e7404.png)



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


```python
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)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/1c9226bf-5374-49fe-5f39-72c4f398e0d3.png)


```python
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
    

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/32c7ec36-7465-cb38-9789-61a924bcc4be.png)


![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/fc7b5ece-c0a9-74a8-20fd-35a5c03da43e.png)


```python
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
    

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/e20111b4-3887-e622-4190-cd02abe3629d.png)


```python
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
    

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/d92f30f6-1fb9-2a13-d55b-c81e471e195b.png)



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

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


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

plot_y_acf_pacf(y, k=20)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/a5cd6878-7088-0311-65f5-9f8859533e01.png)



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


```python
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)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/d9207694-024f-88f0-a655-6e8a02324f23.png)



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


```python
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)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/ffde319c-4415-53f1-7454-e9a3781143a7.png)


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


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

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/b7f31b05-9282-9fc6-ef86-25d43828de7a.png)



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


```python
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)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/78df5061-3e43-fcb9-91ae-7504a077d893.png)


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


```python
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)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/184f30c5-1b18-043e-e98e-e0c086c2d87f.png)

以上となります。

### 次回

予測