【第2章】pythonで「経済・ファイナンスデータの計量時系列分析」の章末問題を解く

自分の勉強用に、[沖本竜義「経済・ファイナンスデータの計量経済分析」](https://www.asakura.co.jp/G_12.php?isbn=ISBN978-4-254-12792-8)の章末問題をpythonで解いてみました。
jupyter notebook上で記述したものをほぼそのまま載せてあります。
ソースコードの細かい説明は省いてありますが、今後余裕があれば説明を加えようと思います。

こちらの記事は第2章の章末問題を解いたものになります。
[第1章](https://qiita.com/mckeeeen/items/666e14c94e1548130646#15-5)、[第4章](https://qiita.com/mckeeeen/items/4afa1167008c1f315d0b)、[第5章](https://qiita.com/mckeeeen/items/d293f8005fee5452ecbb)、[第6章](https://qiita.com/mckeeeen/items/55cb0130e5a7fa035e1e)、[第7章](https://qiita.com/mckeeeen/items/ce000b00b0c725503072)も公開しています。



各種設定・モジュールのインポート


```python
%matplotlib inline
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import chi2
from statsmodels.tsa.stattools import acf
from statsmodels.tsa.stattools import pacf
from statsmodels.stats.diagnostic import acorr_ljungbox
from statsmodels.tsa.arima_model import ARMA

import matplotlib as mpl
font = {"family":"IPAexGothic"}
mpl.rc('font', **font)
plt.rcParams["font.size"] = 12
```




# 2.5 
## (1)

ファイルをインポートする


```python
eco = pd.read_csv('./input/economicdata.csv')
print(eco.shape)
eco.head()
```

    (364, 7)





<div>

<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>date</th>
      <th>topix</th>
      <th>exrate</th>
      <th>indprod</th>
      <th>cpi</th>
      <th>saunemp</th>
      <th>intrate</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>Jan-75</td>
      <td>276.09</td>
      <td>29.13</td>
      <td>47.33</td>
      <td>52.625</td>
      <td>1.7</td>
      <td>12.67</td>
    </tr>
    <tr>
      <th>1</th>
      <td>Feb-75</td>
      <td>299.81</td>
      <td>29.70</td>
      <td>46.86</td>
      <td>52.723</td>
      <td>1.8</td>
      <td>13.00</td>
    </tr>
    <tr>
      <th>2</th>
      <td>Mar-75</td>
      <td>313.50</td>
      <td>29.98</td>
      <td>46.24</td>
      <td>53.114</td>
      <td>1.8</td>
      <td>12.92</td>
    </tr>
    <tr>
      <th>3</th>
      <td>Apr-75</td>
      <td>320.57</td>
      <td>29.80</td>
      <td>47.33</td>
      <td>54.092</td>
      <td>1.8</td>
      <td>12.02</td>
    </tr>
    <tr>
      <th>4</th>
      <td>May-75</td>
      <td>329.65</td>
      <td>29.79</td>
      <td>47.33</td>
      <td>54.385</td>
      <td>1.8</td>
      <td>11.06</td>
    </tr>
  </tbody>
</table>
</div>



対数差分系列を計算する


```python
eco['indprod_dlog'] = np.log(eco.indprod).diff()
eco.head()
```




<div>

<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>date</th>
      <th>topix</th>
      <th>exrate</th>
      <th>indprod</th>
      <th>cpi</th>
      <th>saunemp</th>
      <th>intrate</th>
      <th>indprod_dlog</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>Jan-75</td>
      <td>276.09</td>
      <td>29.13</td>
      <td>47.33</td>
      <td>52.625</td>
      <td>1.7</td>
      <td>12.67</td>
      <td>NaN</td>
    </tr>
    <tr>
      <th>1</th>
      <td>Feb-75</td>
      <td>299.81</td>
      <td>29.70</td>
      <td>46.86</td>
      <td>52.723</td>
      <td>1.8</td>
      <td>13.00</td>
      <td>-0.009980</td>
    </tr>
    <tr>
      <th>2</th>
      <td>Mar-75</td>
      <td>313.50</td>
      <td>29.98</td>
      <td>46.24</td>
      <td>53.114</td>
      <td>1.8</td>
      <td>12.92</td>
      <td>-0.013319</td>
    </tr>
    <tr>
      <th>3</th>
      <td>Apr-75</td>
      <td>320.57</td>
      <td>29.80</td>
      <td>47.33</td>
      <td>54.092</td>
      <td>1.8</td>
      <td>12.02</td>
      <td>0.023299</td>
    </tr>
    <tr>
      <th>4</th>
      <td>May-75</td>
      <td>329.65</td>
      <td>29.79</td>
      <td>47.33</td>
      <td>54.385</td>
      <td>1.8</td>
      <td>11.06</td>
      <td>0.000000</td>
    </tr>
  </tbody>
</table>
</div>



標本自己相関関数と標本偏自己相関関数を計算する
ただし、標本偏自己相関関数はOLS係数推定量として求める（教科書に沿った形です）


```python
acf_ind=acf(eco.indprod_dlog[1:]) 
pacf_ind=pacf(eco.indprod_dlog[1:],method='ols') 
lag=range(41)
```

プロットする


```python
y=[acf_ind,pacf_ind]
fig,ax=plt.subplots(1,2,figsize=[14,6])
for i in range(2):   
    ax[i].bar(lag[1:21],y[i][1:21])
    ax[i].set_ylim(-0.5,0.5)
    ax[i].set_xticks([1,5,10,15,20])
    ax[i].set_yticks(np.asarray(range(-5,6))*0.1)
    ax[i].hlines(0,0,21,linewidth=1)
    ax[i].hlines([1.96/np.sqrt(eco.shape[0]-1),-1.96/np.sqrt(eco.shape[0]-1)],0,21,linestyles='dashed',linewidth=1)
    ax[i].set_xlabel('ラグ',fontsize=12)
    ax[i].set_ylabel('自己相関',fontsize=12)
ax[0].set_title('標本自己相関（鉱工業生産指数の成長率）')
ax[1].set_title('偏標本自己相関（鉱工業生産指数の成長率）')

plt.show()
```


![output_10_0.png](https://qiita-image-store.s3.amazonaws.com/0/235450/c363d29e-c248-6126-5fad-19243d88afaf.png)



AR(4),MA(3),ARMA(1,1),ARMA(2,1),ARMA(1,2),ARMA(2,2)でモデル化する
ただし、statsmodelsにはMAモデルのクラスが見当たらなかったのでARMAモデルのAR成分を0にすることで対応
今回は、ARモデルもARMAモデルのMA成分を0にすることで対応してみる
参考: http://www.blackarbs.com/blog/time-series-analysis-in-python-linear-models-to-garch/11/1/2016#AR


```python
ar_4=ARMA(eco.indprod_dlog[1:].values,[4,0]).fit()
ma_3=ARMA(eco.indprod_dlog[1:].values,[0,3]).fit()
arma_1_1=ARMA(eco.indprod_dlog[1:].values,[1,1]).fit()
arma_2_1=ARMA(eco.indprod_dlog[1:].values,[2,1]).fit()
arma_1_2=ARMA(eco.indprod_dlog[1:].values,[1,2]).fit()
arma_2_2=ARMA(eco.indprod_dlog[1:].values,[2,2]).fit()
```

各モデルのAIC,SICを計算する
教科書の値と異なる値となったが、選択されるモデルは同じ


```python
model=[ar_4,ma_3,arma_1_1,arma_2_1,arma_1_2,arma_2_2]
aic=[]
sic=[]
T=eco.indprod_dlog[1:].shape[0]
for i in range(len(model)):
    aic.append(model[i].aic/T)
    sic.append(model[i].bic/T)
model_sel=pd.DataFrame({'AIC':aic,'SIC':sic},index=['AR(4)','MA(3)','ARMA(1,1)','ARMA(2,1)','ARMA(1,2)','ARMA(2,2)'])
model_sel.T
```




<div>

<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>AR(4)</th>
      <th>MA(3)</th>
      <th>ARMA(1,1)</th>
      <th>ARMA(2,1)</th>
      <th>ARMA(1,2)</th>
      <th>ARMA(2,2)</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>AIC</th>
      <td>-5.969405</td>
      <td>-5.961817</td>
      <td>-5.876225</td>
      <td>-5.923565</td>
      <td>-5.968149</td>
      <td>-5.954818</td>
    </tr>
    <tr>
      <th>SIC</th>
      <td>-5.905035</td>
      <td>-5.908175</td>
      <td>-5.833311</td>
      <td>-5.869923</td>
      <td>-5.914507</td>
      <td>-5.890448</td>
    </tr>
  </tbody>
</table>
</div>



AR(4)モデルとARMA(1,2)モデルの残差を計算し、標本自己相関関数を計算する


```python
acf_ar_4=acf(ar_4.resid)
acf_arma_1_2=acf(arma_1_2.resid)
```

プロットする


```python
y=[acf_ar_4,acf_arma_1_2]
fig,ax=plt.subplots(1,2,figsize=[14,6])
for i in range(2):   
    ax[i].bar(lag[1:21],y[i][1:21])
    ax[i].set_ylim(-0.5,0.5)
    ax[i].set_xticks([1,5,10,15,20])
    ax[i].set_yticks(np.asarray(range(-5,6))*0.1)
    ax[i].hlines(0,0,21,linewidth=1)
    ax[i].set_xlabel('ラグ',fontsize=12)
    ax[i].set_ylabel('自己相関',fontsize=12)
ax[0].hlines([1.96/np.sqrt(ar_4.resid.shape[0]),-1.96/np.sqrt(ar_4.resid.shape[0])],0,21,linestyles='dashed',linewidth=1)
ax[1].hlines([1.96/np.sqrt(arma_1_2.resid.shape[0]),-1.96/np.sqrt(arma_1_2.resid.shape[0])],0,21,linestyles='dashed',linewidth=1)
ax[0].set_title('標本自己相関（AR(4)モデルの残差）')
ax[1].set_title('標本自己相関（ARMA(1,2)モデルの残差）')

plt.show()
```


![output_18_0.png](https://qiita-image-store.s3.amazonaws.com/0/235450/4b7f6fe4-855e-f3d8-e22b-7ff7957f3b2a.png)



## (2)

AR(4)モデルの残差について自己相関検定を行う

Q(m)やP値が教科書と異なるのは、学習に用いたmethodが異なるからだと考えられる（残差も正確には教科書と異なっている）

結果、AR(4)モデルの残差が自己相関を持つとは言えない（適切なモデルであると言える）


```python
q_m_ar_4=acorr_ljungbox(ar_4.resid)[0]
p_value_ar_4=chi2.sf(q_m_ar_4,6)
lag=range(1,41)
test_ar_4=pd.DataFrame({'lag':lag,'autocorr':acf_ar_4[1:],'Q(m)':q_m_ar_4,'P_value':p_value_ar_4})
test_ar_4[['lag','autocorr','Q(m)','P_value']].iloc[0:10,:].T
```




<div>

<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>0</th>
      <th>1</th>
      <th>2</th>
      <th>3</th>
      <th>4</th>
      <th>5</th>
      <th>6</th>
      <th>7</th>
      <th>8</th>
      <th>9</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>lag</th>
      <td>1.000000</td>
      <td>2.000000</td>
      <td>3.000000</td>
      <td>4.000000</td>
      <td>5.000000</td>
      <td>6.000000</td>
      <td>7.000000</td>
      <td>8.000000</td>
      <td>9.000000</td>
      <td>10.000000</td>
    </tr>
    <tr>
      <th>autocorr</th>
      <td>-0.009783</td>
      <td>-0.017589</td>
      <td>-0.004028</td>
      <td>0.012414</td>
      <td>0.027854</td>
      <td>-0.007964</td>
      <td>0.031076</td>
      <td>0.127591</td>
      <td>0.058008</td>
      <td>-0.080286</td>
    </tr>
    <tr>
      <th>Q(m)</th>
      <td>0.035029</td>
      <td>0.148575</td>
      <td>0.154546</td>
      <td>0.211425</td>
      <td>0.498562</td>
      <td>0.522104</td>
      <td>0.881525</td>
      <td>6.957480</td>
      <td>8.216885</td>
      <td>10.636289</td>
    </tr>
    <tr>
      <th>P_value</th>
      <td>0.999999</td>
      <td>0.999935</td>
      <td>0.999927</td>
      <td>0.999818</td>
      <td>0.997856</td>
      <td>0.997559</td>
      <td>0.989708</td>
      <td>0.324797</td>
      <td>0.222641</td>
      <td>0.100289</td>
    </tr>
  </tbody>
</table>
</div>



ARMA(1,2)モデルの残差について自己相関検定を行う
Q(m)やP値が教科書と異なるのは、学習に用いたmethodが異なるからだと考えられる（残差も正確には教科書と異なっている）
結果、AR(4)モデルの残差が自己相関を持つとは言えない（適切なモデルであると言える）


```python
q_m_arma_1_2=acorr_ljungbox(arma_1_2.resid)[0]
p_value_arma_1_2=chi2.sf(q_m_arma_1_2,7)
lag=range(1,41)
test_arma_1_2=pd.DataFrame({'lag':lag,'autocorr':acf_arma_1_2[1:],'Q(m)':q_m_arma_1_2,'P_value':p_value_arma_1_2})
test_arma_1_2[['lag','autocorr','Q(m)','P_value']].iloc[0:10,:].T
```




<div>

<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>0</th>
      <th>1</th>
      <th>2</th>
      <th>3</th>
      <th>4</th>
      <th>5</th>
      <th>6</th>
      <th>7</th>
      <th>8</th>
      <th>9</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>lag</th>
      <td>1.000000</td>
      <td>2.000000</td>
      <td>3.000000</td>
      <td>4.000000</td>
      <td>5.000000</td>
      <td>6.000000</td>
      <td>7.000000</td>
      <td>8.000000</td>
      <td>9.000000</td>
      <td>10.000000</td>
    </tr>
    <tr>
      <th>autocorr</th>
      <td>-0.011002</td>
      <td>0.005265</td>
      <td>0.060629</td>
      <td>-0.034401</td>
      <td>-0.070794</td>
      <td>0.010740</td>
      <td>0.020772</td>
      <td>0.116308</td>
      <td>0.054963</td>
      <td>-0.062463</td>
    </tr>
    <tr>
      <th>Q(m)</th>
      <td>0.044306</td>
      <td>0.054481</td>
      <td>1.407365</td>
      <td>1.844126</td>
      <td>3.698984</td>
      <td>3.741790</td>
      <td>3.902374</td>
      <td>8.951175</td>
      <td>10.081840</td>
      <td>11.546259</td>
    </tr>
    <tr>
      <th>P_value</th>
      <td>1.000000</td>
      <td>1.000000</td>
      <td>0.985345</td>
      <td>0.967962</td>
      <td>0.813722</td>
      <td>0.808991</td>
      <td>0.790951</td>
      <td>0.256181</td>
      <td>0.183983</td>
      <td>0.116503</td>
    </tr>
  </tbody>
</table>
</div>



# 2.6
## (1)

ファイルをインポートする


```python
arma = pd.read_csv('./input/arma.csv')
print(arma.shape)
arma.head()
```

    (300, 3)





<div>

<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>y1</th>
      <th>y2</th>
      <th>y3</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>-0.091682</td>
      <td>0.104645</td>
      <td>-1.350264</td>
    </tr>
    <tr>
      <th>1</th>
      <td>0.533702</td>
      <td>1.253242</td>
      <td>-0.324967</td>
    </tr>
    <tr>
      <th>2</th>
      <td>-0.501343</td>
      <td>0.189352</td>
      <td>0.291934</td>
    </tr>
    <tr>
      <th>3</th>
      <td>0.486156</td>
      <td>-0.189821</td>
      <td>-1.472085</td>
    </tr>
    <tr>
      <th>4</th>
      <td>0.489030</td>
      <td>0.451640</td>
      <td>-1.335709</td>
    </tr>
  </tbody>
</table>
</div>



y1の標本自己相関関数と標本偏自己相関関数を計算する


```python
acf_y1=acf(arma.y1)
pacf_y1=pacf(arma.y1)
```

プロットする


```python
y=[acf_y1,pacf_y1]
fig,ax=plt.subplots(1,2,figsize=[14,6])
for i in range(2):   
    ax[i].bar(range(1,21),y[i][1:21])
    ax[i].set_ylim(-1,1)
    ax[i].set_xticks([1,5,10,15,20])
    ax[i].set_yticks(np.asarray(range(-10,11))*0.1)
    ax[i].hlines(0,0,21,linewidth=1)
    ax[i].hlines([1.96/np.sqrt(arma.shape[0]),-1.96/np.sqrt(arma.shape[0])],0,21,linestyles='dashed',linewidth=1)
    ax[i].set_xlabel('ラグ',fontsize=12)
    ax[i].set_ylabel('自己相関',fontsize=12)
ax[0].set_title('標本自己相関（y1）')
ax[1].set_title('偏標本自己相関（y1）')

plt.show()
```


![output_30_0.png](https://qiita-image-store.s3.amazonaws.com/0/235450/dd1bca9c-0322-17c6-be8b-8e95188bf8ba.png)



## (2)

上記プロットの結果から、AR(2)が最有力で、ARMAの可能性もある
AR(2)、ARMA(1,1)、ARMA(2,1)、ARMA(1,2)をモデル候補とする（ARMA(2,2)は定常性の過程が成り立たないとして却下された）


```python
ar_2=ARMA(arma.y1.values,[2,0]).fit()
arma_1_1=ARMA(arma.y1.values,[1,1]).fit()
arma_2_1=ARMA(arma.y1.values,[2,1]).fit()
arma_1_2=ARMA(arma.y1.values,[1,2]).fit()
# arma_2_2=ARMA(arma.y1.values,[2,2]).fit()
```


```python
model=[ar_2,arma_1_1,arma_2_1,arma_1_2]
aic=[]
sic=[]
T=eco.indprod_dlog[1:].shape[0]
for i in range(len(model)):
    aic.append(model[i].aic/T)
    sic.append(model[i].bic/T)
model_sel=pd.DataFrame({'AIC':aic,'SIC':sic},index=['AR(2)','ARMA(1,1)','ARMA(2,1)','ARMA(1,2)'])
model_sel.T
```




<div>

<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>AR(2)</th>
      <th>ARMA(1,1)</th>
      <th>ARMA(2,1)</th>
      <th>ARMA(1,2)</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>AIC</th>
      <td>2.499130</td>
      <td>2.507680</td>
      <td>2.497055</td>
      <td>2.503431</td>
    </tr>
    <tr>
      <th>SIC</th>
      <td>2.539943</td>
      <td>2.548493</td>
      <td>2.548072</td>
      <td>2.554448</td>
    </tr>
  </tbody>
</table>
</div>



## (4)

AICではARMA(2,1)、SICではAR(2)が選択される

まず、ARMA(2,1)モデルの残差について自己相関検定を行う（ただし、m=10とする）

結果、ARMA(2,1)モデルの残差が自己相関を持つとは言えない（適切なモデルであると言える）


```python
q_m_arma_2_1=acorr_ljungbox(arma_2_1.resid)[0]
p_value_arma_2_1=chi2.sf(q_m_arma_2_1,7)
lag=range(1,41)
test_arma_2_1=pd.DataFrame({'lag':lag,'autocorr':acf_y1[1:],'Q(m)':q_m_arma_2_1,'P_value':p_value_arma_2_1})
test_arma_2_1[['lag','autocorr','Q(m)','P_value']].iloc[0:10,:].T
```




<div>

<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>0</th>
      <th>1</th>
      <th>2</th>
      <th>3</th>
      <th>4</th>
      <th>5</th>
      <th>6</th>
      <th>7</th>
      <th>8</th>
      <th>9</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>lag</th>
      <td>1.000000</td>
      <td>2.000000</td>
      <td>3.000000</td>
      <td>4.000000</td>
      <td>5.000000</td>
      <td>6.000000</td>
      <td>7.000000</td>
      <td>8.000000</td>
      <td>9.000000</td>
      <td>10.000000</td>
    </tr>
    <tr>
      <th>autocorr</th>
      <td>0.830991</td>
      <td>0.761174</td>
      <td>0.652177</td>
      <td>0.584840</td>
      <td>0.503205</td>
      <td>0.450812</td>
      <td>0.383598</td>
      <td>0.328185</td>
      <td>0.278296</td>
      <td>0.229232</td>
    </tr>
    <tr>
      <th>Q(m)</th>
      <td>0.002502</td>
      <td>0.036143</td>
      <td>0.112850</td>
      <td>0.128471</td>
      <td>0.190753</td>
      <td>0.915685</td>
      <td>0.920768</td>
      <td>0.941598</td>
      <td>1.040412</td>
      <td>3.827893</td>
    </tr>
    <tr>
      <th>P_value</th>
      <td>1.000000</td>
      <td>1.000000</td>
      <td>0.999996</td>
      <td>0.999995</td>
      <td>0.999979</td>
      <td>0.996076</td>
      <td>0.996007</td>
      <td>0.995716</td>
      <td>0.994150</td>
      <td>0.799374</td>
    </tr>
  </tbody>
</table>
</div>



次に、AR(2)モデルの残差について自己相関検定を行う（ただし、m=10とする）

結果、AR(2)モデルの残差が自己相関を持つとは言えない（適切なモデルであると言える）


```python
q_m_ar_2=acorr_ljungbox(ar_2.resid)[0]
p_value_ar_2=chi2.sf(q_m_ar_2,8)
lag=range(1,41)
test_ar_2=pd.DataFrame({'lag':lag,'autocorr':acf_y1[1:],'Q(m)':q_m_ar_2,'P_value':p_value_ar_2})
test_ar_2[['lag','autocorr','Q(m)','P_value']].iloc[0:10,:].T
```




<div>

<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>0</th>
      <th>1</th>
      <th>2</th>
      <th>3</th>
      <th>4</th>
      <th>5</th>
      <th>6</th>
      <th>7</th>
      <th>8</th>
      <th>9</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>lag</th>
      <td>1.000000</td>
      <td>2.000000</td>
      <td>3.000000</td>
      <td>4.000000</td>
      <td>5.000000</td>
      <td>6.000000</td>
      <td>7.000000</td>
      <td>8.000000</td>
      <td>9.000000</td>
      <td>10.000000</td>
    </tr>
    <tr>
      <th>autocorr</th>
      <td>0.830991</td>
      <td>0.761174</td>
      <td>0.652177</td>
      <td>0.584840</td>
      <td>0.503205</td>
      <td>0.450812</td>
      <td>0.383598</td>
      <td>0.328185</td>
      <td>0.278296</td>
      <td>0.229232</td>
    </tr>
    <tr>
      <th>Q(m)</th>
      <td>0.138707</td>
      <td>1.082496</td>
      <td>2.001496</td>
      <td>2.244934</td>
      <td>2.769646</td>
      <td>3.477106</td>
      <td>3.482886</td>
      <td>3.487488</td>
      <td>3.695579</td>
      <td>6.845535</td>
    </tr>
    <tr>
      <th>P_value</th>
      <td>0.999999</td>
      <td>0.997672</td>
      <td>0.980966</td>
      <td>0.972567</td>
      <td>0.947971</td>
      <td>0.900959</td>
      <td>0.900514</td>
      <td>0.900159</td>
      <td>0.883492</td>
      <td>0.553386</td>
    </tr>
  </tbody>
</table>
</div>



## (5)

まず、y2に関して同様の分析を行う

y2の標本自己相関関数と標本偏自己相関関数を計算する


```python
acf_y2=acf(arma.y2)
pacf_y2=pacf(arma.y2)
```




プロットする


```python
y=[acf_y2,pacf_y2]
fig,ax=plt.subplots(1,2,figsize=[14,6])
for i in range(2):   
    ax[i].bar(range(1,21),y[i][1:21])
    ax[i].set_ylim(-1,1)
    ax[i].set_xticks([1,5,10,15,20])
    ax[i].set_yticks(np.asarray(range(-10,11))*0.1)
    ax[i].hlines(0,0,21,linewidth=1)
    ax[i].hlines([1.96/np.sqrt(arma.shape[0]),-1.96/np.sqrt(arma.shape[0])],0,21,linestyles='dashed',linewidth=1)
    ax[i].set_xlabel('ラグ',fontsize=12)
    ax[i].set_ylabel('自己相関',fontsize=12)
ax[0].set_title('標本自己相関（y2）')
ax[1].set_title('偏標本自己相関（y2）')

plt.show()
```

![output_45_0.png](https://qiita-image-store.s3.amazonaws.com/0/235450/d84a6f3c-a1b6-696a-36ce-17552c7abd0f.png)


上記プロットの結果から、MA(2)、MA(3)が有力で、ARMAの可能性もある
季節性が残っているようにも見える
MA(2)、ARMA(2,1)、ARMA(2,2)、ARMA(2,3)、ARMA(3,1)、ARMA(3,2)、ARMA(3,3)をモデル候補とする（MA(3)、ARMA(1,1)、ARMA(1,2)、ARMA(1,3)は定常性の過程が成り立たないとして却下された）


```python
ma_2=ARMA(arma.y2.values,[0,2]).fit()
# ma_3=ARMA(arma.y2.values,[0,3]).fit()
# arma_1_1=ARMA(arma.y2.values,[1,1]).fit()
# arma_1_2=ARMA(arma.y2.values,[1,2]).fit()
# arma_1_3=ARMA(arma.y2.values,[1,3]).fit()
arma_2_1=ARMA(arma.y2.values,[2,1]).fit()
arma_2_2=ARMA(arma.y2.values,[2,2]).fit()
arma_2_3=ARMA(arma.y2.values,[2,3]).fit()
arma_3_1=ARMA(arma.y2.values,[3,1]).fit()
arma_3_2=ARMA(arma.y2.values,[3,2]).fit()
arma_3_3=ARMA(arma.y2.values,[3,3]).fit()
```


```python
model=[ma_2,arma_2_1,arma_2_2,arma_2_3,arma_3_1,arma_3_2,arma_3_3]
aic=[]
sic=[]
T=eco.indprod_dlog[1:].shape[0]
for i in range(len(model)):
    aic.append(model[i].aic/T)
    sic.append(model[i].bic/T)
model_sel=pd.DataFrame({'AIC':aic,'SIC':sic},index=['MA(2)','ARMA(2,1)','ARMA(2,2)','ARMA(2,3)','ARMA(3,1)','ARMA(3,2)','ARMA(3,3)'])
model_sel.T
```




<div>

<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>MA(2)</th>
      <th>ARMA(2,1)</th>
      <th>ARMA(2,2)</th>
      <th>ARMA(2,3)</th>
      <th>ARMA(3,1)</th>
      <th>ARMA(3,2)</th>
      <th>ARMA(3,3)</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>AIC</th>
      <td>2.356034</td>
      <td>2.293800</td>
      <td>2.298995</td>
      <td>2.296927</td>
      <td>2.29918</td>
      <td>2.302264</td>
      <td>2.302285</td>
    </tr>
    <tr>
      <th>SIC</th>
      <td>2.396847</td>
      <td>2.344816</td>
      <td>2.360215</td>
      <td>2.368350</td>
      <td>2.36040</td>
      <td>2.373687</td>
      <td>2.383911</td>
    </tr>
  </tbody>
</table>
</div>



AICとSICどちらでもARMA(2,1)が選択される
ARMA(2,1)モデルの残差について自己相関検定を行う（ただし、m=10とする）

結果、ARMA(2,1)モデルの残差が自己相関を持つとは言えない（適切なモデルであると言える）


```python
q_m_arma_2_1=acorr_ljungbox(arma_2_1.resid)[0]
p_value_arma_2_1=chi2.sf(q_m_arma_2_1,7)
lag=range(1,41)
test_arma_2_1=pd.DataFrame({'lag':lag,'autocorr':acf_y2[1:],'Q(m)':q_m_arma_2_1,'P_value':p_value_arma_2_1})
test_arma_2_1[['lag','autocorr','Q(m)','P_value']].iloc[0:10,:].T
```




<div>

<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>0</th>
      <th>1</th>
      <th>2</th>
      <th>3</th>
      <th>4</th>
      <th>5</th>
      <th>6</th>
      <th>7</th>
      <th>8</th>
      <th>9</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>lag</th>
      <td>1.000000</td>
      <td>2.000000</td>
      <td>3.000000</td>
      <td>4.000000</td>
      <td>5.000000</td>
      <td>6.000000</td>
      <td>7.000000</td>
      <td>8.000000</td>
      <td>9.000000</td>
      <td>10.000000</td>
    </tr>
    <tr>
      <th>autocorr</th>
      <td>0.301309</td>
      <td>-0.330950</td>
      <td>-0.110256</td>
      <td>0.026641</td>
      <td>0.001279</td>
      <td>-0.001494</td>
      <td>0.013197</td>
      <td>0.035491</td>
      <td>0.036527</td>
      <td>-0.034517</td>
    </tr>
    <tr>
      <th>Q(m)</th>
      <td>0.000551</td>
      <td>0.236224</td>
      <td>0.296538</td>
      <td>2.742206</td>
      <td>2.891474</td>
      <td>2.917825</td>
      <td>3.211499</td>
      <td>3.232176</td>
      <td>3.393568</td>
      <td>3.562742</td>
    </tr>
    <tr>
      <th>P_value</th>
      <td>1.000000</td>
      <td>0.999956</td>
      <td>0.999904</td>
      <td>0.907783</td>
      <td>0.894857</td>
      <td>0.892499</td>
      <td>0.864772</td>
      <td>0.862726</td>
      <td>0.846367</td>
      <td>0.828536</td>
    </tr>
  </tbody>
</table>
</div>



y3の標本自己相関関数と標本偏自己相関関数を計算する


```python
acf_y3=acf(arma.y3)
pacf_y3=pacf(arma.y3)
```




プロットする


```python
y=[acf_y3,pacf_y3]
fig,ax=plt.subplots(1,2,figsize=[14,6])
for i in range(2):   
    ax[i].bar(range(1,21),y[i][1:21])
    ax[i].set_ylim(-1,1)
    ax[i].set_xticks([1,5,10,15,20])
    ax[i].set_yticks(np.asarray(range(-10,11))*0.1)
    ax[i].hlines(0,0,21,linewidth=1)
    ax[i].hlines([1.96/np.sqrt(arma.shape[0]),-1.96/np.sqrt(arma.shape[0])],0,21,linestyles='dashed',linewidth=1)
    ax[i].set_xlabel('ラグ',fontsize=12)
    ax[i].set_ylabel('自己相関',fontsize=12)
ax[0].set_title('標本自己相関（y3）')
ax[1].set_title('偏標本自己相関（y3）')

plt.show()
```


![output_54_0.png](https://qiita-image-store.s3.amazonaws.com/0/235450/1d2356a7-edb9-2208-39f5-667723f74a68.png)



上記プロットの結果から、季節性があり、ARMAは当てはまらなさそう
ARMA(1,1)、ARMA(1,2)、ARMA(1,3)、ARMA(1,4)、ARMA(2,1)、ARMA(2,2)、ARMA(2,3)、ARMA(2,4)、ARMA(3,1)、ARMA(3,2)、ARMA(3,3)、ARMA(3,4)、ARMA(4,1)、ARMA(4,2)、ARMA(4,3)、ARMA(4,4)をモデル候補とする


```python
arma_1_1=ARMA(arma.y3.values,[1,1]).fit()
arma_1_2=ARMA(arma.y3.values,[1,2]).fit()
arma_1_3=ARMA(arma.y3.values,[1,3]).fit()
arma_1_4=ARMA(arma.y3.values,[1,4]).fit()
arma_2_1=ARMA(arma.y3.values,[2,1]).fit()
arma_2_2=ARMA(arma.y3.values,[2,2]).fit()
arma_2_3=ARMA(arma.y3.values,[2,3]).fit()
arma_2_4=ARMA(arma.y3.values,[2,4]).fit()
arma_3_1=ARMA(arma.y3.values,[3,1]).fit()
arma_3_2=ARMA(arma.y3.values,[3,2]).fit()
arma_3_3=ARMA(arma.y3.values,[3,3]).fit()
arma_3_4=ARMA(arma.y3.values,[3,4]).fit()
arma_4_1=ARMA(arma.y3.values,[4,1]).fit()
arma_4_2=ARMA(arma.y3.values,[4,2]).fit()
arma_4_3=ARMA(arma.y3.values,[4,3]).fit()
arma_4_4=ARMA(arma.y3.values,[4,4]).fit()
```





```python
model=[arma_1_1,arma_1_2,arma_1_3,arma_1_4,arma_2_1,arma_2_2,arma_2_3,arma_2_4,arma_3_1,arma_3_2,arma_3_3,arma_3_4,arma_4_1,arma_4_2,arma_4_3,arma_4_4]
aic=[]
sic=[]
T=eco.indprod_dlog[1:].shape[0]
for i in range(len(model)):
    aic.append(model[i].aic/T)
    sic.append(model[i].bic/T)
model_sel=pd.DataFrame({'AIC':aic,'SIC':sic},index=['ARMA(1,1)','ARMA(1,2)','ARMA(1,3)','ARMA(1,4)','ARMA(2,1)','ARMA(2,2)','ARMA(2,3)','ARMA(2,4)','ARMA(3,1)','ARMA(3,2)','ARMA(3,3)','ARMA(4,4)','ARMA(4,1)','ARMA(4,2)','ARMA(4,3)','ARMA(4,4)'])
model_sel.T
```




<div>

<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>ARMA(1,1)</th>
      <th>ARMA(1,2)</th>
      <th>ARMA(1,3)</th>
      <th>ARMA(1,4)</th>
      <th>ARMA(2,1)</th>
      <th>ARMA(2,2)</th>
      <th>ARMA(2,3)</th>
      <th>ARMA(2,4)</th>
      <th>ARMA(3,1)</th>
      <th>ARMA(3,2)</th>
      <th>ARMA(3,3)</th>
      <th>ARMA(4,4)</th>
      <th>ARMA(4,1)</th>
      <th>ARMA(4,2)</th>
      <th>ARMA(4,3)</th>
      <th>ARMA(4,4)</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>AIC</th>
      <td>2.537023</td>
      <td>2.533205</td>
      <td>2.538552</td>
      <td>2.543976</td>
      <td>2.533802</td>
      <td>2.538541</td>
      <td>2.531302</td>
      <td>2.535771</td>
      <td>2.538775</td>
      <td>2.544050</td>
      <td>2.533076</td>
      <td>2.545031</td>
      <td>2.543609</td>
      <td>2.536865</td>
      <td>2.540552</td>
      <td>2.546328</td>
    </tr>
    <tr>
      <th>SIC</th>
      <td>2.577836</td>
      <td>2.584221</td>
      <td>2.599771</td>
      <td>2.615399</td>
      <td>2.584818</td>
      <td>2.599761</td>
      <td>2.602725</td>
      <td>2.617397</td>
      <td>2.599994</td>
      <td>2.615473</td>
      <td>2.614702</td>
      <td>2.636860</td>
      <td>2.615032</td>
      <td>2.618491</td>
      <td>2.632381</td>
      <td>2.648361</td>
    </tr>
  </tbody>
</table>
</div>



AICではARMA(2,3)、SICではARMA(1,1)が選択される
ARMA(2,3)モデルの残差について自己相関検定を行う（ただし、m=10とする）

結果、ARMA(2,3)モデルの残差が自己相関を持つとは言えない（適切なモデルであると言える）


```python
q_m_arma_2_3=acorr_ljungbox(arma_2_3.resid)[0]
p_value_arma_2_3=chi2.sf(q_m_arma_2_3,5)
lag=range(1,41)
test_arma_2_3=pd.DataFrame({'lag':lag,'autocorr':acf_y3[1:],'Q(m)':q_m_arma_2_3,'P_value':p_value_arma_2_3})
test_arma_2_3[['lag','autocorr','Q(m)','P_value']].iloc[0:10,:].T
```




<div>

<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>0</th>
      <th>1</th>
      <th>2</th>
      <th>3</th>
      <th>4</th>
      <th>5</th>
      <th>6</th>
      <th>7</th>
      <th>8</th>
      <th>9</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>lag</th>
      <td>1.000000</td>
      <td>2.000000</td>
      <td>3.000000</td>
      <td>4.000000</td>
      <td>5.000000</td>
      <td>6.000000</td>
      <td>7.000000</td>
      <td>8.000000</td>
      <td>9.000000</td>
      <td>10.000000</td>
    </tr>
    <tr>
      <th>autocorr</th>
      <td>0.867000</td>
      <td>0.655612</td>
      <td>0.519535</td>
      <td>0.419404</td>
      <td>0.327677</td>
      <td>0.223622</td>
      <td>0.136581</td>
      <td>0.079381</td>
      <td>0.017358</td>
      <td>-0.049272</td>
    </tr>
    <tr>
      <th>Q(m)</th>
      <td>0.012603</td>
      <td>0.163411</td>
      <td>0.333079</td>
      <td>0.397876</td>
      <td>2.769192</td>
      <td>2.770274</td>
      <td>4.432651</td>
      <td>5.903803</td>
      <td>5.942976</td>
      <td>6.232675</td>
    </tr>
    <tr>
      <th>P_value</th>
      <td>0.999999</td>
      <td>0.999458</td>
      <td>0.996974</td>
      <td>0.995388</td>
      <td>0.735517</td>
      <td>0.735351</td>
      <td>0.488945</td>
      <td>0.315692</td>
      <td>0.311807</td>
      <td>0.284232</td>
    </tr>
  </tbody>
</table>
</div>



ARMA(1,1)モデルの残差について自己相関検定を行う（ただし、m=10とする）
結果、ARMA(1,1)モデルの残差が自己相関を持つとは言えない（適切なモデルであると言える）


```python
q_m_arma_1_1=acorr_ljungbox(arma_1_1.resid)[0]
p_value_arma_1_1=chi2.sf(q_m_arma_1_1,8)
lag=range(1,41)
test_arma_1_1=pd.DataFrame({'lag':lag,'autocorr':acf_y3[1:],'Q(m)':q_m_arma_1_1,'P_value':p_value_arma_1_1})
test_arma_1_1[['lag','autocorr','Q(m)','P_value']].iloc[0:10,:].T
```




<div>

<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>0</th>
      <th>1</th>
      <th>2</th>
      <th>3</th>
      <th>4</th>
      <th>5</th>
      <th>6</th>
      <th>7</th>
      <th>8</th>
      <th>9</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>lag</th>
      <td>1.000000</td>
      <td>2.000000</td>
      <td>3.000000</td>
      <td>4.000000</td>
      <td>5.000000</td>
      <td>6.000000</td>
      <td>7.000000</td>
      <td>8.000000</td>
      <td>9.000000</td>
      <td>10.000000</td>
    </tr>
    <tr>
      <th>autocorr</th>
      <td>0.867000</td>
      <td>0.655612</td>
      <td>0.519535</td>
      <td>0.419404</td>
      <td>0.327677</td>
      <td>0.223622</td>
      <td>0.136581</td>
      <td>0.079381</td>
      <td>0.017358</td>
      <td>-0.049272</td>
    </tr>
    <tr>
      <th>Q(m)</th>
      <td>1.293730</td>
      <td>1.472988</td>
      <td>2.071876</td>
      <td>2.197028</td>
      <td>6.618453</td>
      <td>6.624761</td>
      <td>7.779293</td>
      <td>8.929929</td>
      <td>8.950020</td>
      <td>9.727831</td>
    </tr>
    <tr>
      <th>P_value</th>
      <td>0.995627</td>
      <td>0.993147</td>
      <td>0.978729</td>
      <td>0.974368</td>
      <td>0.578301</td>
      <td>0.577605</td>
      <td>0.455321</td>
      <td>0.348242</td>
      <td>0.346530</td>
      <td>0.284651</td>
    </tr>
  </tbody>
</table>
</div>

# 参考サイト
http://www.blackarbs.com/blog/time-series-analysis-in-python-linear-models-to-garch/11/1/2016#AR