In [ ]:
# 必要なライブラリをインストール
!pip install japanize_matplotlib
!pip install linearmodels
In [3]:
# ライブラリを起動
from dateutil.relativedelta import relativedelta
import japanize_matplotlib
import math as math
import numpy as np
import pandas as pd
import scipy as sp
from scipy import stats
import random
import statsmodels.api as sm
from statsmodels.regression.linear_model import OLS
from statsmodels.tsa.api import VAR
from linearmodels.iv import IV2SLS
from sklearn.linear_model import LinearRegression
from regex import R
from datetime import datetime as dt
import datetime
import os
from matplotlib import pyplot as plt
import matplotlib.dates as mdates
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

(9) 局所予測操作変数法と反転可能性¶

(9-1) 局所予測操作変数法による金融政策の評価¶

  • Kubota and Shintani (2022, 2025)のデータを用いて、局所予測操作変数法を行う
  • 変数：鉱工業生産指数、消費者物価指数（除く生鮮食品、消費税率の影響を調整済み）、1年物国債利回り、東証株価指数、社債スプレッドの5変数
  • 期間は1992年7月から2020年1月
  • ラグ次数は12か月
  • インパルス応答関数の信頼区間はデルタ法により計算（信頼区間の幅は90%ile）
  • 操作変数は、ターゲット因子
In [4]:
# データの読み込み
DATASET_VAR = pd.read_excel('/content/drive/My Drive/Data/data_VAR9.xlsx',
                            sheet_name='data_VAR', header=[0], index_col=0)
dti = pd.date_range('1992-07-01', '2020-01-01', freq='MS')
endo=np.zeros([len(dti),5])
endo[:, 0]=DATASET_VAR.IP
endo[:, 1]=DATASET_VAR.cCPItaxSA
endo[:, 2]=DATASET_VAR.JGB1Y
endo[:, 3]=DATASET_VAR.TOPIXMA
endo[:, 4]=DATASET_VAR.CBSm
iv=np.zeros([len(dti),1])
iv[:,0]=DATASET_VAR.TARGET
dlen = len(endo)

L = 12 + 1
MAX_H = 49
CONFIDENCE_LEVEL = 0.90
coef_result, std_result = np.zeros([MAX_H, len(endo.T)]), np.zeros(
    [MAX_H, len(endo.T)])
lbound, ubound = sp.stats.norm.interval(confidence=CONFIDENCE_LEVEL, loc=0,
                                        scale=1)
ssize = len(endo)
for h in range(0, MAX_H, 1):
    for dep in range(0, len(endo.T), 1):
        target_y = pd.DataFrame(endo[L + h - 1 : ssize, dep])
        policy_x = pd.DataFrame(endo[L - 1 : L - 1 + len(target_y), 2])
        instrument_x = pd.DataFrame(iv[L - 1 : L - 1 + len(target_y)])
        control_x = pd.DataFrame()
        for l in range(1, L, 1):
            control_x_lagged1 = pd.DataFrame(endo[L - 1 - l :
                                                  L - 1 - l + len(target_y)])
            control_x_lagged2 = pd.DataFrame(iv[L - 1 - l :
                                                L - 1 - l + len(target_y)])
            control_x = pd.concat([control_x, control_x_lagged1,
                                   control_x_lagged2], axis=1)
        model1 = LinearRegression()
        model1.fit(control_x, target_y)
        target_resid = target_y - model1.predict(control_x)
        model2 = LinearRegression()
        model2.fit(control_x, policy_x)
        policy_resid = policy_x - model2.predict(control_x)
        model3 = LinearRegression()
        model3.fit(control_x, instrument_x)
        instrument_resid = instrument_x - model3.predict(control_x)
        alpha = (target_resid.values).flatten()
        beta = (policy_resid.values).flatten()
        gamma = (instrument_resid.values).flatten()
        data = pd.DataFrame({'y': alpha, 'x': beta, 'z': gamma})
        iv_model = IV2SLS.from_formula('y ~ 1 + [x ~ z]', data=data)
        iv_res = iv_model.fit(cov_type='kernel', kernel='bartlett',
                              bandwidth = 4*((len(target_resid)/100)**(2/9)))
        coef_result[h, dep] = np.array(iv_res.params[1])
        std_result[h, dep] = np.array(iv_res.std_errors[1])
