In [ ]:
# 必要なライブラリをインストール
!pip install japanize_matplotlib
In [ ]:
# ライブラリを起動
from dateutil.relativedelta import relativedelta
import japanize_matplotlib
import math as math
import numpy as np
import pandas as pd
import scipy as sp
import random
import statsmodels.api as sm
from statsmodels.tsa.api import VAR
from sklearn.linear_model import LinearRegression
from matplotlib import pyplot as plt
import matplotlib.dates as mdates
import matplotlib.ticker as ticker
from matplotlib.ticker import MultipleLocator
from matplotlib.pylab import rcParams
rcParams["figure.figsize"] = 15, 6
# Colabでファイルを読み込むために、Google Driveをマウント
from google.colab import drive
drive.mount('/content/drive')
Mounted at /content/drive
(6) 構造ショックの部分識別による政策評価¶
インフレ率、輸入物価インフレ率、失業率、金利、M2、ドル円レートの6変量VARモデル¶
- 変数：インフレ率（除く生鮮食品・エネルギー、消費税率引き上げの影響は調整済み）、輸入物価インフレ率、失業率（季節調整済み）、金利（公定歩合とシャドーレートの接続系列）、M2（マネーストック統計以前の期間については、マネーサプライ統計におけるM2+CDの前年比伸び率を用いてM2の水準を補外）、ドル円レート
- 期間は1973年第1四半期から2023年第4四半期
- ラグ次数はBICにより1四半期が選択
インパルス応答関数のエラーバンドは$90\%$タイル
結果の解釈：
- 金融政策ショック→インフレ率：金融引き締めショックは、インフレ率を押し下げる
- 金融政策ショック→輸入物価インフレ率：金融引き締めショックは、円高などを通じて、輸入物価インフレ率が低下する
- 金融政策ショック→失業率：金融引き締めショックは、失業率を押し上げる
- 金融政策ショック→金利：金融引き締めショックは、金利が上昇する
- 金融政策ショック→M2：金融引き締めショックは、M2が減少する
- 金融政策ショック→ドル円レート：金融引き締めショックは、円高方向に作用する
In [ ]:
# データの読み込み
dinput = pd.read_excel('/content/drive/My Drive/Data/data_VAR6.xlsx',
sheet_name='data_q', header=[0])
data = np.zeros([len(dinput), 6])
for i in range(1, len(dinput)):
data[:, 0] = dinput.cpi / dinput.cpi.shift(4) * 100 - 100
data[:, 1] = dinput.ipi / dinput.ipi.shift(4) * 100 - 100
data[:, 2] = dinput.u
data[:, 3] = dinput.RSR
data[:, 4] = np.log(dinput.m2)
data[:, 5] = dinput.JPYUSD
data = data[8:]
dlen = len(data)
dti = pd.date_range('1973/01/01', periods=dlen, freq='QS')
fig, ax = plt.subplots(3, 2, figsize=[14, 16])
names = ['インフレ率', '輸入物価インフレ率', '失業率', '金利', 'M2',
'ドル円レート']
counter = 0
for i in range(3):
for j in range(2):
ax[i,j].plot(data[:, counter], linewidth=1, color='black')
ax[i,j].set_title('%s' %(names[counter]), fontsize=14)
ax[i,j].hlines([0], 0, dlen, color='black', linewidth=1)
ax[i,j].set_xlim(0, dlen)
axx = ax[i, j]
xaxis_ = axx.xaxis
new_xticks = [0, 20, 40, 60, 80, 100, 120, 140, 160, 180, 200]
xaxis_.set_major_locator(ticker.FixedLocator(new_xticks))
xaxis_.set_ticklabels(['1973', '1978', '1983', '1988', '1993', '1998',
'2003', '2008', '2013', '2018', '2023'])
ax[i,j].set_xlabel('年', fontsize=14, loc='center')
counter = counter + 1
ax[0,0].set_ylim(-5, 25)
ax[0,1].set_ylim(-50, 90)
ax[0,1].set_yticks(np.arange(-50,91,10))
ax[1,0].set_ylim(1, 6)
ax[1,1].set_ylim(-6, 10)
ax[2,0].set_ylim(13, 17)
ax[2,1].set_ylim(70, 320)
ax[2,1].set_yticks(np.arange(70,321,50))
In [ ]:
# VARモデルの推定
PERIOD = 50
CONFIDENCE_LEVEL = 0.9
lbound, ubound = sp.stats.norm.interval(confidence=CONFIDENCE_LEVEL, loc=0,
scale=1)
models = VAR(data)
results = models.fit(maxlags=12, ic='bic')
irf = results.irf(PERIOD)
impulse = irf.orth_irfs
lband = impulse + irf.stderr(orth=True) * lbound
uband = impulse + irf.stderr(orth=True) * ubound
# グラフの描画
names = ['インフレ率', '輸入物価インフレ率', '失業率', '金利', 'M2',
'ドル円レート']
fig, ax=plt.subplots(3, 2, figsize=[14,16])
counter = 0
for i in range(3):
for j in range(2):
ax[i, j].plot(impulse[:, counter, 3], linewidth=2, color='black')
ax[i, j].plot(lband[:, counter, 3], linewidth=1, linestyle='dashed', color='black')
ax[i, j].plot(uband[:, counter, 3], linewidth=1, linestyle='dashed', color='black')
ax[i, j].hlines([0], 0, PERIOD, color='black', linestyle='solid')
ax[i, j].set_xlim(0, PERIOD)
ax[i, j].set_xlabel('四半期', fontsize=14)
ax[i, j].set_ylabel('インパルス応答', fontsize=14)
ax[i, j].set_title('金融政策ショック→%s' %(names[counter]), fontsize=18)
counter = counter + 1
ax[0, 0].set_ylim(-0.4, 0.1)
ax[0, 1].set_ylim(-3, 2)
ax[1, 0].set_ylim(-0.04, 0.1)
ax[1, 1].set_ylim(-0.1, 0.5)
ax[2, 0].set_ylim(-0.015, 0.0025)
ax[2, 1].set_ylim(-4, 1)
plt.subplots_adjust(hspace=0.4)
plt.show()
In [ ]:
# VARモデルの誤差項の分散共分散行列
print('誤差項の分散共分散行列\n', results.sigma_u)
# VARモデルのショックのサイズの計算
L = np.linalg.cholesky(results.sigma_u)
Linv = np.linalg.inv(L)
print('ショックの1標準偏差 \n', np.diag(L))
print('コレスキー分解を行った行列C \n', L)
print('行列Cの逆行列A0 \n', Linv)
print('phi_pi \n', -1* Linv[3,0] * L[3,3])
print('phi_piim \n', -1* Linv[3,1] * L[3,3])
print('phi_u \n', -1* Linv[3,2] * L[3,3])
誤差項の分散共分散行列 [[ 6.16880427e-01 2.14040898e+00 -9.36363633e-03 8.81032908e-02 -5.79323437e-04 6.22864920e-01] [ 2.14040898e+00 7.43869849e+01 -2.25467215e-01 1.16034778e+00 -1.28654674e-03 3.01824889e+01] [-9.36363633e-03 -2.25467215e-01 1.54955012e-02 -5.07566795e-03 3.23180887e-05 -5.75252660e-02] [ 8.81032908e-02 1.16034778e+00 -5.07566795e-03 1.82409711e-01 4.09004258e-04 4.51708889e-01] [-5.79323437e-04 -1.28654674e-03 3.23180887e-05 4.09004258e-04 4.93812604e-05 2.34880016e-03] [ 6.22864920e-01 3.01824889e+01 -5.75252660e-02 4.51708889e-01 2.34880016e-03 5.83335852e+01]] ショックの1標準偏差 [7.85417359e-01 8.18292992e+00 1.21643796e-01 3.98507726e-01 6.87790772e-03 6.76523060e+00] コレスキー分解を行った行列C [[ 7.85417359e-01 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00] [ 2.72518675e+00 8.18292992e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00] [-1.19218607e-02 -2.35829855e-02 1.21643796e-01 0.00000000e+00 0.00000000e+00 0.00000000e+00] [ 1.12173852e-01 1.04443407e-01 -1.04835560e-02 3.98507726e-01 0.00000000e+00 0.00000000e+00] [-7.37599483e-04 8.84218248e-05 2.10530926e-04 1.21632691e-03 6.87790772e-03 0.00000000e+00] [ 7.93036865e-01 3.42436213e+00 2.68701674e-01 1.98635770e-02 3.70784966e-01 6.76523060e+00]] 行列Cの逆行列A0 [[ 1.27320843e+00 4.07393375e-17 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00] [-4.24020585e-01 1.22205617e-01 -0.00000000e+00 -0.00000000e+00 -0.00000000e+00 -0.00000000e+00] [ 4.25779393e-02 2.36919053e-02 8.22072337e+00 -7.74713511e-17 3.70024481e-14 0.00000000e+00] [-2.46138694e-01 -3.14051517e-02 2.16262842e-01 2.50936164e+00 -3.48103119e-16 0.00000000e+00] [ 1.84217606e-01 3.25759076e-03 -2.89879262e-01 -4.43769268e-01 1.45393053e+02 -3.31346424e-16] [ 5.43133189e-02 -6.28842407e-02 -3.11258423e-01 1.69540525e-02 -7.96862093e+00 1.47814621e-01]] phi_pi 0.09808817119755878 phi_piim 0.01251519561129239 phi_u -0.08618241350493404
インフレ率、輸入物価インフレ率、失業率、金利、M2、ドル円レートの6変量の局所予測法¶
- 変数：インフレ率、輸入物価インフレ率、失業率、金利、M2、ドル円レートの6変数
- 期間は1973年第1四半期から2023年第4四半期まで
- ラグ次数は4四半期とする
- 信頼区間の幅は90%ile
In [ ]:
#データの読み込み
dinput = pd.read_excel('/content/drive/My Drive/Data/data_VAR6.xlsx',
sheet_name='data_q', header=[0])
data = np.zeros([len(dinput), 6])
for i in range(1, len(dinput)):
data[:, 0] = dinput.cpi / dinput.cpi.shift(4) * 100 - 100
data[:, 1] = dinput.ipi / dinput.ipi.shift(4) * 100 - 100
data[:, 2] = dinput.u
data[:, 3] = dinput.RSR
data[:, 4] = np.log(dinput.m2)
data[:, 5] = dinput.JPYUSD
data = data[8:]
dlen = len(data)
dti = pd.date_range('1973/01/01', periods=dlen, freq='QS')
# LPの推定
L = 4 + 1
MAX_H = 21
CONFIDENCE_LEVEL = 0.90
coef_result, std_result = np.zeros([MAX_H, len(data.T)]), np.zeros(
[MAX_H, len(data.T)])
lbound, ubound = sp.stats.norm.interval(confidence=CONFIDENCE_LEVEL, loc=0,
scale=1)
ssize = len(data)
counter = 0
h = 0
for dep in range(0, len(data.T), 1):
for indep in range(3, 4, 1):
if dep == 0 and indep == 3:
coef_result[h, dep] = np.array(0)
elif dep == 1 and indep == 3:
coef_result[h, dep] = np.array(0)
elif dep == 2 and indep == 3:
coef_result[h, dep] = np.array(0)
elif dep == 3 and indep == 3:
coef_result[h, dep] = np.array(1)
elif dep == 4 and indep == 3:
sample_y = pd.DataFrame(data[L + h - 1 : ssize, dep])
sample_x = pd.DataFrame(data[L - 1 : L - 1 + len(sample_y), indep])
for l in range(1, L, 1):
sample_x_lagged = pd.DataFrame(data[L - 1 - l :
L - 1 - l + len(sample_y)])
sample_x = pd.concat([sample_x, sample_x_lagged], axis=1)
sample_x = sm.add_constant(sample_x)
model = sm.OLS(sample_y, sample_x)
result = model.fit(cov_type='HAC', cov_kwds={'maxlags':1,
'kernel':'bartlett'})
coef_result[h, dep] = np.array(result.params.iloc[indep + 1])
std_result[h, dep] = np.array(result.bse.iloc[indep + 1])
elif dep == 5 and indep == 3:
sample_y = pd.DataFrame(data[L + h - 1 : ssize, dep])
sample_x = pd.DataFrame(data[L - 1 : L - 1 + len(sample_y), indep])
for l in range(1, L, 1):
sample_x_lagged = pd.DataFrame(data[L - 1 - l :
L - 1 - l + len(sample_y)])
sample_x = pd.concat([sample_x, sample_x_lagged], axis=1)
sample_x = sm.add_constant(sample_x)
model = sm.OLS(sample_y, sample_x)
result = model.fit(cov_type='HAC', cov_kwds={'maxlags':1,
'kernel':'bartlett'})
coef_result[h, dep] = np.array(result.params.iloc[indep + 1])
std_result[h, dep] = np.array(result.bse.iloc[indep + 1])
for dep in range(0, len(data.T), 1):
for indep in range(3, 4, 1):
for h in range(1, MAX_H, 1):
if dep == 0 and indep == 3:
sample_y = pd.DataFrame(data[L + h - 1 : ssize, dep])
sample_x = pd.DataFrame(data[L - 1 : L - 1 + len(sample_y)])
for l in range(1, L, 1):
sample_x_lagged = pd.DataFrame(data[L - 1 - l :
L - 1 - l +
len(sample_y)])
sample_x = pd.concat([sample_x, sample_x_lagged], axis=1)
sample_x = sm.add_constant(sample_x)
model = sm.OLS(sample_y, sample_x)
result = model.fit(cov_type='HAC',
cov_kwds={'maxlags':h + 1,
'kernel':'bartlett'})
coef_result[h, counter] = np.array(result.params.iloc[indep+1])
std_result[h, counter] = np.array(result.bse.iloc[indep+1])
if dep == 1 and indep == 3:
sample_y = pd.DataFrame(data[L + h - 1 : ssize, dep])
sample_x = pd.DataFrame(data[L - 1 : L - 1 + len(sample_y)])
for l in range(1, L, 1):
sample_x_lagged = pd.DataFrame(data[L - 1 - l :
L - 1 - l +
len(sample_y)])
sample_x = pd.concat([sample_x, sample_x_lagged], axis=1)
sample_x = sm.add_constant(sample_x)
model = sm.OLS(sample_y, sample_x)
result = model.fit(cov_type='HAC',
cov_kwds={'maxlags':h + 1,
'kernel':'bartlett'})
coef_result[h, counter] = np.array(result.params.iloc[indep+1])
std_result[h, counter] = np.array(result.bse.iloc[indep+1])
if dep == 2 and indep == 3:
sample_y = pd.DataFrame(data[L + h - 1 : ssize, dep])
sample_x = pd.DataFrame(data[L - 1 : L - 1 + len(sample_y)])
for l in range(1, L, 1):
sample_x_lagged = pd.DataFrame(data[L - 1 - l :
L - 1 - l +
len(sample_y)])
sample_x = pd.concat([sample_x, sample_x_lagged], axis=1)
sample_x = sm.add_constant(sample_x)
model = sm.OLS(sample_y, sample_x)
result = model.fit(cov_type='HAC',
cov_kwds={'maxlags':h + 1,
'kernel':'bartlett'})
coef_result[h, counter] = np.array(result.params.iloc[indep+1])
std_result[h, counter] = np.array(result.bse.iloc[indep+1])
if dep == 3 and indep == 3:
sample_y = pd.DataFrame(data[L + h - 1 : ssize, dep])
sample_x = pd.DataFrame(data[L - 1 : L - 1 + len(sample_y)])
for l in range(1, L, 1):
sample_x_lagged = pd.DataFrame(data[L - 1 - l :
L - 1 - l +
len(sample_y)])
sample_x = pd.concat([sample_x, sample_x_lagged], axis=1)
sample_x = sm.add_constant(sample_x)
model = sm.OLS(sample_y, sample_x)
result = model.fit(cov_type='HAC',
cov_kwds={'maxlags':h + 1,
'kernel':'bartlett'})
coef_result[h, counter] = np.array(result.params.iloc[indep+1])
std_result[h, counter] = np.array(result.bse.iloc[indep+1])
if dep == 4 and indep == 3:
sample_y = pd.DataFrame(data[L + h - 1 : ssize, dep])
sample_x = pd.DataFrame(data[L - 1 : L - 1 + len(sample_y)])
for l in range(1, L, 1):
sample_x_lagged = pd.DataFrame(data[L - 1 - l :
L - 1 - l +
len(sample_y)])
sample_x = pd.concat([sample_x, sample_x_lagged], axis=1)
sample_x = sm.add_constant(sample_x)
model = sm.OLS(sample_y, sample_x)
result = model.fit(cov_type='HAC',
cov_kwds={'maxlags':h + 1,
'kernel':'bartlett'})
coef_result[h, counter] = np.array(result.params.iloc[indep+1])
std_result[h, counter] = np.array(result.bse.iloc[indep+1])
if dep == 5 and indep == 3:
sample_y = pd.DataFrame(data[L + h - 1 : ssize, dep])
sample_x = pd.DataFrame(data[L - 1 : L - 1 + len(sample_y)])
for l in range(1, L, 1):
sample_x_lagged = pd.DataFrame(data[L - 1 - l :
L - 1 - l +
len(sample_y)])
sample_x = pd.concat([sample_x, sample_x_lagged], axis=1)
sample_x = sm.add_constant(sample_x)
model = sm.OLS(sample_y, sample_x)
result = model.fit(cov_type='HAC',
cov_kwds={'maxlags':h + 1,
'kernel':'bartlett'})
coef_result[h, counter] = np.array(result.params.iloc[indep+1])
std_result[h, counter] = np.array(result.bse.iloc[indep+1])
counter = counter + 1
# VAR(1)モデルのインパルス応答関数
models = VAR(data)
results = models.fit(1)
irf = results.irf(MAX_H)
impulse = irf.orth_irfs
# グラフの描画
vnames = ['インフレ率', '輸入物価インフレ率', '失業率', '金利', 'M2',
'ドル円レート']
snames = ['金融政策ショック']
fig, ax=plt.subplots(6, 1, figsize=[6,30])
gcounter = 0
for i in range(6):
ax[i].plot(coef_result[:, gcounter], linewidth=2, color='black',
label='LPのインパルス応答関数')
ax[i].plot(coef_result[:, gcounter] + std_result[:, gcounter] * lbound,
linewidth=1, linestyle='dashed', color='black',
label='信頼区間の上限下限')
ax[i].plot(coef_result[:, gcounter] + std_result[:, gcounter] * ubound,
linewidth=1, linestyle='dashed', color='black')
ax[i].hlines([0], 0, MAX_H, color='black', linestyle='solid')
ax[i].set_xlim(0, MAX_H - 1)
ax[i].xaxis.set_major_locator(MultipleLocator(4))
ax[i].set_xlabel('四半期', fontsize=14)
ax[i].set_ylabel('インパルス応答', fontsize=14)
ax[i].set_title('%s→%s' %(snames[0],vnames[i]), fontsize=14)
ax[0].legend(loc='lower left', frameon=False, fontsize=11)
ax[0].set_ylim(-1.5, 0.5)
ax[1].set_ylim(-20, 15)
ax[2].set_ylim(-0.2, 0.6)
ax[3].set_ylim(-1.0, 1.5)
ax[4].set_ylim(-0.06, 0.01)
ax[5].set_ylim(-20, 5)
gcounter = gcounter + 1
plt.subplots_adjust(wspace=0.4, hspace=0.3)
plt.show()
