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

自分の勉強用に、沖本竜義「経済・ファイナンスデータの計量経済分析」の章末問題を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()