std_result[0, 2] = np.array(0)

VAR_figname = ['鉱工業生産', '消費者物価指数', '1年物国債利回り', '東証株価指数',
               '社債スプレッド']
nvars = len(VAR_figname)
fig_hrzn = np.arange(0, MAX_H)
fig = plt.figure(figsize=(nvars*4, 4))
ylim=([-40, 40], [-4, 3], [-1.5, 2.5], [-150, 100], [-2, 1.5])
for i in range(nvars):
    ax = fig.add_subplot(1, 5, i+1)
    ax.set_title(VAR_figname[i], fontsize=14)
    ax.axhline(y=0, xmin=0, xmax=1, color='black', linewidth=0.75)
    ax.plot(fig_hrzn, coef_result[:, i], color='black', linewidth=2)
    ax.plot(fig_hrzn, coef_result[:, i] + std_result[:, i] * lbound,
            color='black', linewidth=1, linestyle='dashed')
    ax.plot(fig_hrzn, coef_result[:, i] + std_result[:, i] * ubound,
            color='black', linewidth=1, linestyle='dashed')
    ax.set_xlim([0, MAX_H])
    ax.set_xlabel('月', fontsize=14)
    ax.set_xticks(range(0, MAX_H+1, 6))
    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)
    ax.set_ylim(ylim[i])
fig.tight_layout()
No description has been provided for this image

(9-2) 内部操作変数推定量によるインパルス応答関数¶

  • Kubota and Shintani (2022, 2025)のデータを用いて、内部操作変数推定量によるインパルス応答関数を計算する
  • 変数：鉱工業生産指数、消費者物価指数（除く生鮮食品、消費税率の影響を調整済み）、1年物国債利回り、東証株価指数、社債スプレッドの5変数
  • 期間は1992年7月から2020年1月
  • ラグ次数は12か月
  • インパルス応答関数の信頼区間はブートストラップにより計算（信頼区間の幅は90%ile）
  • 操作変数は、ターゲット因子
In [5]:
# データの読み込み
DATASET_VAR = pd.read_excel('/content/drive/My Drive/Data/data_VAR9.xlsx',
                            sheet_name='data_VAR', header=[0], index_col=0)
dti = pd.date_range('1992-07-01', '2020-01-01', freq='MS')
data=np.zeros([len(dti),6])
data[:, 0]=DATASET_VAR.TARGET
data[:, 1]=DATASET_VAR.JGB1Y
data[:, 2]=DATASET_VAR.IP
data[:, 3]=DATASET_VAR.cCPItaxSA
data[:, 4]=DATASET_VAR.TOPIXMA
data[:, 5]=DATASET_VAR.CBSm
dlen = len(data)

# VARモデルの推定
SIGNIFICANCE_LEVEL = 0.1
NDRAWS = 2000
PERIOD = 49
model = VAR(data)
results = model.fit(12)
lag = results.k_ar
irf = results.irf(PERIOD)
impulse = irf.orth_irfs
adj = irf.orth_irfs[0, 1, 0]
impulse_1 = impulse[:, 0][:, 0] /adj
impulse_2 = impulse[:, 1][:, 0] /adj
impulse_3 = impulse[:, 2][:, 0] /adj
impulse_4 = impulse[:, 3][:, 0] /adj
impulse_5 = impulse[:, 4][:, 0] /adj
impulse_6 = impulse[:, 5][:, 0] /adj
impulse_bs_1 = np.zeros([NDRAWS, PERIOD + 1])
impulse_bs_2 = np.zeros([NDRAWS, PERIOD + 1])
impulse_bs_3 = np.zeros([NDRAWS, PERIOD + 1])
impulse_bs_4 = np.zeros([NDRAWS, PERIOD + 1])
impulse_bs_5 = np.zeros([NDRAWS, PERIOD + 1])
impulse_bs_6 = np.zeros([NDRAWS, PERIOD + 1])

