In [ ]:
# 必要なライブラリをインストール
!pip install japanize_matplotlib
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 numpy.linalg import inv, cholesky, qr
from scipy import stats
from scipy.stats import invwishart, multivariate_normal
import random
import statsmodels.api as sm
from statsmodels.regression.linear_model import OLS
from statsmodels.tsa.api import VAR
from sklearn.linear_model import LinearRegression
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

(10) 符号制約とベイズ推論¶

符号制約構造VARモデルによる金融政策の評価¶

  • 変数：実質GDP（対数値）、インフレ率（除く食料・エネルギー、消費税率の影響を調整済み）、10年物国債利回り、名目実効為替レート（対数値）、日経平均株価（対数値）、銀行貸出残高（前年比）
  • 期間は2007年1月から2024年12月
  • ラグ次数はAICにより2か月
  • 事前分布にはミネソタ事前分布を用いる
  • ショックの識別は、金融引締ショックの発生後12か月間にわたり、10年物国債利回りの正（上昇）の反応が継続し、かつインフレ率の負（低下）の反応が継続する、という符号制約
  • 金融引締めショックに対するインパルス応答関数の事後メジアンと 68%（等裾）信用区間が示されている（対数変換されていた実質GDP、名目実効為替レート、日経平均株価は％表示に変換済み）
In [5]:
# データの読み込み
dinput = pd.read_excel('/content/drive/My Drive/Data/data_VAR10.xlsx',
                       sheet_name='data_m', header=[0], index_col=0)
Y = dinput[['GDP','CPI','LTR','EXR','NIKKEI','LENDING']].values
T, k = Y.shape
model = VAR(Y)
res = model.fit(maxlags=12, ic='aic')
p = res.k_ar

# VARモデルの推定
X = []
for i in range(p):
    X.append(Y[p-1-i:T-1-i])
X = np.hstack(X)
X = np.hstack([np.ones((T-p,1)), X])
Y_trim = Y[p:]
B_ols = inv(X.T @ X) @ X.T @ Y_trim
U = Y_trim - X @ B_ols
S_post = U.T @ U
nu_post = T - p

# ミネソタ事前分布の設定
sigma_sq = np.var(U, axis=0)
lambda_ = 0.1
alpha = 2.0
def minnesota_prior():
    ncoef = 1 + k * p
    b_prior = np.zeros((ncoef, k))
    V_prior = np.zeros((ncoef, ncoef))
    V_prior[0, 0] = 1e6
    for l in range(1, p+1):
        for i in range(k):
            for j in range(k):
                idx = 1 + (l-1)*k + j
                if i == j and l == 1:
                    b_prior[idx, i] = 1.0
                var_ij = (lambda_**2) / (l**(2*alpha)) * (sigma_sq[i]
                                                          / sigma_sq[j])
                V_prior[idx, idx] = var_ij
    return b_prior, V_prior
b_prior, V_prior = minnesota_prior()

# 関数の設定
def draw_sigma():
    return invwishart.rvs(df=nu_post, scale=S_post)

def draw_beta(Sigma):
    ncoef = X.shape[1]
    beta_draw = np.zeros((ncoef, k))
    XX = X.T @ X
    XY = X.T @ Y_trim
    for i in range(k):
        sigma_i = Sigma[i, i]
        V_post_i = inv(inv(V_prior) + XX / sigma_i)
        b_post_i = V_post_i @ (
            inv(V_prior) @ b_prior[:, i] + XY[:, i] / sigma_i
        )
        beta_draw[:, i] = (
            b_post_i
            + cholesky(V_post_i) @ np.random.randn(ncoef)
        )
    return beta_draw

def beta_to_A(beta):
    A = np.zeros((p, k, k))
    for i in range(p):
        A[i] = beta[1+i*k:1+(i+1)*k, :].T
    return A

def random_orthogonal(k):
    Z = np.random.randn(k, k)
    Q, R = qr(Z)
    return Q

def compute_irfs(A, B, horizons):
    p, k, _ = A.shape
    Phi = np.zeros((horizons+1, k, k))
    Phi[0] = np.eye(k)
    for h in range(1, horizons+1):
        for j in range(1, min(p, h)+1):
            Phi[h] += A[j-1] @ Phi[h-j]
    return np.array([Phi[h] @ B for h in range(horizons+1)])

def check_sign_restrictions(irfs, horizon=12):
    shock = 2
    for h in range(1, horizon+1):
        if irfs[h, 1, shock] <= 0:
            return False
        if irfs[h, 2, shock] >= 0:
            return False
    return True

def scale_irfs_to_10bps(irfs, shock=2, policy_index=2, target=0.10):
    impact_response = irfs[0, policy_index, shock]
    scale = target / impact_response
    return irfs * scale

