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

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

こちらの記事は第4章の章末問題を解いたものになります。
[第1章](https://qiita.com/mckeeeen/items/666e14c94e1548130646#15-5)、[第2章](https://qiita.com/mckeeeen/items/666e14c94e1548130646)、[第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 statsmodels.tsa.vector_ar.var_model import VAR

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




# 4.5

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


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

    (1391, 8)





<div>

<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>Date</th>
      <th>ca</th>
      <th>fr</th>
      <th>ge</th>
      <th>it</th>
      <th>jp</th>
      <th>uk</th>
      <th>us</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>2003/1/1</td>
      <td>560.099</td>
      <td>902.081</td>
      <td>724.932</td>
      <td>290.187</td>
      <td>1593.175</td>
      <td>791.076</td>
      <td>824.583</td>
    </tr>
    <tr>
      <th>1</th>
      <td>2003/1/2</td>
      <td>574.701</td>
      <td>927.206</td>
      <td>768.150</td>
      <td>296.963</td>
      <td>1578.214</td>
      <td>797.813</td>
      <td>852.219</td>
    </tr>
    <tr>
      <th>2</th>
      <td>2003/1/3</td>
      <td>580.212</td>
      <td>929.297</td>
      <td>768.411</td>
      <td>298.757</td>
      <td>1578.411</td>
      <td>800.175</td>
      <td>851.935</td>
    </tr>
    <tr>
      <th>3</th>
      <td>2003/1/6</td>
      <td>589.619</td>
      <td>943.002</td>
      <td>788.164</td>
      <td>303.273</td>
      <td>1619.700</td>
      <td>803.966</td>
      <td>871.515</td>
    </tr>
    <tr>
      <th>4</th>
      <td>2003/1/7</td>
      <td>585.822</td>
      <td>923.785</td>
      <td>774.054</td>
      <td>297.892</td>
      <td>1590.951</td>
      <td>793.625</td>
      <td>865.992</td>
    </tr>
  </tbody>
</table>
</div>



Date列をdatetime型に変換してmsciのDatetimeIndexとする


```python
msci.index=pd.to_datetime(msci.Date.values)
msci.drop('Date',axis=1,inplace=True)
msci.head()
```




<div>

<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>ca</th>
      <th>fr</th>
      <th>ge</th>
      <th>it</th>
      <th>jp</th>
      <th>uk</th>
      <th>us</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>2003-01-01</th>
      <td>560.099</td>
      <td>902.081</td>
      <td>724.932</td>
      <td>290.187</td>
      <td>1593.175</td>
      <td>791.076</td>
      <td>824.583</td>
    </tr>
    <tr>
      <th>2003-01-02</th>
      <td>574.701</td>
      <td>927.206</td>
      <td>768.150</td>
      <td>296.963</td>
      <td>1578.214</td>
      <td>797.813</td>
      <td>852.219</td>
    </tr>
    <tr>
      <th>2003-01-03</th>
      <td>580.212</td>
      <td>929.297</td>
      <td>768.411</td>
      <td>298.757</td>
      <td>1578.411</td>
      <td>800.175</td>
      <td>851.935</td>
    </tr>
    <tr>
      <th>2003-01-06</th>
      <td>589.619</td>
      <td>943.002</td>
      <td>788.164</td>
      <td>303.273</td>
      <td>1619.700</td>
      <td>803.966</td>
      <td>871.515</td>
    </tr>
    <tr>
      <th>2003-01-07</th>
      <td>585.822</td>
      <td>923.785</td>
      <td>774.054</td>
      <td>297.892</td>
      <td>1590.951</td>
      <td>793.625</td>
      <td>865.992</td>
    </tr>
  </tbody>
</table>
</div>



欠損データ(e.g. 2005-01-01)があるので、グラフの目盛用にその年で1番最初の日かどうかをカラムに加える


```python
msci['is_year_start']=[msci.index[i].year!=msci.index[i-1].year for i in range(msci.shape[0])]
```

JP,UK,USのMSCIをプロットする


```python
title=['JP','UK','US']
fig,ax=plt.subplots(1,3,figsize=[20,5])
for i in range(3):
    ax[i].plot(msci.iloc[:,i+4])
    ax[i].set_title(title[i])
    ax[i].set_xticks(msci.index[msci.is_year_start])
    ax[i].set_xlabel('year',fontsize=12)
    ax[i].set_ylabel('MSCI',fontsize=12)
plt.show()
```


![output_10_0.png](https://qiita-image-store.s3.amazonaws.com/0/235450/06c9f041-cd66-95af-63b5-a4735cf327da.png)



対数差分系列を計算して株式収益率（PER）のDataFrameを作る


```python
per=pd.DataFrame(np.log(msci.iloc[:,:7]).diff(),index=msci.index[1:])
per.head()
```




<div>

<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>ca</th>
      <th>fr</th>
      <th>ge</th>
      <th>it</th>
      <th>jp</th>
      <th>uk</th>
      <th>us</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>2003-01-02</th>
      <td>0.025736</td>
      <td>0.027471</td>
      <td>0.057907</td>
      <td>0.023082</td>
      <td>-0.009435</td>
      <td>0.008480</td>
      <td>0.032966</td>
    </tr>
    <tr>
      <th>2003-01-03</th>
      <td>0.009544</td>
      <td>0.002253</td>
      <td>0.000340</td>
      <td>0.006023</td>
      <td>0.000125</td>
      <td>0.002956</td>
      <td>-0.000333</td>
    </tr>
    <tr>
      <th>2003-01-06</th>
      <td>0.016083</td>
      <td>0.014640</td>
      <td>0.025381</td>
      <td>0.015003</td>
      <td>0.025822</td>
      <td>0.004727</td>
      <td>0.022723</td>
    </tr>
    <tr>
      <th>2003-01-07</th>
      <td>-0.006461</td>
      <td>-0.020589</td>
      <td>-0.018065</td>
      <td>-0.017902</td>
      <td>-0.017909</td>
      <td>-0.012946</td>
      <td>-0.006357</td>
    </tr>
    <tr>
      <th>2003-01-08</th>
      <td>-0.016002</td>
      <td>-0.025601</td>
      <td>-0.041735</td>
      <td>-0.013862</td>
      <td>-0.016723</td>
      <td>-0.012756</td>
      <td>-0.014646</td>
    </tr>
  </tbody>
</table>
</div>



欠損データ(e.g. 2005-01-01)があるので、グラフの目盛用にその年で1番最初の日かどうかをカラムに加える


```python
per['is_year_start']=[per.index[i].year!=per.index[i-1].year for i in range(per.shape[0])]
```

JP,UK,USのPERをプロットする


```python
names=['JP','UK','US']
fig,ax=plt.subplots(1,3,figsize=[20,5])
for i in range(3):
    ax[i].plot(per.iloc[:,i+4])
    ax[i].set_title(names[i])
    ax[i].set_xticks(per.index[per.is_year_start])
    ax[i].set_xlabel('year',fontsize=12)
    ax[i].set_ylabel('PER',fontsize=12)
plt.show()
```


![output_16_0.png](https://qiita-image-store.s3.amazonaws.com/0/235450/bf6c02bc-1988-6a3b-5b20-b20bd92f7c41.png)



0期から10期までのVARモデルを推定する


```python
var=[]
for i in range(11):
     var.append(VAR(per[['jp','uk','us']]).fit(maxlags=i))
```

各次数のモデルのAICを確認


```python
aic=[]
for i in range(11):
    aic.append(var[i].aic)
pd.Series(aic,index=range(11),name='AIC')
```




    0    -27.559452
    1    -27.967305
    2    -27.995273
    3    -28.001670
    4    -28.003014
    5    -28.000596
    6    -28.002719
    7    -27.993201
    8    -27.986282
    9    -27.980579
    10   -27.981139
    Name: AIC, dtype: float64



教科書ではVAR(3)がAIC最小だったが、今回はVAR(4)がAIC最小となった
（これは、フィッティングに使うmethodが異なるからだと考えていますが、もし理由が分かる方いればご教授ください...）
今回は教科書に合わせてVAR(3)モデルを使う

グレンジャー因果性検定を行う


```python
names=['JP','UK','US']
stat=[]
pval=[]
null=[]
for i in range(3):
    for j in range(3):
        if i!=j:
            test=var[3].test_causality(i,j,kind='f',signif=0.05,verbose=False)
            stat.append(test['statistic'].round(1))
            pval.append(test['pvalue'].round(3))
            null.append('%s->%s'%(names[j],names[i]))
        else:
            pass
pd.DataFrame({'Statistics':stat,'P_value':pval},index=null,columns=['Statistics','P_value']).T
```




<div>

<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>UK-&gt;JP</th>
      <th>US-&gt;JP</th>
      <th>JP-&gt;UK</th>
      <th>US-&gt;UK</th>
      <th>JP-&gt;US</th>
      <th>UK-&gt;US</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>Statistics</th>
      <td>12.0</td>
      <td>53.5</td>
      <td>0.800</td>
      <td>71.9</td>
      <td>0.200</td>
      <td>0.600</td>
    </tr>
    <tr>
      <th>P_value</th>
      <td>0.0</td>
      <td>0.0</td>
      <td>0.519</td>
      <td>0.0</td>
      <td>0.882</td>
      <td>0.617</td>
    </tr>
  </tbody>
</table>
</div>



インパルス応答関数（IRF）を計算する
信頼区間はブートストラップ法によるものを採用


```python
irf=var[3].irf(10).orth_irfs
interval=var[3].irf_errband_mc(orth=True,T=10,signif=0.05)
```

プロットする


```python
names=['JP','UK','US']
fig,ax=plt.subplots(3,3,figsize=[15,15])
for i in range(3):
    for j in range(3):
        ax[i,j].plot(irf[:,i,j]*100,linewidth=2)
        ax[i,j].plot(interval[0][:,i,j]*100,linestyle='dashed',color='green',linewidth=1)
        ax[i,j].plot(interval[1][:,i,j]*100,linestyle='dashed',color='green',linewidth=1)
        ax[i,j].hlines(0,0,10)
        ax[i,j].set_xlim(0,10)
        ax[i,j].set_ylim(-0.2,1.2)
        ax[i,j].set_xlabel('lag',fontsize=12)
        ax[i,j].set_ylabel('IRF',fontsize=12)
        ax[i,j].set_title('%s->%s'%(names[j],names[i]))
plt.subplots_adjust(wspace=0.3, hspace=0.3)
plt.show()
```


![output_27_0.png](https://qiita-image-store.s3.amazonaws.com/0/235450/1c1641a1-f153-6883-263b-0cf1f1a3904a.png)



分散分解分析を行うため、相対的分散寄与率（RVC）を計算する


```python
decomp=var[3].fevd(11).decomp
```

プロットする


```python
names=['JP','UK','US']
fig,ax=plt.subplots(3,3,figsize=[15,15])
for i in range(3):
    for j in range(3):
        ax[i,j].plot(decomp[i,:,j]*100,linewidth=2)
        ax[i,j].set_xlim(0,10)
        ax[i,j].set_ylim(0,100)
        ax[i,j].set_xlabel('lag',fontsize=12)
        ax[i,j].set_ylabel('RVC',fontsize=12)
        ax[i,j].set_title('%s->%s'%(names[j],names[i]))
plt.subplots_adjust(wspace=0.3, hspace=0.3)
plt.show()
```


![output_31_0.png](https://qiita-image-store.s3.amazonaws.com/0/235450/1a0af6e2-b07a-2811-c3f6-b274447f5946.png)


# 4.6

## (1)

4.5で計算済みなので、ここでは表示だけする
ただし、PERは%表示に直していない


```python
per[['jp','fr','ca']].head()
```




<div>

<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>jp</th>
      <th>fr</th>
      <th>ca</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>2003-01-02</th>
      <td>-0.009435</td>
      <td>0.027471</td>
      <td>0.025736</td>
    </tr>
    <tr>
      <th>2003-01-03</th>
      <td>0.000125</td>
      <td>0.002253</td>
      <td>0.009544</td>
    </tr>
    <tr>
      <th>2003-01-06</th>
      <td>0.025822</td>
      <td>0.014640</td>
      <td>0.016083</td>
    </tr>
    <tr>
      <th>2003-01-07</th>
      <td>-0.017909</td>
      <td>-0.020589</td>
      <td>-0.006461</td>
    </tr>
    <tr>
      <th>2003-01-08</th>
      <td>-0.016723</td>
      <td>-0.025601</td>
      <td>-0.016002</td>
    </tr>
  </tbody>
</table>
</div>



## (2)

日本、フランス、カナダの順に外生性が高いような再帰的構造が仮定されている
各国の証券取引所の取引時間を調べてみると、フランスはイギリスと、カナダはアメリカと同じで、UTC基準では以下のようになる

日本 0:00-6:00
フランス 8:00-16:30
カナダ 14:30-21:00

取引時間を見れば分かる通り、日本、フランス、カナダの順で取引開始時間が早く、しかもフランスとカナダの2時間以外、2つ以上の取引所が同時に開いている時間帯がない
PERのデータが日次データであるため、同時点で日本のPERは他国の影響を受けることはなく、フランスのPERは日本の影響は受けるがカナダの影響は小さく、カナダのPERは日本とフランスの影響を受けると考えられる
つまり、日本、フランス、カナダの順にが高いような再帰的構造が仮定されていることは妥当だと言える

## (3)

時価総額ベースでみても、月間売買代金ベースでみても、取引所の規模は日本、フランス（ユーロネクスト）、カナダの順で大きい
よって、週次や年次のデータであっても日本、フランス、カナダの順に外生性が高いと考えられる
しかし、日次データと違って週次や年次のデータでは再帰的構造の仮定を支持する根拠はないため、実際は他の順に並べて試して結果を比較して妥当な順番を見つけ出すことになると考えられる
（筆者は計量経済学の専門ではなく、教科書の知識を基に考察していますので、ここはあまり参考にしないでください）

## (4)

JP,FR,CAのPERをプロットする


```python
names=['JP','FR','CA']
fig,ax=plt.subplots(1,3,figsize=[20,5])
for i in range(3):
    ax[i].plot(per.loc[:,names[i].lower()])
    ax[i].set_title(names[i])
    ax[i].set_xticks(per.index[per.is_year_start])
    ax[i].set_xlabel('year',fontsize=12)
    ax[i].set_ylabel('PER',fontsize=12)
plt.show()
```


![output_42_0.png](https://qiita-image-store.s3.amazonaws.com/0/235450/0116d62a-c8a4-1195-ec73-c2d3f1418577.png)



0期から10期までのAICを計算する


```python
var=[]
for i in range(11):
     var.append(VAR(per[['jp','fr','ca']]).fit(maxlags=i))
```

各次数のモデルのAICを確認


```python
aic=[]
for i in range(11):
    aic.append(var[i].aic)
pd.Series(aic,index=range(11),name='AIC')
```




    0    -27.177659
    1    -27.385471
    2    -27.392535
    3    -27.390516
    4    -27.383707
    5    -27.386572
    6    -27.383015
    7    -27.378343
    8    -27.369757
    9    -27.364752
    10   -27.359088
    Name: AIC, dtype: float64



AICが最小となったのはVAR(2)モデル

## (5)

グレンジャー因果性検定を行う


```python
names=['JP','FR','CA']
stat=[]
pval=[]
null=[]
for i in range(3):
    for j in range(3):
        if i!=j:
            test=var[2].test_causality(i,j,kind='f',signif=0.05,verbose=False)
            stat.append(test['statistic'].round(1))
            pval.append(test['pvalue'].round(3))
            null.append('%s->%s'%(names[j],names[i]))
        else:
            pass
pd.DataFrame({'Statistics':stat,'P_value':pval},index=null,columns=['Statistics','P_value']).T
```




<div>

<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>FR-&gt;JP</th>
      <th>CA-&gt;JP</th>
      <th>JP-&gt;FR</th>
      <th>CA-&gt;FR</th>
      <th>JP-&gt;CA</th>
      <th>FR-&gt;CA</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>Statistics</th>
      <td>39.7</td>
      <td>24.7</td>
      <td>4.100</td>
      <td>20.9</td>
      <td>2.300</td>
      <td>2.700</td>
    </tr>
    <tr>
      <th>P_value</th>
      <td>0.0</td>
      <td>0.0</td>
      <td>0.017</td>
      <td>0.0</td>
      <td>0.104</td>
      <td>0.069</td>
    </tr>
  </tbody>
</table>
</div>



フランスに対する日本とカナダのグレンシャー因果検定では、有意水準5%で帰無仮説が棄却され、グレンジャー因果性が存在することが示唆された
つまり、フランスのPERを予測する上では、日本とカナダのPERは有用であるといえる
ただし、因果の方向性や大きさは分からない

## (6)

インパルス応答関数（IRF）を計算する
信頼区間はブートストラップ法によるものを採用


```python
irf=var[2].irf(10).orth_irfs
interval=var[2].irf_errband_mc(orth=True,T=10,signif=0.05)
```

プロットする


```python
names=['JP','FR','CA']
fig,ax=plt.subplots(3,3,figsize=[15,15])
for i in range(3):
    for j in range(3):
        ax[i,j].plot(irf[:,i,j]*100,linewidth=2)
        ax[i,j].plot(interval[0][:,i,j]*100,linestyle='dashed',color='green',linewidth=1)
        ax[i,j].plot(interval[1][:,i,j]*100,linestyle='dashed',color='green',linewidth=1)
        ax[i,j].hlines(0,0,10)
        ax[i,j].set_xlim(0,10)
        ax[i,j].set_ylim(-0.2,1.2)
        ax[i,j].set_xlabel('lag',fontsize=12)
        ax[i,j].set_ylabel('IRF',fontsize=12)
        ax[i,j].set_title('%s->%s'%(names[j],names[i]))
plt.subplots_adjust(wspace=0.3, hspace=0.3)
plt.show()
```


![output_56_0.png](https://qiita-image-store.s3.amazonaws.com/0/235450/d5e84bc1-56d9-13a9-14d2-6c20e086e558.png)



フランスの株式市場に対するショックが日本の株式市場に与える影響（FR->JP）をみると、IRFは0期では0%程度、1期では0.45%程度、2期以降は0という結果になった
また、FR->FRのIRFは0期では1.1%程度である
つまり、フランスのPERが1.1%上昇したら、1日後の日本のPERは0.45%上昇し、2日後以降には影響を与えないということが示唆される

## (7)

分散分解分析を行うため、相対的分散寄与率（RVC）を計算する


```python
decomp=var[2].fevd(11).decomp
```

プロットする


```python
names=['JP','FR','CA']
fig,ax=plt.subplots(3,3,figsize=[15,15])
for i in range(3):
    for j in range(3):
        ax[i,j].plot(decomp[i,:,j]*100,linewidth=2)
        ax[i,j].set_xlim(0,10)
        ax[i,j].set_ylim(0,100)
        ax[i,j].set_xlabel('lag',fontsize=12)
        ax[i,j].set_ylabel('RVC',fontsize=12)
        ax[i,j].set_title('%s->%s'%(names[j],names[i]))
plt.subplots_adjust(wspace=0.3, hspace=0.3)
plt.show()
```


![output_62_0.png](https://qiita-image-store.s3.amazonaws.com/0/235450/9f0b47b2-936d-18a7-2a57-95fc2dd57791.png)



カナダの株式市場の変動（JP,FR,CA->CA）をみると、日本の株式市場の変動の寄与が3%程度、フランスの株式市場の変動の寄与が27%程度という結果になった
つまり、カナダの株式市場の変動に対して、日本の株式市場の変動のは3%の説明力しか持たないが、フランスの株式市場の変動は27%と比較的大きな説明力を持つことが示唆される


# 参考サイト

statsmodelsのドキュメントを読んだだけでは目的の変数がどこに格納されているのか分からなかったので、ソースコードをかなり参照しました
: https://github.com/statsmodels
取引所データ
: https://ja.wikipedia.org/wiki/%E8%A8%BC%E5%88%B8%E5%8F%96%E5%BC%95%E6%89%80