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.ticker as ticker
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
(7) 分散分解と歴史分解¶
(7-1) 平均2乗予測誤差とグレンジャー因果性¶
(7-1-1) 物価と賃金間のグレンジャー因果性¶
- 変数：インフレ率（消費者物価指数、除く生鮮食品・エネルギー、消費税率の影響を調整済み）、名目賃金上昇率（毎月勤労統計における現金給与総額の前年比）の2変数
- 期間は1972年1月から2023年12月
- ラグ次数はBICにより13か月が選択
- 「賃金から物価」と「物価から賃金」のいずれにおいても、「グレンジャー因果性はない」との帰無仮説が棄却され、グレンジャーの意味での因果性が確認された。
In [ ]:
# データの読み込み
dinput = pd.read_excel('/content/drive/My Drive/Data/data_VAR7.xlsx',
sheet_name='data_m', header=[0])
data = np.zeros([len(dinput), 2])
for i in range(1, len(dinput)):
data[:, 0] = dinput.infl
data[:, 1] = dinput.wage
data = data[12:]
dlen = len(data)
dti = np.arange(0, dlen, 1)
# グラフの図示
plt.figure(figsize=(10, 6))
plt.plot(dti, data[:, 0], label='インフレ率', color='black')
plt.plot(dti, data[:, 1], label='賃金上昇率', color='grey')
plt.hlines([0], 0, dlen, color='black', linewidth=0.75)
plt.legend(frameon=False, ncol=2, fontsize=14)
plt.tick_params(labelsize=14)
plt.xlim(0, dlen)
plt.ylim(-20, 40)
plt.xlabel('年', fontsize=14, loc='center')
plt.ylabel('％', fontsize=14)
axx = plt.gca()
xaxis_ = axx.xaxis
new_xticks = [0, 60, 120, 180, 240, 300, 360, 420, 480, 540, 600]
xaxis_.set_major_locator(ticker.FixedLocator(new_xticks))
xaxis_.set_ticklabels(['1972', '1977', '1982', '1987', '1992', '1997', '2002',
'2007', '2012', '2017', '2022'])
plt.xlabel('年', fontsize=14, loc='center')
plt.ylabel('％', fontsize=14, loc='center')
Out[ ]:
Text(0, 0.5, '％')
In [ ]:
# VARを推定
model = VAR(data)
result = model.fit(maxlags=24, ic="bic")
print(f"ラグ次数: {result.k_ar}")
# グレンジャーの因果性検定（ワルド統計量で検定、有意水準5％）
# 帰無仮説「y2（賃金上昇率）からy1（インフレ率）へグレンジャー因果性がない」
test12_outcome = result.test_causality("y1", "y2", kind="wald", signif=0.05)
print(test12_outcome.summary())
# 帰無仮説「y1（インフレ率）からy2（賃金上昇率）へグレンジャー因果性がない」
test21_outcome = result.test_causality("y2", "y1", kind="wald", signif=0.05)
print(test21_outcome.summary())
ラグ次数: 13
Granger causality Wald-test. H_0: y2 does not Granger-cause y1. Conclusion: reject H_0 at 5% significance level.
========================================
Test statistic Critical value p-value df
----------------------------------------
85.71 22.36 0.000 13
----------------------------------------
Granger causality Wald-test. H_0: y1 does not Granger-cause y2. Conclusion: reject H_0 at 5% significance level.
========================================
Test statistic Critical value p-value df
----------------------------------------
126.3 22.36 0.000 13
----------------------------------------
(7-2) 分散分解¶
(7-2-1) 3変量VARモデルの分散分解¶
- 変数：インフレ率、失業率、金利の3変数（各系列の平均が0となるように調整）
- 期間は1972年第1四半期から2023年第4四半期
- ラグ次数はBICにより2四半期が選択
In [ ]:
# データの読み込み
dinput = pd.read_excel('/content/drive/My Drive/Data/data_VAR7.xlsx',
sheet_name='data_q', header=[0])
data = np.zeros([len(dinput), 3])
for i in range(1, len(dinput)):
data[:, 0] = dinput.cpi / dinput.cpi.shift(4) * 100 - 100
data[:, 1] = dinput.u
data[:, 2] = dinput.RSR
data = data[4:]
data[:, 0] = data[:, 0] - data[:, 0].mean()
data[:, 1] = data[:, 1] - data[:, 1].mean()
data[:, 2] = data[:, 2] - data[:, 2].mean()
dlen = len(data)
# VARモデルの推定
model = VAR(data)
PERIOD = 100
results = model.fit(ic='bic', maxlags=8)
lag = results.k_ar
fit = results.fittedvalues
resid = results.resid
sigma = np.dot((resid.T), resid) / ((len(data) - lag) - (len(data.T) * lag))
inva = np.linalg.cholesky(sigma)
for mm in range(0, len(data.T), 1):
fit_raw = results.params.T
fit_comp_before = np.append(np.eye(len(data.T) * (lag - 1)),
np.zeros([len(data.T) * (lag - 1),
len(data.T)]), axis=1)
fit_comp = np.append(fit_raw[:, 1:(len(data.T) * lag + 1)],
fit_comp_before, axis=0)
In [ ]:
# 分散分解
psi = np.zeros([PERIOD, len(data.T), len(data.T)])
sigma_fevd = np.eye(len(data.T))
inva_fevd = np.linalg.cholesky(sigma_fevd)
impulse_fevd = np.ones([len(data.T), 1])
response_fevd = np.zeros(shape=[len(data.T), len(data.T), PERIOD])
for mm in range(0, len(data.T), 1):
response_temp_fevd = inva_fevd * impulse_fevd
response_fevd[mm, :, 0] = response_temp_fevd[:, mm]
impulse_big_fevd = (np.append(response_fevd[mm, :, 0].T, np.zeros(
[1, len(data.T) * (lag - 1)])))
fit_raw_fevd = results.params.T
fit_comp_before_fevd = np.append(np.eye(len(data.T) * (lag - 1)),
np.zeros([len(data.T) * (lag - 1),
len(data.T)]), axis=1)
fit_comp_fevd = np.append(fit_raw_fevd[:, 1:(len(data.T) * lag + 1)],
fit_comp_before_fevd, axis=0)
fit_comp_eye_fevd = np.eye(len(fit_comp_fevd))
for recur in range(1, PERIOD, 1):
fit_comp_eye_fevd = np.dot(fit_comp_fevd, fit_comp_eye_fevd)
response_big_fevd = np.dot(fit_comp_eye_fevd, impulse_big_fevd)
response_fevd[mm, :, recur] = response_big_fevd[0:len(data.T)]
for i in range(0, len(data.T), 1):
for j in range(0, PERIOD, 1):
psi[j, :, i] = response_fevd[i, :, j]
mse = np.zeros([PERIOD, len(data.T), len(data.T)])
mse_j = np.zeros([PERIOD, len(data.T), len(data.T)])
fevd = np.zeros([len(data.T), PERIOD, len(data.T)])
se = np.zeros([PERIOD, len(data.T)])
for i in range(0, len(data.T), 1):
mse[0, :, :] = sigma
for j in range(1, PERIOD, 1):
mse[j, :, :] = mse[j - 1, :, :] + \
np.dot(np.dot(psi[j, :, :], sigma), psi[j, :, :].T)
col = pd.DataFrame(inva[:, i])
mse_j[0, :, :] = np.dot(col, col.T)
for j in range(1, PERIOD, 1):
mse_j[j, :, :] = mse_j[j - 1, :, :] + \
np.dot(np.dot(psi[j, :, :], np.dot(col, col.T)), psi[j, :, :].T)
fecd = mse_j / mse
for j in range(0, PERIOD, 1):
for k in range(0, len(data.T), 1):
fevd[k, j, i] = fecd[j, k, k]
se[j, :] = np.sqrt(np.diag(mse[j, :, :]).T)
# 分散分解の結果一覧
color_dic = {"black":"\033[30m", "red":"\033[31m", "green":"\033[32m",
"yellow":"\033[33m", "blue":"\033[34m", "end":"\033[0m"}
def print_color(text, color="red"):
print(color_dic[color] + text + color_dic["end"])
vnames = ['インフレ率', '失業率', '金利']
for i in range(0, len(data.T), 1):
print_color(text='%s' %(vnames[i]), color='red')
for j in (0, 1, 3, 7, 99):
print_color('%s 期' %(j+1), 'blue')
print('供 給 ショック:', "{:.1f}".format(fevd[i,j,0] * 100))
print('需 要 ショック:', "{:.1f}".format(fevd[i,j,1] * 100))
print('金融政策ショック:', "{:.1f}".format(fevd[i,j,2] * 100))
インフレ率 1 期 供 給 ショック: 100.0 需 要 ショック: 0.0 金融政策ショック: 0.0 2 期 供 給 ショック: 98.4 需 要 ショック: 0.8 金融政策ショック: 0.8 4 期 供 給 ショック: 93.6 需 要 ショック: 3.4 金融政策ショック: 3.0 8 期 供 給 ショック: 86.4 需 要 ショック: 8.5 金融政策ショック: 5.0 100 期 供 給 ショック: 55.3 需 要 ショック: 35.4 金融政策ショック: 9.3 失業率 1 期 供 給 ショック: 0.9 需 要 ショック: 99.1 金融政策ショック: 0.0 2 期 供 給 ショック: 0.8 需 要 ショック: 99.1 金融政策ショック: 0.2 4 期 供 給 ショック: 0.6 需 要 ショック: 99.0 金融政策ショック: 0.4 8 期 供 給 ショック: 0.3 需 要 ショック: 99.4 金融政策ショック: 0.2 100 期 供 給 ショック: 3.0 需 要 ショック: 82.4 金融政策ショック: 14.6 金利 1 期 供 給 ショック: 5.5 需 要 ショック: 0.4 金融政策ショック: 94.1 2 期 供 給 ショック: 7.8 需 要 ショック: 1.4 金融政策ショック: 90.8 4 期 供 給 ショック: 10.1 需 要 ショック: 3.6 金融政策ショック: 86.3 8 期 供 給 ショック: 10.8 需 要 ショック: 7.6 金融政策ショック: 81.5 100 期 供 給 ショック: 6.4 需 要 ショック: 43.2 金融政策ショック: 50.4
(7-2-2) 6変量VARモデルの分散分解¶
- 変数：インフレ率、輸入物価インフレ率、失業率、金利、M2、ドル円レートの6変数（各系列の平均が0となるように調整）
- 期間は1973年第1四半期から2023年第4四半期
- ラグ次数はBICにより1四半期が選択
- 構造ショックは部分識別（金融政策ショックのみ）のため、ショックは「金融政策ショック」と「その他のショック」となっている
In [ ]:
# データの読み込み
dinput = pd.read_excel('/content/drive/My Drive/Data/data_VAR7.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:]
data[:, 0] = data[:, 0] - data[:, 0].mean()
data[:, 1] = data[:, 1] - data[:, 1].mean()
data[:, 2] = data[:, 2] - data[:, 2].mean()
data[:, 3] = data[:, 3] - data[:, 3].mean()
data[:, 4] = data[:, 4] - data[:, 4].mean()
data[:, 5] = data[:, 5] - data[:, 5].mean()
dlen = len(data)
# VARモデルの推定
model = VAR(data)
PERIOD = 100
results = model.fit(ic='bic', maxlags=12)
lag = results.k_ar
fit = results.fittedvalues
resid = results.resid
sigma = np.dot((resid.T), resid) / ((len(data) - lag) - (len(data.T) * lag))
inva = np.linalg.cholesky(sigma)
for mm in range(0, len(data.T), 1):
fit_raw = results.params.T
fit_comp_before = np.append(np.eye(len(data.T) * (lag - 1)),
np.zeros([len(data.T) * (lag - 1),
len(data.T)]), axis=1)
fit_comp = np.append(fit_raw[:, 1:(len(data.T) * lag + 1)],
fit_comp_before, axis=0)
In [ ]:
# 分散分解
psi = np.zeros([PERIOD, len(data.T), len(data.T)])
sigma_fevd = np.eye(len(data.T))
inva_fevd = np.linalg.cholesky(sigma_fevd)
impulse_fevd = np.ones([len(data.T), 1])
response_fevd = np.zeros(shape=[len(data.T), len(data.T), PERIOD])
for mm in range(0, len(data.T), 1):
response_temp_fevd = inva_fevd * impulse_fevd
response_fevd[mm, :, 0] = response_temp_fevd[:, mm]
impulse_big_fevd = (np.append(response_fevd[mm, :, 0].T, np.zeros(
[1, len(data.T) * (lag - 1)])))
fit_raw_fevd = results.params.T
fit_comp_before_fevd = np.append(np.eye(len(data.T) * (lag - 1)),
np.zeros([len(data.T) * (lag - 1),
len(data.T)]), axis=1)
fit_comp_fevd = np.append(fit_raw_fevd[:, 1:(len(data.T) * lag + 1)],
fit_comp_before_fevd, axis=0)
fit_comp_eye_fevd = np.eye(len(fit_comp_fevd))
for recur in range(1, PERIOD, 1):
fit_comp_eye_fevd = np.dot(fit_comp_fevd, fit_comp_eye_fevd)
response_big_fevd = np.dot(fit_comp_eye_fevd, impulse_big_fevd)
response_fevd[mm, :, recur] = response_big_fevd[0:len(data.T)]
for i in range(0, len(data.T), 1):
for j in range(0, PERIOD, 1):
psi[j, :, i] = response_fevd[i, :, j]
mse = np.zeros([PERIOD, len(data.T), len(data.T)])
mse_j = np.zeros([PERIOD, len(data.T), len(data.T)])
fevd = np.zeros([len(data.T), PERIOD, len(data.T)])
se = np.zeros([PERIOD, len(data.T)])
for i in range(0, len(data.T), 1):
mse[0, :, :] = sigma
for j in range(1, PERIOD, 1):
mse[j, :, :] = mse[j - 1, :, :] + \
np.dot(np.dot(psi[j, :, :], sigma), psi[j, :, :].T)
col = pd.DataFrame(inva[:, i])
mse_j[0, :, :] = np.dot(col, col.T)
for j in range(1, PERIOD, 1):
mse_j[j, :, :] = mse_j[j - 1, :, :] + \
np.dot(np.dot(psi[j, :, :], np.dot(col, col.T)), psi[j, :, :].T)
fecd = mse_j / mse
for j in range(0, PERIOD, 1):
for k in range(0, len(data.T), 1):
fevd[k, j, i] = fecd[j, k, k]
se[j, :] = np.sqrt(np.diag(mse[j, :, :]).T)
# 分散分解の結果一覧
color_dic = {"black":"\033[30m", "red":"\033[31m", "green":"\033[32m",
"yellow":"\033[33m", "blue":"\033[34m", "end":"\033[0m"}
def print_color(text, color="red"):
print(color_dic[color] + text + color_dic["end"])
vnames = ['インフレ率', '輸入物価インフレ率', '失業率', '金利', 'M2', 'ドル円']
for i in range(0, len(data.T), 1):
print_color('%s' %(vnames[i]), 'red')
for j in (0, 1, 3, 7, 99):
print_color('%s 期' %(j+1), 'blue')
print('金融政策ショック:', "{:.1f}".format(fevd[i,j,3] * 100))
print('その他のショック:', "{:.1f}".format(100 - fevd[i,j,3] * 100))
インフレ率 1 期 金融政策ショック: 0.0 その他のショック: 100.0 2 期 金融政策ショック: 0.1 その他のショック: 99.9 4 期 金融政策ショック: 1.0 その他のショック: 99.0 8 期 金融政策ショック: 5.3 その他のショック: 94.7 100 期 金融政策ショック: 15.6 その他のショック: 84.4 輸入物価インフレ率 1 期 金融政策ショック: 0.0 その他のショック: 100.0 2 期 金融政策ショック: 0.3 その他のショック: 99.7 4 期 金融政策ショック: 1.9 その他のショック: 98.1 8 期 金融政策ショック: 6.8 その他のショック: 93.2 100 期 金融政策ショック: 9.8 その他のショック: 90.2 失業率 1 期 金融政策ショック: 0.0 その他のショック: 100.0 2 期 金融政策ショック: 0.1 その他のショック: 99.9 4 期 金融政策ショック: 0.5 その他のショック: 99.5 8 期 金融政策ショック: 2.0 その他のショック: 98.0 100 期 金融政策ショック: 11.9 その他のショック: 88.1 金利 1 期 金融政策ショック: 87.1 その他のショック: 12.9 2 期 金融政策ショック: 79.4 その他のショック: 20.6 4 期 金融政策ショック: 63.3 その他のショック: 36.7 8 期 金融政策ショック: 42.5 その他のショック: 57.5 100 期 金融政策ショック: 29.5 その他のショック: 70.5 M2 1 期 金融政策ショック: 3.0 その他のショック: 97.0 2 期 金融政策ショック: 1.9 その他のショック: 98.1 4 期 金融政策ショック: 1.0 その他のショック: 99.0 8 期 金融政策ショック: 2.2 その他のショック: 97.8 100 期 金融政策ショック: 26.8 その他のショック: 73.2 ドル円 1 期 金融政策ショック: 0.0 その他のショック: 100.0 2 期 金融政策ショック: 0.0 その他のショック: 100.0 4 期 金融政策ショック: 0.2 その他のショック: 99.8 8 期 金融政策ショック: 1.5 その他のショック: 98.5 100 期 金融政策ショック: 14.9 その他のショック: 85.1
(7-3) 歴史分解¶
(7-3-1) GDPギャップのAR(1)モデルの分解¶
GDPギャップ（日本銀行が推計している需給ギャップ）のAR(1)モデルを推定し、その歴史分解を行う。
In [ ]:
# データの読み込み
dinput = pd.read_excel('/content/drive/My Drive/Data/data_VAR7.xlsx',
sheet_name='data_q', header=[0])
data = np.zeros([len(dinput), 1])
for i in range(1, len(dinput)):
data[:, 0] = dinput.gapboj
data = pd.DataFrame(data[48:])
dlen = len(data)
dep_y = data[1:]
ind_x = data[:-1]
model = LinearRegression()
model.fit(ind_x, dep_y)
y_hat = model.predict(ind_x)
resid = np.array(dep_y - y_hat)
# ショックの加重和
shocks = np.zeros(len(dep_y))
shocks[0] = resid[0]
for i in range(1, len(dep_y), 1):
shocks[i] = resid[i] + shocks[i -1] * model.coef_
# 平均
con = np.zeros(len(dep_y))
conv = model.intercept_ / (1 - model.coef_)
for i in range(0, len(dep_y), 1):
con[i] = conv
# 初期値の影響（残差として計算）
init = np.squeeze(dep_y.to_numpy()) - con - shocks
# 歴史分解を図示
plt.figure(figsize=(10, 6))
dti1 = np.arange(1, dlen, 1)
dti2 = np.arange(0, dlen, 1)
temp = [con[i] if (shocks[i] > 0) else init[i] for i in range(0, len(dep_y), 1)]
plt.bar(dti1, con, label="平均（μ）", width=1, color="black", edgecolor='black')
plt.bar(dti1, init, label="初期値の影響", width=1, color="grey", edgecolor=
'black')
plt.bar(dti1, shocks, bottom=temp, label="ショックの加重和", width=1,
color="white", edgecolor='black')
plt.plot(dti2, data, color='black', marker='o', markerfacecolor='white',
label='GDPギャップ')
plt.plot(dti1, shocks + con, color='lightgrey', linewidth=2,
label='GDPギャップ（初期値の影響を除く）')
plt.hlines([0], -1, dlen, color='black', linewidth=0.75)
plt.legend(frameon=False, ncol = 2, fontsize=12)
plt.tick_params(labelsize=14)
plt.xlim(-1, dlen)
plt.ylim(-6, 6)
plt.xlabel('年', fontsize=14, loc='center')
plt.ylabel('％', fontsize=14)
axx = plt.gca()
xaxis_ = axx.xaxis
new_xticks = [0, 20, 40, 60, 80, 100, 120, 140, 160]
xaxis_.set_major_locator(ticker.FixedLocator(new_xticks))
xaxis_.set_ticklabels(['1983', '1988', '1993', '1998', '2003', '2008', '2013',
'2018', '2023'])
plt.xlabel('年', fontsize=14, loc='center')
plt.ylabel('％', fontsize=14, loc='center')
Out[ ]:
Text(0, 0.5, '％')
(7-3-2) 3変数VARモデルの歴史分解¶
- 変数：インフレ率、失業率、金利の3変数（各系列の平均が0となるように調整）
- 期間は1972年第1四半期から2023年第4四半期
- ラグ次数はBICにより2四半期が選択
In [ ]:
# データの読み込み
dinput = pd.read_excel('/content/drive/My Drive/Data/data_VAR7.xlsx',
sheet_name='data_q', header=[0])
data = np.zeros([len(dinput), 3])
for i in range(1, len(dinput)):
data[:, 0] = dinput.cpi / dinput.cpi.shift(4) * 100 - 100
data[:, 1] = dinput.u
data[:, 2] = dinput.RSR
data = data[4:]
data[:, 0] = data[:, 0] - data[:, 0].mean()
data[:, 1] = data[:, 1] - data[:, 1].mean()
data[:, 2] = data[:, 2] - data[:, 2].mean()
dlen = len(data)
# VARモデルの推定
model = VAR(data)
PERIOD = 25
results = model.fit(ic='bic', maxlags=8)
lag = results.k_ar
fit = results.fittedvalues
resid = results.resid
sigma = np.dot((resid.T), resid) / ((len(data) - lag) - (len(data.T) * lag))
inva = np.linalg.cholesky(sigma)
for mm in range(0, len(data.T), 1):
fit_raw = results.params.T
fit_comp_before = np.append(np.eye(len(data.T) * (lag - 1)),
np.zeros([len(data.T) * (lag - 1),
len(data.T)]), axis=1)
fit_comp = np.append(fit_raw[:, 1:(len(data.T) * lag + 1)],
fit_comp_before, axis=0)
In [ ]:
# 歴史分解
inva_big = np.zeros([len(data.T) * lag, len(data.T)])
inva_big[0:len(data.T), :] = inva
icomp = np.append(np.eye(len(data.T)), np.zeros(
[len(data.T), (lag - 1) * len(data.T)]), axis=1)
hdshock_big = np.zeros([len(data.T), len(data.T) * lag, len(data) - lag + 1])
hdshock = np.zeros([len(data.T), len(data.T), len(data) - lag + 1])
for i in range(0, len(data.T), 1):
eps_big = np.zeros([len(data.T), len(data) - lag + 1])
eps_big[i, 1:] = (np.dot(np.linalg.inv(inva), resid.T))[i, :]
for j in range(1, len(data) - lag + 1, 1):
hdshock_big[i, :, j] = np.dot(inva_big, eps_big[:, j]) + np.dot(
fit_comp, hdshock_big[i, :, j - 1])
hdshock[i, :, j] = np.dot(icomp, hdshock_big[i, :, j])
shock = np.zeros([len(data.T), len(data), len(data.T)])
for i in range(0, len(data.T), 1):
for j in range(0, len(data.T), 1):
shock[i, :, j] = np.append(np.zeros(lag), hdshock[j, i, 1:], axis=0)
x = np.zeros([len(data) - lag, len(data.T)])
for i in range(0, lag, 1):
x = np.append(data[i:len(data) - lag + i, :], x, axis=1)
x = x[:, 0:len(x.T) - len(data.T)]
hdinit_big = np.zeros([lag * len(data.T), len(data) - lag + 1])
hdinit = np.zeros([len(data.T), len(data) - lag + 1])
hdinit_big[:, 0] = x[0, :].T
hdinit[:, 0] = np.dot(icomp, hdinit_big[:, 0])
for i in range(1, len(data) - lag + 1, 1):
hdinit_big[:, i] = np.dot(fit_comp, hdinit_big[:, i - 1])
hdinit[:, i] = np.dot(icomp, hdinit_big[:, i])
hdconst_big = np.zeros([lag * len(data.T), len(data) - lag + 1])
hdconst = np.zeros([len(data.T), len(data) - lag + 1])
cc = np.zeros(len(data.T) * lag)
cc[0:len(data.T)] = fit_raw[:, 0]
for i in range(1, len(data) - lag + 1):
hdconst_big[:, i] = cc + np.dot(fit_comp, hdconst_big[:, i - 1])
hdconst[:, i] = np.dot(icomp, hdconst_big[:, i])
hdendo = np.sum(hdshock, axis=0) + hdinit + hdconst
# 歴史分解を図示
dti = np.arange(dlen)
vnames = ['インフレ率', '失業率', '金利']
location = ['upper right', 'lower right', 'upper right']
fig, ax = plt.subplots(3, 1, figsize=[10, 20])
initvalue = np.zeros([len(data), len(data.T)])
zerovalue = np.zeros([len(data), len(data.T)])
for vl in range(0, len(data.T), 1):
initvalue[:, vl] = data[:, vl] - shock[vl, :, 0] - shock[
vl, :, 1] - shock[vl, :, 2]
for vl in range(0, len(data.T), 1):
ax[vl].set_title('%s' %(vnames[vl]), fontsize=18)
const = zerovalue[:, vl]
temp1 = [const[i] if (const[i] > 0 and shock[vl, i, 0] > 0) \
else const[i] if (const[i] <= 0 and shock[vl, i, 0] <= 0) \
else 0 for i in range(0, len(shock[vl, :, 0]))]
ax[vl].bar(dti, shock[vl, :, 0], bottom=temp1, label='供給ショック',
width=1, color='lightgrey', edgecolor='lightgrey')
temp2 = [const[i] + shock[vl, i, 0] \
if (const[i] > 0 and shock[vl, i, 0] > 0 and shock[vl, i, 1] > 0) \
else const[i] \
if (const[i] > 0 and shock[vl, i, 0] <= 0 and shock[vl, i, 1] > 0) \
else shock[vl, i, 0] \
if (const[i] <= 0 and shock[vl, i, 0] > 0 and shock[vl, i, 1] > 0) \
else 0 \
if (const[i] <= 0 and shock[vl, i, 0] <= 0 and shock[vl, i, 1] > 0) \
else const[i] + shock[vl, i, 0] \
if (const[i] <= 0 and shock[vl, i, 0] <= 0 and shock[vl, i, 1] <= 0) \
else const[i] \
if (const[i] <= 0 and shock[vl, i, 0] > 0 and shock[vl, i, 1] <= 0) \
else shock[vl, i, 0] \
if (const[i] > 0 and shock[vl, i, 0] <= 0 and shock[vl, i, 1] <= 0) \
else 0 for i in range(0, len(shock[vl, :, 0]))]
ax[vl].bar(dti, shock[vl, :, 1], bottom=temp2, label='需要ショック',
width=1, color='grey', edgecolor='grey')
temp3 = [const[i] + shock[vl, i, 0] + shock[vl, i, 1] \
if (const[i] > 0 and shock[vl, i, 0] > 0 and shock[vl, i, 1] > 0
and shock[vl, i, 2] > 0) \
else const[i] + shock[vl, i, 0] \
if (const[i] > 0 and shock[vl, i, 0] > 0 and shock[vl, i, 1] <= 0
and shock[vl, i, 2] > 0) \
else const[i] + shock[vl, i, 1] \
if (const[i] > 0 and shock[vl, i, 0] <= 0 and shock[vl, i, 1] > 0
and shock[vl, i, 2] > 0) \
else shock[vl, i, 0] + shock[vl, i, 1] \
if (const[i] <= 0 and shock[vl, i, 0] > 0 and shock[vl, i, 1] > 0
and shock[vl, i, 2] > 0) \
else const[i] \
if (const[i] > 0 and shock[vl, i, 0] <= 0 and shock[vl, i, 1] <= 0
and shock[vl, i, 2] > 0) \
else shock[vl, i, 0] \
if (const[i] <= 0 and shock[vl, i, 0] > 0 and shock[vl, i, 1] <= 0
and shock[vl, i, 2] > 0) \
else shock[vl, i, 1] \
if (const[i] <= 0 and shock[vl, i, 0] <= 0 and shock[vl, i, 1] > 0
and shock[vl, i, 2] > 0) \
else 0 \
if (const[i] <= 0 and shock[vl, i, 0] <= 0 and shock[vl, i, 1] <= 0
and shock[vl, i, 2] > 0) \
else const[i] + shock[vl, i, 0] + shock[vl, i, 1] \
if (const[i] <= 0 and shock[vl, i, 0] <= 0 and shock[vl, i, 1] <= 0
and shock[vl, i, 2] <= 0) \
else const[i] + shock[vl, i, 0] \
if (const[i] <= 0 and shock[vl, i, 0] <= 0 and shock[vl, i, 1] > 0
and shock[vl, i, 2] <= 0) \
else const[i] + shock[vl, i, 1] \
if (const[i] <= 0 and shock[vl, i, 0] > 0 and shock[vl, i, 1] <= 0
and shock[vl, i, 2] <= 0) \
else shock[vl, i, 0] + shock[vl, i, 1] \
if (const[i] > 0 and shock[vl, i, 0] <= 0 and shock[vl, i, 1] <= 0
and shock[vl, i, 2] <= 0) \
else const[i] \
if (const[i] <= 0 and shock[vl, i, 0] > 0 and shock[vl, i, 1] > 0
and shock[vl, i, 2] <= 0) \
else shock[vl, i, 0] \
if (const[i] > 0 and shock[vl, i, 0] <= 0 and shock[vl, i, 1] > 0
and shock[vl, i, 2] <= 0) \
else shock[vl, i, 1] \
if (const[i] > 0 and shock[vl, i, 0] > 0 and shock[vl, i, 1] <= 0
and shock[vl, i, 2] <= 0) \
else 0 for i in range(0, len(shock[vl, :, 0]))]
ax[vl].bar(dti, shock[vl, :, 2], bottom=temp3, label='金融政策ショック',
width=1, color='black', edgecolor='black')
ax[vl].hlines([0], 0, dlen, color='black', linewidth=1)
ax[vl].legend(loc='%s' %(location[vl]), frameon=False, ncol = 2,
fontsize=14)
ax[vl].tick_params(labelsize=14)
ax[vl].set_xlim(0, dlen)
ax[vl].set_xlabel('年', fontsize=14, loc='center')
ax[vl].set_ylabel('％', fontsize=14)
axx = ax[vl]
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(['1972', '1977', '1982', '1987', '1992', '1997',
'2002', '2007', '2012', '2017', '2022'])
ax[0].set_ylim(-10, 20)
ax[1].set_ylim(-3, 3)
ax[2].set_ylim(-8, 10)
(7-3-3) 6変量VARモデルの歴史分解¶
- 変数：インフレ率、輸入物価インフレ率、失業率、金利、M2、ドル円レートの6変数（各系列の平均が0となるように調整）
- 期間は1973年第1四半期から2023年第4四半期
- ラグ次数はBICにより1四半期が選択
- 構造ショックは部分識別（金融政策ショックのみ）のため、ショックは「金融政策ショック」と「その他のショック」となっている
In [ ]:
# データの読み込み
dinput = pd.read_excel('/content/drive/My Drive/Data/data_VAR7.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:]
data[:, 0] = data[:, 0] - data[:, 0].mean()
data[:, 1] = data[:, 1] - data[:, 1].mean()
data[:, 2] = data[:, 2] - data[:, 2].mean()
data[:, 3] = data[:, 3] - data[:, 3].mean()
data[:, 4] = data[:, 4] - data[:, 4].mean()
data[:, 5] = data[:, 5] - data[:, 5].mean()
dlen = len(data)
# VARモデルの推定
model = VAR(data)
PERIOD = 25
results = model.fit(ic='bic', maxlags=12)
lag = results.k_ar
fit = results.fittedvalues
resid = results.resid
sigma = np.dot((resid.T), resid) / ((len(data) - lag) - (len(data.T) * lag))
inva = np.linalg.cholesky(sigma)
for mm in range(0, len(data.T), 1):
fit_raw = results.params.T
fit_comp_before = np.append(np.eye(len(data.T) * (lag - 1)),
np.zeros([len(data.T) * (lag - 1),
len(data.T)]), axis=1)
fit_comp = np.append(fit_raw[:, 1:(len(data.T) * lag + 1)],
fit_comp_before, axis=0)
In [ ]:
# 歴史分解
inva_big = np.zeros([len(data.T) * lag, len(data.T)])
inva_big[0:len(data.T), :] = inva
icomp = np.append(np.eye(len(data.T)), np.zeros(
[len(data.T), (lag - 1) * len(data.T)]), axis=1)
hdshock_big = np.zeros([len(data.T), len(data.T) * lag, len(data) - lag + 1])
hdshock = np.zeros([len(data.T), len(data.T), len(data) - lag + 1])
for i in range(0, len(data.T), 1):
eps_big = np.zeros([len(data.T), len(data) - lag + 1])
eps_big[i, 1:] = (np.dot(np.linalg.inv(inva), resid.T))[i, :]
for j in range(1, len(data) - lag + 1, 1):
hdshock_big[i, :, j] = np.dot(inva_big, eps_big[:, j]) + np.dot(
fit_comp, hdshock_big[i, :, j - 1])
hdshock[i, :, j] = np.dot(icomp, hdshock_big[i, :, j])
shock = np.zeros([len(data.T), len(data), len(data.T)])
for i in range(0, len(data.T), 1):
for j in range(0, len(data.T), 1):
shock[i, :, j] = np.append(np.zeros(lag), hdshock[j, i, 1:], axis=0)
x = np.zeros([len(data) - lag, len(data.T)])
for i in range(0, lag, 1):
x = np.append(data[i:len(data) - lag + i, :], x, axis=1)
x = x[:, 0:len(x.T) - len(data.T)]
hdinit_big = np.zeros([lag * len(data.T), len(data) - lag + 1])
hdinit = np.zeros([len(data.T), len(data) - lag + 1])
hdinit_big[:, 0] = x[0, :].T
hdinit[:, 0] = np.dot(icomp, hdinit_big[:, 0])
for i in range(1, len(data) - lag + 1, 1):
hdinit_big[:, i] = np.dot(fit_comp, hdinit_big[:, i - 1])
hdinit[:, i] = np.dot(icomp, hdinit_big[:, i])
hdconst_big = np.zeros([lag * len(data.T), len(data) - lag + 1])
hdconst = np.zeros([len(data.T), len(data) - lag + 1])
cc = np.zeros(len(data.T) * lag)
cc[0:len(data.T)] = fit_raw[:, 0]
for i in range(1, len(data) - lag + 1):
hdconst_big[:, i] = cc + np.dot(fit_comp, hdconst_big[:, i - 1])
hdconst[:, i] = np.dot(icomp, hdconst_big[:, i])
hdendo = np.sum(hdshock, axis=0) + hdinit + hdconst
# 歴史分解を図示
dti = np.arange(dlen)
vnames = ['インフレ率', '輸入物価インフレ率', '失業率', '金利', 'M2', 'ドル円']
location = ['upper right', 'upper center', 'upper left', 'upper right',
'upper left', 'upper right']
initvalue = np.zeros([len(data), len(data.T)])
zerovalue = np.zeros([len(data), len(data.T)])
tshock = shock[:, :, 0] + shock[:, :, 1] + shock[:, :, 2] + shock[
:, :, 3] + shock[:, :, 4] + shock[:, :, 5]
for vl in range(0, 1, 1):
initvalue[:, vl] = data[:, vl] - tshock[vl, :]
plt.figure(figsize=(10, 6))
plt.title('%s' %(vnames[vl]), fontsize=18)
const = zerovalue[:, vl]
temp1 = [const[i] if (const[i] > 0 and shock[vl, i, 3] > 0) \
else const[i] if (const[i] <= 0 and shock[vl, i, 3] <= 0) \
else 0 for i in range(0, len(shock[vl, :, 3]))]
plt.bar(dti, shock[vl, :, 3], bottom=temp1, label='金融政策ショック',
width=1, color='black', edgecolor='black')
soos = shock[:, :, 0] + shock[:, :, 1] + shock[
:, :, 2] + shock[:, :, 4] + shock[:, :, 5]
temp2 = [const[i] + shock[vl, i, 3] \
if (const[i] > 0 and shock[vl, i, 3] > 0 and soos[vl, i] > 0) \
else const[i] \
if (const[i] > 0 and shock[vl, i, 3] <= 0 and soos[vl, i] > 0) \
else shock[vl, i, 3] \
if (const[i] <= 0 and shock[vl, i, 3] > 0 and soos[vl, i] > 0) \
else 0 \
if (const[i] <= 0 and shock[vl, i, 3] <= 0 and soos[vl, i] > 0) \
else const[i] + shock[vl, i, 3] \
if (const[i] <= 0 and shock[vl, i, 3] <= 0 and soos[vl, i] <= 0) \
else const[i] \
if (const[i] <= 0 and shock[vl, i, 3] > 0 and soos[vl, i] <= 0) \
else shock[vl, i, 3] \
if (const[i] > 0 and shock[vl, i, 3] <= 0 and soos[vl, i] <= 0) \
else 0 for i in range(0, len(shock[vl, :, 3]))]
plt.bar(dti, soos[vl, :], bottom=temp2, label="その他のショック",
width=1, color="lightgrey", edgecolor='lightgrey')
plt.hlines([0], 0, dlen, color='black', linewidth=1)
plt.legend(loc='%s' %(location[vl]), frameon=False, ncol = 2,
fontsize=14)
plt.tick_params(labelsize=14)
plt.xlim(0, dlen)
plt.yticks([-5, 0, 5, 10, 15])
axx = plt.gca()
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'])
plt.xlabel('年', fontsize=14, loc='center')
plt.ylabel('％', fontsize=14, loc='center')
In [ ]:
# データの読み込み
dinput = pd.read_excel('/content/drive/My Drive/Data/data_VAR7.xlsx',
sheet_name='data_q', header=[0])
data = np.zeros([len(dinput), 3])
for i in range(1, len(dinput)):
data[:, 0] = dinput.cpi / dinput.cpi.shift(4) * 100 - 100
data[:, 1] = dinput.u
data[:, 2] = dinput.RSR
data = data[4:]
data[:, 0] = data[:, 0] - data[:, 0].mean()
data[:, 1] = data[:, 1] - data[:, 1].mean()
data[:, 2] = data[:, 2] - data[:, 2].mean()
dlen = len(data)
# VARモデルの推定
model = VAR(data)
PERIOD = 25
results = model.fit(ic='bic', maxlags=8)
lag = results.k_ar
fit = results.fittedvalues
resid = results.resid
sigma = np.dot((resid.T), resid) / ((len(data) - lag) - (len(data.T) * lag))
inva = np.linalg.cholesky(sigma)
for mm in range(0, len(data.T), 1):
fit_raw = results.params.T
fit_comp_before = np.append(np.eye(len(data.T) * (lag - 1)),
np.zeros([len(data.T) * (lag - 1),
len(data.T)]), axis=1)
fit_comp = np.append(fit_raw[:, 1:(len(data.T) * lag + 1)],
fit_comp_before, axis=0)
# 歴史分解
inva_big = np.zeros([len(data.T) * lag, len(data.T)])
inva_big[0:len(data.T), :] = inva
icomp = np.append(np.eye(len(data.T)), np.zeros(
[len(data.T), (lag - 1) * len(data.T)]), axis=1)
hdshock_big = np.zeros([len(data.T), len(data.T) * lag, len(data) - lag + 1])
hdshock = np.zeros([len(data.T), len(data.T), len(data) - lag + 1])
for i in range(0, len(data.T), 1):
eps_big = np.zeros([len(data.T), len(data) - lag + 1])
eps_big[i, 1:] = (np.dot(np.linalg.inv(inva), resid.T))[i, :]
for j in range(1, len(data) - lag + 1, 1):
hdshock_big[i, :, j] = np.dot(inva_big, eps_big[:, j]) + np.dot(
fit_comp, hdshock_big[i, :, j - 1])
hdshock[i, :, j] = np.dot(icomp, hdshock_big[i, :, j])
shock = np.zeros([len(data.T), len(data), len(data.T)])
for i in range(0, len(data.T), 1):
for j in range(0, len(data.T), 1):
shock[i, :, j] = np.append(np.zeros(lag), hdshock[j, i, 1:], axis=0)
x = np.zeros([len(data) - lag, len(data.T)])
for i in range(0, lag, 1):
x = np.append(data[i:len(data) - lag + i, :], x, axis=1)
x = x[:, 0:len(x.T) - len(data.T)]
hdinit_big = np.zeros([lag * len(data.T), len(data) - lag + 1])
hdinit = np.zeros([len(data.T), len(data) - lag + 1])
hdinit_big[:, 0] = x[0, :].T
hdinit[:, 0] = np.dot(icomp, hdinit_big[:, 0])
for i in range(1, len(data) - lag + 1, 1):
hdinit_big[:, i] = np.dot(fit_comp, hdinit_big[:, i - 1])
hdinit[:, i] = np.dot(icomp, hdinit_big[:, i])
hdconst_big = np.zeros([lag * len(data.T), len(data) - lag + 1])
hdconst = np.zeros([len(data.T), len(data) - lag + 1])
cc = np.zeros(len(data.T) * lag)
cc[0:len(data.T)] = fit_raw[:, 0]
for i in range(1, len(data) - lag + 1):
hdconst_big[:, i] = cc + np.dot(fit_comp, hdconst_big[:, i - 1])
hdconst[:, i] = np.dot(icomp, hdconst_big[:, i])
hdendo = np.sum(hdshock, axis=0) + hdinit + hdconst
# 歴史分解を図示
dti = np.arange(dlen)
vnames = ['(a) インフレ率', '(b) 失業率', '(c) 金利']
location = ['upper right', 'lower right', 'upper right']
fig, ax = plt.subplots(2, 2, figsize=[20, 15])
initvalue = np.zeros([len(data), len(data.T)])
zerovalue = np.zeros([len(data), len(data.T)])
for vl in range(0, len(data.T), 1):
initvalue[:, vl] = data[:, vl] - shock[vl, :, 0] - shock[
vl, :, 1] - shock[vl, :, 2]
for vl in range(0, len(data.T), 1):
if vl == 0:
gx = 0
gy = 0
elif vl == 1:
gx = 0
gy = 1
elif vl == 2:
gx = 1
gy = 0
ax[gx,gy].set_title('%s（3変量モデル）' %(vnames[vl]), fontsize=18)
const = zerovalue[:, vl]
temp1 = [const[i] if (const[i] > 0 and shock[vl, i, 0] > 0) \
else const[i] if (const[i] <= 0 and shock[vl, i, 0] <= 0) \
else 0 for i in range(0, len(shock[vl, :, 0]))]
ax[gx,gy].bar(dti, shock[vl, :, 0], bottom=temp1, label='供給ショック',
width=1, color='lightgrey', edgecolor='lightgrey')
temp2 = [const[i] + shock[vl, i, 0] \
if (const[i] > 0 and shock[vl, i, 0] > 0 and shock[vl, i, 1] > 0) \
else const[i] \
if (const[i] > 0 and shock[vl, i, 0] <= 0 and shock[vl, i, 1] > 0) \
else shock[vl, i, 0] \
if (const[i] <= 0 and shock[vl, i, 0] > 0 and shock[vl, i, 1] > 0) \
else 0 \
if (const[i] <= 0 and shock[vl, i, 0] <= 0 and shock[vl, i, 1] > 0) \
else const[i] + shock[vl, i, 0] \
if (const[i] <= 0 and shock[vl, i, 0] <= 0 and shock[vl, i, 1] <= 0) \
else const[i] \
if (const[i] <= 0 and shock[vl, i, 0] > 0 and shock[vl, i, 1] <= 0) \
else shock[vl, i, 0] \
if (const[i] > 0 and shock[vl, i, 0] <= 0 and shock[vl, i, 1] <= 0) \
else 0 for i in range(0, len(shock[vl, :, 0]))]
ax[gx,gy].bar(dti, shock[vl, :, 1], bottom=temp2, label='需要ショック',
width=1, color='grey', edgecolor='grey')
temp3 = [const[i] + shock[vl, i, 0] + shock[vl, i, 1] \
if (const[i] > 0 and shock[vl, i, 0] > 0 and shock[vl, i, 1] > 0
and shock[vl, i, 2] > 0) \
else const[i] + shock[vl, i, 0] \
if (const[i] > 0 and shock[vl, i, 0] > 0 and shock[vl, i, 1] <= 0
and shock[vl, i, 2] > 0) \
else const[i] + shock[vl, i, 1] \
if (const[i] > 0 and shock[vl, i, 0] <= 0 and shock[vl, i, 1] > 0
and shock[vl, i, 2] > 0) \
else shock[vl, i, 0] + shock[vl, i, 1] \
if (const[i] <= 0 and shock[vl, i, 0] > 0 and shock[vl, i, 1] > 0
and shock[vl, i, 2] > 0) \
else const[i] \
if (const[i] > 0 and shock[vl, i, 0] <= 0 and shock[vl, i, 1] <= 0
and shock[vl, i, 2] > 0) \
else shock[vl, i, 0] \
if (const[i] <= 0 and shock[vl, i, 0] > 0 and shock[vl, i, 1] <= 0
and shock[vl, i, 2] > 0) \
else shock[vl, i, 1] \
if (const[i] <= 0 and shock[vl, i, 0] <= 0 and shock[vl, i, 1] > 0
and shock[vl, i, 2] > 0) \
else 0 \
if (const[i] <= 0 and shock[vl, i, 0] <= 0 and shock[vl, i, 1] <= 0
and shock[vl, i, 2] > 0) \
else const[i] + shock[vl, i, 0] + shock[vl, i, 1] \
if (const[i] <= 0 and shock[vl, i, 0] <= 0 and shock[vl, i, 1] <= 0
and shock[vl, i, 2] <= 0) \
else const[i] + shock[vl, i, 0] \
if (const[i] <= 0 and shock[vl, i, 0] <= 0 and shock[vl, i, 1] > 0
and shock[vl, i, 2] <= 0) \
else const[i] + shock[vl, i, 1] \
if (const[i] <= 0 and shock[vl, i, 0] > 0 and shock[vl, i, 1] <= 0
and shock[vl, i, 2] <= 0) \
else shock[vl, i, 0] + shock[vl, i, 1] \
if (const[i] > 0 and shock[vl, i, 0] <= 0 and shock[vl, i, 1] <= 0
and shock[vl, i, 2] <= 0) \
else const[i] \
if (const[i] <= 0 and shock[vl, i, 0] > 0 and shock[vl, i, 1] > 0
and shock[vl, i, 2] <= 0) \
else shock[vl, i, 0] \
if (const[i] > 0 and shock[vl, i, 0] <= 0 and shock[vl, i, 1] > 0
and shock[vl, i, 2] <= 0) \
else shock[vl, i, 1] \
if (const[i] > 0 and shock[vl, i, 0] > 0 and shock[vl, i, 1] <= 0
and shock[vl, i, 2] <= 0) \
else 0 for i in range(0, len(shock[vl, :, 0]))]
ax[gx,gy].bar(dti, shock[vl, :, 2], bottom=temp3, label='金融政策ショック',
width=1, color='black', edgecolor='black')
ax[gx,gy].hlines([0], 0, dlen, color='black', linewidth=1)
ax[gx,gy].legend(loc='%s' %(location[vl]), frameon=False,
ncol = 2, fontsize=14)
ax[gx,gy].tick_params(labelsize=14)
ax[gx,gy].set_xlim(0, dlen)
ax[gx,gy].set_xlabel('年', fontsize=14, loc='center')
ax[gx,gy].set_ylabel('％', fontsize=14)
axx = ax[gx,gy]
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(['1972', '1977', '1982', '1987', '1992', '1997',
'2002', '2007', '2012', '2017', '2022'])
ax[0,0].set_ylim(-10, 20)
ax[0,1].set_ylim(-3, 3)
ax[1,0].set_ylim(-8, 10)
# データの読み込み
dinput = pd.read_excel('/content/drive/My Drive/Data/data_VAR7.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:]
data[:, 0] = data[:, 0] - data[:, 0].mean()
data[:, 1] = data[:, 1] - data[:, 1].mean()
data[:, 2] = data[:, 2] - data[:, 2].mean()
data[:, 3] = data[:, 3] - data[:, 3].mean()
data[:, 4] = data[:, 4] - data[:, 4].mean()
data[:, 5] = data[:, 5] - data[:, 5].mean()
dlen = len(data)
# VARモデルの推定
model = VAR(data)
PERIOD = 25
results = model.fit(ic='bic', maxlags=12)
lag = results.k_ar
fit = results.fittedvalues
resid = results.resid
sigma = np.dot((resid.T), resid) / ((len(data) - lag) - (len(data.T) * lag))
inva = np.linalg.cholesky(sigma)
for mm in range(0, len(data.T), 1):
fit_raw = results.params.T
fit_comp_before = np.append(np.eye(len(data.T) * (lag - 1)),
np.zeros([len(data.T) * (lag - 1),
len(data.T)]), axis=1)
fit_comp = np.append(fit_raw[:, 1:(len(data.T) * lag + 1)],
fit_comp_before, axis=0)
# 歴史分解
inva_big = np.zeros([len(data.T) * lag, len(data.T)])
inva_big[0:len(data.T), :] = inva
icomp = np.append(np.eye(len(data.T)), np.zeros(
[len(data.T), (lag - 1) * len(data.T)]), axis=1)
hdshock_big = np.zeros([len(data.T), len(data.T) * lag, len(data) - lag + 1])
hdshock = np.zeros([len(data.T), len(data.T), len(data) - lag + 1])
for i in range(0, len(data.T), 1):
eps_big = np.zeros([len(data.T), len(data) - lag + 1])
eps_big[i, 1:] = (np.dot(np.linalg.inv(inva), resid.T))[i, :]
for j in range(1, len(data) - lag + 1, 1):
hdshock_big[i, :, j] = np.dot(inva_big, eps_big[:, j]) + np.dot(
fit_comp, hdshock_big[i, :, j - 1])
hdshock[i, :, j] = np.dot(icomp, hdshock_big[i, :, j])
shock = np.zeros([len(data.T), len(data), len(data.T)])
for i in range(0, len(data.T), 1):
for j in range(0, len(data.T), 1):
shock[i, :, j] = np.append(np.zeros(lag), hdshock[j, i, 1:], axis=0)
x = np.zeros([len(data) - lag, len(data.T)])
for i in range(0, lag, 1):
x = np.append(data[i:len(data) - lag + i, :], x, axis=1)
x = x[:, 0:len(x.T) - len(data.T)]
hdinit_big = np.zeros([lag * len(data.T), len(data) - lag + 1])
hdinit = np.zeros([len(data.T), len(data) - lag + 1])
hdinit_big[:, 0] = x[0, :].T
hdinit[:, 0] = np.dot(icomp, hdinit_big[:, 0])
for i in range(1, len(data) - lag + 1, 1):
hdinit_big[:, i] = np.dot(fit_comp, hdinit_big[:, i - 1])
hdinit[:, i] = np.dot(icomp, hdinit_big[:, i])
hdconst_big = np.zeros([lag * len(data.T), len(data) - lag + 1])
hdconst = np.zeros([len(data.T), len(data) - lag + 1])
cc = np.zeros(len(data.T) * lag)
cc[0:len(data.T)] = fit_raw[:, 0]
for i in range(1, len(data) - lag + 1):
hdconst_big[:, i] = cc + np.dot(fit_comp, hdconst_big[:, i - 1])
hdconst[:, i] = np.dot(icomp, hdconst_big[:, i])
hdendo = np.sum(hdshock, axis=0) + hdinit + hdconst
# 歴史分解を図示
dti = np.arange(dlen)
vnames = ['(d) インフレ率', '輸入物価インフレ率', '失業率', '金利', 'M2', 'ドル円']
location = ['upper right', 'upper center', 'upper left', 'upper right',
'upper left', 'upper right']
initvalue = np.zeros([len(data), len(data.T)])
zerovalue = np.zeros([len(data), len(data.T)])
tshock = shock[:, :, 0] + shock[:, :, 1] + shock[:, :, 2] + shock[
:, :, 3] + shock[:, :, 4] + shock[:, :, 5]
for vl in range(0, 1, 1):
initvalue[:, vl] = data[:, vl] - tshock[vl, :]
const = zerovalue[:, vl]
temp1 = [const[i] if (const[i] > 0 and shock[vl, i, 3] > 0) \
else const[i] if (const[i] <= 0 and shock[vl, i, 3] <= 0) \
else 0 for i in range(0, len(shock[vl, :, 3]))]
soos = shock[:, :, 0] + shock[:, :, 1] + shock[
:, :, 2] + shock[:, :, 4] + shock[:, :, 5]
temp2 = [const[i] + shock[vl, i, 3] \
if (const[i] > 0 and shock[vl, i, 3] > 0 and soos[vl, i] > 0) \
else const[i] \
if (const[i] > 0 and shock[vl, i, 3] <= 0 and soos[vl, i] > 0) \
else shock[vl, i, 3] \
if (const[i] <= 0 and shock[vl, i, 3] > 0 and soos[vl, i] > 0) \
else 0 \
if (const[i] <= 0 and shock[vl, i, 3] <= 0 and soos[vl, i] > 0) \
else const[i] + shock[vl, i, 3] \
if (const[i] <= 0 and shock[vl, i, 3] <= 0 and soos[vl, i] <= 0) \
else const[i] \
if (const[i] <= 0 and shock[vl, i, 3] > 0 and soos[vl, i] <= 0) \
else shock[vl, i, 3] \
if (const[i] > 0 and shock[vl, i, 3] <= 0 and soos[vl, i] <= 0) \
else 0 for i in range(0, len(shock[vl, :, 3]))]
ax[1,1].set_title('%s（6変量モデル）' %(vnames[vl]), fontsize=18)
ax[1,1].bar(dti, shock[vl, :, 3], bottom=temp1, label='金融政策ショック',
width=1, color='black', edgecolor='black')
ax[1,1].bar(dti, soos[vl, :], bottom=temp2, label="その他のショック",
width=1, color="lightgrey", edgecolor='lightgrey')
ax[1,1].hlines([0], 0, dlen, color='black', linewidth=1)
ax[1,1].legend(loc='%s' %(location[vl]), frameon=False, ncol = 2,
fontsize=14)
ax[1,1].tick_params(labelsize=14)
ax[1,1].set_xlim(0, dlen)
ax[1,1].set_xlabel('年', fontsize=14, loc='center')
ax[1,1].set_ylabel('％', fontsize=14)
axx = ax[1,1]
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[1,1].set_yticks([-5, 0, 5, 10, 15])
（練習問題）グレンジャーの因果性検定¶
- 変数：インフレ率（消費者物価指数、除く生鮮食品・エネルギー、消費税率の影響を調整済み）、名目賃金上昇率（毎月勤労統計における現金給与総額の前年比）の2変数
- 期間は1998年1月から2023年12月
- ラグ次数はBICにより13か月が選択
- 「賃金から物価」では、「グレンジャー因果性はない」との帰無仮説が棄却されたが、「物価から賃金」では、「グレンジャーの因果性はない」との帰無仮説を棄却することができなかった。
In [ ]:
# データの読み込み
dinput = pd.read_excel('/content/drive/My Drive/Data/data_VAR7.xlsx',
sheet_name='data_m', header=[0])
data = np.zeros([len(dinput), 2])
for i in range(1, len(dinput)):
data[:, 0] = dinput.infl
data[:, 1] = dinput.wage
data = data[324:]
# VARを推定
model = VAR(data)
result = model.fit(maxlags=24, ic="bic")
print(f"ラグ次数: {result.k_ar}")
# グレンジャーの因果性検定（ワルド統計量で検定、有意水準5％）
# 帰無仮説「y2（賃金上昇率）からy1（インフレ率）へグレンジャー因果性がない」
test12_outcome = result.test_causality("y1", "y2", kind="wald", signif=0.05)
print(test12_outcome.summary())
# 帰無仮説「y1（インフレ率）からy2（賃金上昇率）へグレンジャー因果性がない」
test21_outcome = result.test_causality("y2", "y1", kind="wald", signif=0.05)
print(test21_outcome.summary())
ラグ次数: 3
Granger causality Wald-test. H_0: y2 does not Granger-cause y1. Conclusion: reject H_0 at 5% significance level.
========================================
Test statistic Critical value p-value df
----------------------------------------
14.97 7.815 0.002 3
----------------------------------------
Granger causality Wald-test. H_0: y1 does not Granger-cause y2. Conclusion: fail to reject H_0 at 5% significance level.
========================================
Test statistic Critical value p-value df
----------------------------------------
6.447 7.815 0.092 3
----------------------------------------