自分の勉強用に、沖本竜義「経済・ファイナンスデータの計量経済分析」の章末問題をpythonで解いてみました。
jupyter notebook上で記述したものをほぼそのまま載せてあります。
ソースコードの細かい説明は省いてありますが、今後余裕があれば説明を加えようと思います。
こちらの記事は第1章の章末問題を解いたものになります。
第2章、第4章、第5章、第6章、第7章も公開しています。
各種設定・モジュールのインポート
%matplotlib inline
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import acf
from statsmodels.stats.diagnostic import acorr_ljungbox
import matplotlib as mpl
font = {"family":"IPAexGothic"}
mpl.rc('font', **font)
plt.rcParams["font.size"] = 12
1.3
(1.8)式のモデルからデータを生成する
gwn_0_1 = pd.Series(np.random.randn(100),name=u'\u03bc=0,\u03c3=1')
gwn_2_1 = pd.Series(2+np.random.randn(100),name=u'\u03bc=2,\u03c3=1')
gwn__2_1 = pd.Series(-2+np.random.randn(100),name=u'\u03bc=-2,\u03c3=1')
gwn_0_2 = pd.Series(2*np.random.randn(100),name=u'\u03bc=0,\u03c3=2')
gwn_0_3 = pd.Series(3*np.random.randn(100),name=u'\u03bc=0,\u03c3=3')
gwn_2_2 = pd.Series(2+2*np.random.randn(100),name=u'\u03bc=2,\u03c3=2')
プロットする
fig, ax = plt.subplots(3,2,figsize=(20,20))
ax[0,0].plot(gwn_0_1)
ax[0,0].set_title(u'\u03bc=0,\u03c3=1')
ax[0,1].plot(gwn_2_1)
ax[0,1].set_title(u'\u03bc=2,\u03c3=1')
ax[1,0].plot(gwn__2_1)
ax[1,0].set_title(u'\u03bc=-2,\u03c3=1')
ax[1,1].plot(gwn_0_2)
ax[1,1].set_title(u'\u03bc=0,\u03c3=2')
ax[2,0].plot(gwn_0_3)
ax[2,0].set_title(u'\u03bc=0,\u03c3=3')
ax[2,1].plot(gwn_2_2)
ax[2,1].set_title(u'\u03bc=2,\u03c3=2')
plt.show()
1.5
##(1)
ファイルをインポートする
eco = pd.read_csv('./input/economicdata.csv')
print(eco.shape)
eco.head()
(364, 7)
| date | topix | exrate | indprod | cpi | saunemp | intrate | |
|---|---|---|---|---|---|---|---|
| 0 | Jan-75 | 276.09 | 29.13 | 47.33 | 52.625 | 1.7 | 12.67 |
| 1 | Feb-75 | 299.81 | 29.70 | 46.86 | 52.723 | 1.8 | 13.00 |
| 2 | Mar-75 | 313.50 | 29.98 | 46.24 | 53.114 | 1.8 | 12.92 |
| 3 | Apr-75 | 320.57 | 29.80 | 47.33 | 54.092 | 1.8 | 12.02 |
| 4 | May-75 | 329.65 | 29.79 | 47.33 | 54.385 | 1.8 | 11.06 |
date列をdatetime型に変換してecoのDatetimeIndexとする
eco.index=pd.to_datetime(eco.date.values,format='%b-%y')
eco.drop('date',axis=1,inplace=True)
eco.head()
| topix | exrate | indprod | cpi | saunemp | intrate | |
|---|---|---|---|---|---|---|
| 1975-01-01 | 276.09 | 29.13 | 47.33 | 52.625 | 1.7 | 12.67 |
| 1975-02-01 | 299.81 | 29.70 | 46.86 | 52.723 | 1.8 | 13.00 |
| 1975-03-01 | 313.50 | 29.98 | 46.24 | 53.114 | 1.8 | 12.92 |
| 1975-04-01 | 320.57 | 29.80 | 47.33 | 54.092 | 1.8 | 12.02 |
| 1975-05-01 | 329.65 | 29.79 | 47.33 | 54.385 | 1.8 | 11.06 |
プロットする
title=['TOPIX','実効為替レート','鉱工業生産指数','CPI','失業率','コールレート']
ylim=[3500,120,120,120,6,14]
fig,ax = plt.subplots(nrows=3,ncols=2,figsize=[10,15])
for i in range(3):
for j in range(2):
ax[i,j].plot(eco.iloc[:,2*i+j])
ax[i,j].set_title(title[2*i+j])
ax[i,j].set_ylim(0,ylim[2*i+j])
ax[i,j].set_xlim(eco.index[0],eco.index[-1])
ax[i,j].hlines(ax[i,j].get_yticks(),ax[i,j].get_xlim()[0],ax[i,j].get_xlim()[1],linewidth=0.5)
ax[i,j].spines['right'].set_visible(False)
ax[i,j].set_xticks(eco.index[eco.index.is_year_start][::3])
ax[i,j].set_xticklabels(eco.index[eco.index.is_year_start][::3].strftime('%y'))
ax[i,j].set_xlabel('年')
plt.subplots_adjust(wspace=0.2,hspace=0.3)
plt.show()
(2)
対数差分系列を計算する
eco['topix_dlog'] = np.log(eco.topix).diff()
eco['exrate_dlog'] = np.log(eco.exrate).diff()
eco['indprod_dlog'] = np.log(eco.indprod).diff()
eco.head()
| topix | exrate | indprod | cpi | saunemp | intrate | topix_dlog | exrate_dlog | indprod_dlog | |
|---|---|---|---|---|---|---|---|---|---|
| 1975-01-01 | 276.09 | 29.13 | 47.33 | 52.625 | 1.7 | 12.67 | NaN | NaN | NaN |
| 1975-02-01 | 299.81 | 29.70 | 46.86 | 52.723 | 1.8 | 13.00 | 0.082422 | 0.019378 | -0.009980 |
| 1975-03-01 | 313.50 | 29.98 | 46.24 | 53.114 | 1.8 | 12.92 | 0.044650 | 0.009383 | -0.013319 |
| 1975-04-01 | 320.57 | 29.80 | 47.33 | 54.092 | 1.8 | 12.02 | 0.022301 | -0.006022 | 0.023299 |
| 1975-05-01 | 329.65 | 29.79 | 47.33 | 54.385 | 1.8 | 11.06 | 0.027931 | -0.000336 | 0.000000 |
(3)
プロットする
title=['TOPIX','実効為替レート','鉱工業生産指数']
ylim_u=[15,12,5]
ylim_l=[-20,-8,-5]
yticks_int=[5,2,1]
fig,ax = plt.subplots(nrows=3,ncols=1,figsize=[10,15])
for i in range(3):
ax[i].plot(eco.iloc[:,6+i]*100)
ax[i].set_title(title[i])
ax[i].set_ylim(ylim_l[i],ylim_u[i])
ax[i].set_xlim(eco.index[0],eco.index[-1])
ax[i].spines['right'].set_visible(False)
ax[i].set_xticks(eco.index[eco.index.is_year_start][::3])
ax[i].set_xticklabels(eco.index[eco.index.is_year_start][::3].strftime('%y'))
ax[i].set_yticks(range(ylim_l[i],ylim_u[i]+1)[::yticks_int[i]])
ax[i].hlines(ax[i].get_yticks(),ax[i].get_xlim()[0],ax[i].get_xlim()[1],linewidth=0.5)
ax[i].set_xlabel('年')
plt.subplots_adjust(wspace=0.2,hspace=0.3)
plt.show()
(4)
indprodの変化率の標本自己相関とQ(m)とP値を求める
autocorr_indprod=acf(eco.indprod_dlog[1:])[1:]
q_m_indprod=acorr_ljungbox(eco.indprod_dlog[1:])[0]
p_value_indprod=acorr_ljungbox(eco.indprod_dlog[1:])[1]
lag=range(1,41)
test_indprod=pd.DataFrame({'lag':lag,'autocorr':autocorr_indprod,'Q(m)':q_m_indprod,'P_value':p_value_indprod})
test_indprod[['lag','autocorr','Q(m)','P_value']].iloc[0:10,:].T
| 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | |
|---|---|---|---|---|---|---|---|---|---|---|
| lag | 1.000000e+00 | 2.000000e+00 | 3.000000e+00 | 4.000000e+00 | 5.000000e+00 | 6.000000e+00 | 7.000000e+00 | 8.000000e+00 | 9.000000e+00 | 1.000000e+01 |
| autocorr | -3.056815e-01 | 1.231154e-01 | 1.699805e-01 | 6.173203e-02 | 1.677640e-02 | 8.291975e-02 | -1.082355e-03 | 8.959420e-02 | 3.584451e-02 | -4.913541e-02 |
| Q(m) | 3.420025e+01 | 3.976335e+01 | 5.039731e+01 | 5.180377e+01 | 5.190793e+01 | 5.445974e+01 | 5.446017e+01 | 5.745610e+01 | 5.793698e+01 | 5.884316e+01 |
| P_value | 4.972263e-09 | 2.320056e-09 | 6.574768e-11 | 1.516132e-10 | 5.634492e-10 | 5.958688e-10 | 1.907042e-09 | 1.466964e-09 | 3.342017e-09 | 5.995413e-09 |
標本自己相関をプロットする
fig,ax=plt.subplots(1,1,figsize=[8,5])
ax.bar(test_indprod.lag[:20],test_indprod.autocorr[:20])
ax.set_ylim(-0.5,0.5)
ax.set_xticks([1,5,10,15,20])
ax.set_yticks(np.asarray(range(-5,6))*0.1)
ax.hlines([1.96/np.sqrt(eco.shape[0]-1),-1.96/np.sqrt(eco.shape[0]-1)],0,21,linestyles='dashed',linewidth=1)
ax.set_xlabel('ラグ')
ax.set_ylabel('自己相関')
ax.set_title('鉱工業生産指数の成長率')
plt.show()
(5)
TOPIXについて自己相関検定を行う
autocorr_topix=acf(eco.topix_dlog[1:])[1:]
q_m_topix=acorr_ljungbox(eco.topix_dlog[1:])[0]
p_value_topix=acorr_ljungbox(eco.topix_dlog[1:])[1]
lag=range(1,41)
test_topix=pd.DataFrame({'lag':lag,'autocorr':autocorr_topix,'Q(m)':q_m_topix,'P_value':p_value_topix})
test_topix[['lag','autocorr','Q(m)','P_value']].iloc[0:10,:].T
| 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | |
|---|---|---|---|---|---|---|---|---|---|---|
| lag | 1.000000e+00 | 2.000000e+00 | 3.000000e+00 | 4.000000e+00 | 5.000000e+00 | 6.000000e+00 | 7.000000 | 8.000000 | 9.000000 | 10.000000 |
| autocorr | 3.077999e-01 | 6.824984e-02 | 3.529433e-02 | 6.801588e-02 | 3.130508e-02 | -5.079917e-02 | 0.015837 | 0.025024 | 0.106257 | 0.064149 |
| Q(m) | 3.467592e+01 | 3.638552e+01 | 3.684399e+01 | 3.855135e+01 | 3.891405e+01 | 3.987178e+01 | 39.965130 | 40.198841 | 44.424647 | 45.969211 |
| P_value | 3.894192e-09 | 1.255984e-08 | 4.964713e-08 | 8.622952e-08 | 2.471338e-07 | 4.827198e-07 | 0.000001 | 0.000003 | 0.000001 | 0.000001 |
標本自己相関をプロットする
fig,ax=plt.subplots(1,1,figsize=[8,5])
ax.bar(test_topix.lag[:20],test_topix.autocorr[:20])
ax.set_ylim(-0.5,0.5)
ax.set_xticks([1,5,10,15,20])
ax.set_yticks(np.asarray(range(-5,6))*0.1)
ax.hlines([1.96/np.sqrt(eco.shape[0]-1),-1.96/np.sqrt(eco.shape[0]-1)],0,21,linestyles='dashed',linewidth=1)
ax.set_xlabel('ラグ')
ax.set_ylabel('自己相関')
ax.set_title('TOPIXの成長率')
plt.show()
実効為替レートについて自己相関検定を行う
autocorr_exrate=acf(eco.exrate_dlog[1:])[1:]
q_m_exrate=acorr_ljungbox(eco.exrate_dlog[1:])[0]
p_value_exrate=acorr_ljungbox(eco.exrate_dlog[1:])[1]
lag=range(1,41)
test_exrate=pd.DataFrame({'lag':lag,'autocorr':autocorr_exrate,'Q(m)':q_m_exrate,'P_value':p_value_exrate})
test_exrate[['lag','autocorr','Q(m)','P_value']].iloc[0:10,:].T
| 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | |
|---|---|---|---|---|---|---|---|---|---|---|
| lag | 1.000000e+00 | 2.000000e+00 | 3.000000e+00 | 4.000000e+00 | 5.000000e+00 | 6.000000e+00 | 7.000000e+00 | 8.000000e+00 | 9.000000e+00 | 1.000000e+01 |
| autocorr | 3.217763e-01 | 4.860159e-02 | 9.068437e-02 | 3.667455e-02 | -3.474773e-02 | -8.257040e-02 | -3.047480e-02 | 1.070266e-01 | 1.015502e-01 | -1.826474e-02 |
| Q(m) | 3.789649e+01 | 3.876344e+01 | 4.179009e+01 | 4.228649e+01 | 4.273335e+01 | 4.526370e+01 | 4.560934e+01 | 4.988452e+01 | 5.374427e+01 | 5.386948e+01 |
| P_value | 7.459915e-10 | 3.824946e-09 | 4.445457e-09 | 1.454942e-08 | 4.184633e-08 | 4.148156e-08 | 1.041334e-07 | 4.300973e-08 | 2.111351e-08 | 5.131560e-08 |
標本自己相関をプロットする
fig,ax=plt.subplots(1,1,figsize=[8,5])
ax.bar(test_exrate.lag[:20],test_exrate.autocorr[:20])
ax.set_ylim(-0.5,0.5)
ax.set_xticks([1,5,10,15,20])
ax.set_yticks(np.asarray(range(-5,6))*0.1)
ax.hlines([1.96/np.sqrt(eco.shape[0]-1),-1.96/np.sqrt(eco.shape[0]-1)],0,21,linestyles='dashed',linewidth=1)
ax.set_xlabel('ラグ')
ax.set_ylabel('自己相関')
ax.set_title('実効為替レートの成長率')
plt.show()