tt = 0
while tt< NDRAWS:
    bs_resid = pd.DataFrame()
    for i in range(len(data) - lag):
            random_resid = random.choice(results.resid)
            random_resid = random_resid.reshape(1, -1)
            df_resid = pd.DataFrame(random_resid)
            bs_resid = pd.concat([bs_resid, df_resid], ignore_index = True)
    bs_resid.columns= ["y1", "y2", "y3", "y4", "y5", "y6"]
    initial_data = data[ : 12]
    input_data = initial_data
    bootstrap_sample = initial_data
    for j in range(0, len(data) - lag):
        forecast_values = results.forecast(y=input_data, steps=1)
        forecast = np.array([forecast_values[0][0] + bs_resid["y1"][j],
                             forecast_values[0][1] + bs_resid["y2"][j],
                             forecast_values[0][2] + bs_resid["y3"][j],
                             forecast_values[0][3] + bs_resid["y4"][j],
                             forecast_values[0][4] + bs_resid["y5"][j],
                             forecast_values[0][5] + bs_resid["y6"][j]])
        forecast = forecast.reshape(1, -1)
        combined_data = np.vstack([input_data, forecast])
        input_data = combined_data[1 : ]
        bootstrap_sample = np.vstack([bootstrap_sample, forecast])
    model_bs = VAR(bootstrap_sample)
    results_bs = model_bs.fit(lag)
    irf_bs = results_bs.irf(PERIOD)
    impulse_bs = irf_bs.orth_irfs
    adj_bs = irf_bs.orth_irfs[0, 1, 0]
    impulse_bs_1[tt, :] = impulse_bs[:, 0][:, 0] / adj_bs
    impulse_bs_2[tt, :] = impulse_bs[:, 1][:, 0] / adj_bs
    impulse_bs_3[tt, :] = impulse_bs[:, 2][:, 0] / adj_bs
    impulse_bs_4[tt, :] = impulse_bs[:, 3][:, 0] / adj_bs
    impulse_bs_5[tt, :] = impulse_bs[:, 4][:, 0] / adj_bs
    impulse_bs_6[tt, :] = impulse_bs[:, 5][:, 0] / adj_bs
    tt=tt+1

VAR_figname = ['鉱工業生産', '消費者物価指数', '1年物国債利回り', '東証株価指数',
               '社債スプレッド']
nvars = len(VAR_figname)
counter = 1
fig_hrzn = np.arange(0, PERIOD + 1)
fig = plt.figure(figsize=(nvars*4, 4))
ylim=([-40, 40], [-4, 3], [-1.5, 2.5], [-150, 50], [-2, 1.5])
order = [1, 2, 0, 3, 4]
for i in range(nvars):
    ax = fig.add_subplot(1, 5, i+1)
    ax.set_title(VAR_figname[i], fontsize=14)
    ax.axhline(y=0, xmin=0, xmax=1, color='black', linewidth=0.75)
    ax.plot(fig_hrzn, globals()[f"impulse_{2+order[i]}"], linewidth=2,
            color='black')
    ax.plot(fig_hrzn, pd.DataFrame(globals()[f"impulse_bs_{2+order[i]}"]
                                   ).quantile(SIGNIFICANCE_LEVEL / 2),
            linewidth=1, linestyle='dashed', color='black')
    ax.plot(fig_hrzn, pd.DataFrame(globals()[f"impulse_bs_{2+order[i]}"]
                                   ).quantile(1 - SIGNIFICANCE_LEVEL / 2),
            linewidth=1, linestyle='dashed', color='black')
    ax.set_xlim([0, PERIOD])
    ax.set_xlabel('月', fontsize=14)
    ax.set_xticks(range(0, PERIOD+1, 6))
    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)
    ax.set_ylim(ylim[i])
fig.tight_layout()
No description has been provided for this image