PERIOD = 120
TARGET_ACCEPT = 5000
MAX_TRIES = 500000
accepted_irfs = []
n_try = 0
np.random.seed(0)
while (len(accepted_irfs) < TARGET_ACCEPT) and (n_try < MAX_TRIES):
    n_try += 1
    # ステップ1: 分散共分散行列と係数行列の抽出
    Sigma = draw_sigma()
    beta = draw_beta(Sigma)
    A = beta_to_A(beta)
    # ステップ2: 直交行列の抽出
    P = cholesky(Sigma)
    Q = random_orthogonal(k)
    B = P @ Q
    # ステップ3: 符号制約のチェック
    irfs = compute_irfs(A, B, PERIOD)
    if check_sign_restrictions(irfs):
        accepted_irfs.append(irfs)
    # ステップ4: 繰り返しと集計
    if n_try % 10000 == 0:
        print(f"tries={n_try}, accepted={len(accepted_irfs)}")

accepted_irfs = np.array(accepted_irfs)
print("Accepted draws:", accepted_irfs.shape[0])

scaled_irfs = []
for irfs in accepted_irfs:
    scaled_irfs.append(
        scale_irfs_to_10bps(irfs)
    )
scaled_irfs = np.array(scaled_irfs)
median_irf = np.median(scaled_irfs, axis=0)
lower_irf  = np.percentile(scaled_irfs, 16, axis=0)
upper_irf  = np.percentile(scaled_irfs, 84, axis=0)

# グラフの描画
VAR_figname = ['実質GDP', 'インフレ率', '10年物国債利回り', '名目実効為替レート',
               '日経平均株価', '銀行貸出残高']
nvars = len(VAR_figname)
fig_hrzn = np.arange(0, PERIOD+1)
fig = plt.figure(figsize=(nvars*2.5, 12))
ylim=([-2.5, 1], [-0.6, 0.1], [-0.15, 0.15], [-3, 5], [-15, 10], [-1, 0.8])
taisu = ([100, 1, 1, 100, 100, 1])
for i in range(nvars):
    ax = fig.add_subplot(3, 3, 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, median_irf[:, i, 2] * taisu[i], color='black',
            linewidth=2)
    ax.plot(fig_hrzn, np.ravel(lower_irf[:, i, 2] * taisu[i]),
            color='black', linewidth=1, linestyle='dashed')
    ax.plot(fig_hrzn, np.ravel(upper_irf[:, i, 2] * taisu[i]),
            color='black', linewidth=1, linestyle='dashed')
    ax.set_xlim([0, PERIOD])
    ax.set_xlabel('月', fontsize=14)
    ax.set_xticks(range(0, PERIOD+1, 12))
    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)
    ax.set_ylim(ylim[i])
fig.tight_layout()

tries=10000, accepted=2653
Accepted draws: 5000
No description has been provided for this image

(練習問題2) 従来の再帰的識別制約（スロー・R・ファスト制約）を課した場合の金融政策ショックに対する6変数のインパルス応答関数¶

  • 用いる変数や推計期間などは符号制約構造VARモデルのケースと同じ
  • ショックの識別を従来の再帰的識別制約に変更
In [6]:
# データの読み込み
dinput = pd.read_excel('/content/drive/My Drive/Data/data_VAR10.xlsx',
                       sheet_name='data_m', header=[0], index_col=0)
data = dinput[['GDP','CPI','LTR','EXR','NIKKEI','LENDING']].values
dlen = len(data)
dti = pd.date_range('2007-07-01', '2024-12-01', freq='MS')

PERIOD = 120
CONFIDENCE_LEVEL = 0.68
lbound, ubound = sp.stats.norm.interval(confidence=CONFIDENCE_LEVEL, loc=0,
                                        scale=1)
model = VAR(data)
result = model.fit(maxlags=12, ic='aic')
irf = result.irf(PERIOD)
impulse = irf.orth_irfs
scaling = 1/impulse[0, 2, 2]*(0.1)
scaled_irfs = impulse * scaling
lband = scaled_irfs + irf.stderr(orth=True) * lbound * scaling
uband = scaled_irfs + irf.stderr(orth=True) * ubound * scaling

# グラフの描画
nvars = len(VAR_figname)
fig_hrzn = np.arange(0, PERIOD+1)
fig = plt.figure(figsize=(nvars*2.5, 12))
ylim=([-2.5, 1], [-0.6, 0.1], [-0.15, 0.15], [-3, 5], [-15, 10], [-1, 0.8])
taisu = ([100, 1, 1, 100, 100, 1])
for i in range(nvars):
    ax = fig.add_subplot(3, 3, 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, scaled_irfs[:, i, 2] * taisu[i], color='black',
            linewidth=2)
    ax.plot(fig_hrzn, np.ravel(lband[:, i, 2] * taisu[i]),
            color='black', linewidth=1, linestyle='dashed')
    ax.plot(fig_hrzn, np.ravel(uband[:, i, 2] * taisu[i]),
            color='black', linewidth=1, linestyle='dashed')
    ax.set_xlim([0, PERIOD])
    ax.set_xlabel('月', fontsize=14)
    ax.set_xticks(range(0, PERIOD+1, 12))
    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