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
from numpy.linalg import inv, cholesky, qr, eigvals
from scipy import stats
from scipy.stats import f, norm, invwishart, multivariate_normal, chi2
from scipy.optimize import minimize
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
import matplotlib.patches as mpatches
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

(12) 変化する経済構造の下での政策評価¶

(12-1) 平均経済成長率の低下¶

  • 変数：実質GDP（前年比、暦年値）
  • 期間は1956年から2025年
In [ ]:
dinput = pd.read_excel('/content/drive/My Drive/Data/data_VAR12.xlsx',
                       sheet_name='data_y', header=0, index_col=0)
y = dinput['GDP'].squeeze().dropna()
series1  = dinput['1956-1973'].dropna()
series2  = dinput['1974-2025'].dropna()
y = y[(y.index >= 1956) & (y.index <= 2025)]
plt.figure(figsize=(12,6))
plt.plot(y.index, y.values, linewidth=3, label='GDP成長率', color='black')
plt.plot(series1.index, series1.values, linestyle='--', linewidth=3,
         label='1956年-1973年平均', color='grey')
plt.plot(series2.index, series2.values, linestyle='dotted', linewidth=3,
         label='1974年-2025年平均', color='grey')
plt.axhline(y=0, linestyle='-', linewidth=1, color='black')
plt.xlim(1956, 2025)
plt.title('GDP成長率', fontsize=16)
plt.legend(fontsize=14, frameon=False)
plt.xlabel('年', fontsize=14)
plt.xticks(range(1956, 2026, 6), fontsize=14)
plt.yticks(range(-8, 15, 2), fontsize=14)
plt.show()
No description has been provided for this image

(12-2) 局所予測法（LP）によるSTVARモデルの推定¶

  • 被説明変数：実質GDP（対数値）
  • 説明変数（金融政策ショック）：高頻度識別による金融政策サプライズ（Kubota and Shintaniのターゲット因子）、
  • 説明変数（コントロール変数）：実質GDP（対数値）、GDPデフレータ（前年同期比）、失業率（季節調整値）、マネタリーベース（対数値）
  • 期間は1994年第1四半期から2025年第4四半期
  • 金融政策ショックに対する実質GDPのインパルス応答関数と90%信頼区間（HAC標準偏差）
  • LPの基準化定数は、金融政策ショックの1標準偏差
  • LPのラグは4四半期
  • 推移変数は内閣府のGDPギャップの1期ラグ
  • 推移変数の閾値：内閣府のGDPギャップ（1期ラグ）が閾値よりも大きい景気拡張期のサンプルが全体の約３割となる（c=-1.3）、$\gamma$=8
In [ ]:
TRANSITION_VAR = "GAPCAO"
TRANSITION_TRANSFORM = "lag1"
THRESHOLD_TYPE = "quantile"
THRESHOLD_VALUE = 0
THRESHOLD_QUANTILE = 0.7
TARGET_VARS = ["GDP"]
SHOCK_VAR    = "TARGET"
variables = ["GDP","UNEMP","MB","GDPdef"]
gamma = 8
lags = 4
H = 12
CONFIDENCE_LEVEL = 0.90

dinput = pd.read_excel('/content/drive/My Drive/Data/data_VAR12.xlsx',
                       sheet_name='data_q', header=0, index_col=0)
df = dinput.copy()

def convert_to_quarter_index(idx):
    idx = idx.astype(str)
    years, quarters = [], []
    for x in idx:
        y, q = x.split(".")
        years.append(int(y))
        quarters.append(int(q))
    return pd.PeriodIndex(year=years, quarter=quarters, freq="Q")
df.index = convert_to_quarter_index(df.index)
df_plot_index = df.index.to_timestamp()

def get_transition_series(df, var, transform):
    s = df[var]
    if transform.startswith("lag"):
        k = int(transform.replace("lag", ""))
        z = s.shift(k)
    elif transform.startswith("ma"):
        k = int(transform.replace("ma", ""))
        z = s.shift(1).rolling(window=k).mean()
    else:
        raise ValueError(f"未対応のtransform: {transform}")
    return z, k
z,k = get_transition_series(df, TRANSITION_VAR, TRANSITION_TRANSFORM)

def get_threshold(z, threshold_type, threshold_value=None,
                  threshold_quantile=None):
    z_valid = z.dropna()
    if threshold_type == "mean":
        return z_valid.mean()
    elif threshold_type == "value":
        return threshold_value
    elif threshold_type == "quantile":
        return z_valid.quantile(1.0 - threshold_quantile)
    else:
        raise ValueError(f"未対応のthreshold_type: {threshold_type}")
c = get_threshold(z, THRESHOLD_TYPE, THRESHOLD_VALUE, THRESHOLD_QUANTILE)

print(f"推移変数: {TRANSITION_VAR} ({TRANSITION_TRANSFORM}), 閾値: {c:.4f}")
print(
    f"閾値より大きいサンプル数（好況期）: {(z > c).sum()} / {z.notna().sum()} "
    f"({(z > c).sum() / z.notna().sum() * 100:.1f}%)"
)
print(
    f"閾値以下のサンプル数（不況期）: {(z <= c).sum()} / {z.notna().sum()} "
    f"({(z <= c).sum() / z.notna().sum() * 100:.1f}%)"
)

z_std = (z - z.mean()) / z.std()
c_std = (c - z.mean()) / z.std()
df["F"] = 1 / (1 + np.exp(-gamma * (z_std - c_std)))
for var in variables:
    for l in range(1, lags + 1):
        df[f"{var}_lag{l}"] = df[var].shift(l)
for l in range(1, lags + 1):
    df[f"SHOCK_lag{l}"] = df[SHOCK_VAR].shift(l)
df["shock_high"] = df[SHOCK_VAR] * df["F"]
df["shock_low"]  = df[SHOCK_VAR] * (1 - df["F"])
controls_base = ([f"{v}_lag{l}" for v in variables for l in range(1, lags + 1)])
for col in controls_base:
    df[f"{col}_high"] = df[col] * df["F"]
    df[f"{col}_low"]  = df[col] * (1 - df["F"])
controls_tv = [f"{col}_high" for col in controls_base] + [
    f"{col}_low"  for col in controls_base]

def estimate_lp_stvar(df_use, y_var, h, controls_tv):
    y_h = df_use[y_var].shift(-h)
    X = df_use[["shock_high", "shock_low"] + controls_tv]
    data = pd.concat([y_h, X], axis=1).dropna()
    y_reg = data.iloc[:, 0]
    X_reg = sm.add_constant(data.iloc[:, 1:])
    res = sm.OLS(y_reg, X_reg).fit(
        cov_type="HAC", cov_kwds={"maxlags": h + 1, "kernel": "bartlett"}
    )
    return res

def estimate_lp(df_use, y_var, h, controls_base):
    y_h = df_use[y_var].shift(-h)
    X = df_use[[SHOCK_VAR] + controls_base]
    data = pd.concat([y_h, X], axis=1).dropna()
    y_reg = data.iloc[:, 0]
    X_reg = sm.add_constant(data.iloc[:, 1:])
    res = sm.OLS(y_reg, X_reg).fit(
        cov_type="HAC", cov_kwds={"maxlags": h + 1, "kernel": "bartlett"}
    )
    return res

irf_results = {}
for y_var in TARGET_VARS:
    print(f"推定中: {y_var}")
    irf_h_pt  = np.zeros(H + 1)
    irf_l_pt  = np.zeros(H + 1)
    irf_n_pt  = np.zeros(H + 1)
    se_high   = np.zeros(H + 1)
    se_low    = np.zeros(H + 1)
    se_normal = np.zeros(H + 1)
    for h in range(H + 1):
        res = estimate_lp_stvar(df, y_var, h, controls_tv)
        irf_h_pt[h] = res.params["shock_high"]
        irf_l_pt[h] = res.params["shock_low"]
        se_high[h]  = res.bse["shock_high"]
        se_low[h]   = res.bse["shock_low"]
        res_s = estimate_lp(df, y_var, h, controls_base)
        irf_n_pt[h]   = res_s.params[SHOCK_VAR]
        se_normal[h]  = res_s.bse[SHOCK_VAR]
    lbound, ubound = sp.stats.norm.interval(confidence=CONFIDENCE_LEVEL,
                                            loc=0, scale=1)
    Z_CRIT = ubound
    irf_results[y_var] = {
        "high":           irf_h_pt,
        "low":            irf_l_pt,
        "normal":         irf_n_pt,
        "ci_high_low":    irf_h_pt - Z_CRIT * se_high,
        "ci_high_high":   irf_h_pt + Z_CRIT * se_high,
        "ci_low_low":     irf_l_pt - Z_CRIT * se_low,
        "ci_low_high":    irf_l_pt + Z_CRIT * se_low,
        "ci_normal_low":  irf_n_pt - Z_CRIT * se_normal,
        "ci_normal_high": irf_n_pt + Z_CRIT * se_normal,
    }

threshold_desc = {
    "mean": "平均",
    "value": f"指定値={THRESHOLD_VALUE}",
    "quantile": f"上位{int(THRESHOLD_QUANTILE*100)}%分位",
}[THRESHOLD_TYPE]
fig_regime, ax_f = plt.subplots(figsize=(12, 5), constrained_layout=True)
f_plot = df["F"].dropna()
f_plot_ts = f_plot.copy()
f_plot_ts.index = f_plot.index.to_timestamp()
ax_f.plot(f_plot_ts.index, 1 - f_plot_ts.shift(-1).values, color="black",
          linewidth=1.5, label="G(z)")
keiki_plot = df["KEIKI"]
keiki_plot.index = keiki_plot.index.to_timestamp()
ax_f.fill_between(keiki_plot.index, 0, 1, where=(keiki_plot == 1),
                  color="gray", alpha=0.25, label="景気基準日付（景気後退期）")
ax_f.set_ylim(-0.01, 1.01)
ax_f.set_ylabel("推移関数の重み G(z)", fontsize=14)
ax_f.tick_params(axis='both', labelsize=12)
ax_f.legend(loc="upper center", bbox_to_anchor=(0.5, -0.08), ncol=2,
            fontsize=14, frameon=False)
ax_f.spines["right"].set_visible(False)
ax_f.spines["top"].set_visible(False)
ax_f.xaxis.set_major_locator(mdates.YearLocator(2))
ax_f.xaxis.set_major_formatter(mdates.DateFormatter('%Y'))
ax_f.set_xlim(pd.Timestamp("1994-01-01"), pd.Timestamp("2025-12-31"))
plt.show()

VAR_figname_all = {"GDP": "実質GDP"}
fig_hrzn = np.arange(0, H + 1)
std = np.std(df[SHOCK_VAR], ddof=1)
colors = {"high": "black", "low": "black", "normal": "black"}
for y_var in TARGET_VARS:
    res = {k: v * std * 100 for k, v in irf_results[y_var].items()}
    var_title = VAR_figname_all.get(y_var, y_var)
    fig, axes = plt.subplots(1, 3, figsize=(12, 5), constrained_layout=True)
    ax = axes[0]
    ax.plot(fig_hrzn, res["normal"], color=colors["normal"], linewidth=2,
            label="LP")
    ax.plot(fig_hrzn, res["ci_normal_low"], color=colors["normal"],
            linestyle="--", linewidth=1)
    ax.plot(fig_hrzn, res["ci_normal_high"], color=colors["normal"],
            linestyle="--", linewidth=1)
    ax.axhline(0, color="black", linewidth=0.75)
    ax.tick_params(axis='both', labelsize=14)
    ax.set_title("線形モデル", fontsize=16)
    ax.set_xlabel("四半期", fontsize=14)
    ax.set_xlim([0, H])
    ax.set_ylim([-1.4, 0.8])
    ax.set_yticks(np.arange(-1.4, 0.9, 0.2))
    ax.set_xticks(range(0, H + 1, 2))
    ax.spines["right"].set_visible(False); ax.spines["top"].set_visible(False)
    ax = axes[1]
    ax.plot(fig_hrzn, res["high"], color=colors["high"], linewidth=2)
    ax.plot(fig_hrzn, res["ci_high_low"], color=colors["high"],
            linestyle="--", linewidth=1)
    ax.plot(fig_hrzn, res["ci_high_high"], color=colors["high"],
            linestyle="--", linewidth=1)
    ax.axhline(0, color="black", linewidth=0.75)
    ax.tick_params(axis='both', labelsize=14)
    ax.set_title("拡張期", fontsize=16)
    ax.set_xlabel("四半期", fontsize=14)
    ax.set_xlim([0, H])
    ax.set_ylim([-1.4, 0.8])
    ax.set_yticks(np.arange(-1.4, 0.9, 0.2))
    ax.set_xticks(range(0, H + 1, 2))
    ax.spines["right"].set_visible(False); ax.spines["top"].set_visible(False)
    ax = axes[2]
    ax.plot(fig_hrzn, res["low"], color=colors["low"], linewidth=2)
    ax.plot(fig_hrzn, res["ci_low_low"], color=colors["high"],
            linestyle="--", linewidth=1)
    ax.plot(fig_hrzn, res["ci_low_high"], color=colors["high"],
            linestyle="--", linewidth=1)
    ax.axhline(0, color="black", linewidth=0.75)
    ax.tick_params(axis='both', labelsize=14)
    ax.set_title("後退期", fontsize=16)
    ax.set_xlabel("四半期", fontsize=14)
    ax.set_xlim([0, H])
    ax.set_ylim([-1.4, 0.8])
    ax.set_yticks(np.arange(-1.4, 0.9, 0.2))
    ax.set_xticks(range(0, H + 1, 2))
    ax.spines["right"].set_visible(False); ax.spines["top"].set_visible(False)
    plt.show()

推移変数: GAPCAO (lag1), 閾値: -1.3000
閾値より大きいサンプル数（好況期）: 85 / 127 (66.9%)
閾値以下のサンプル数（不況期）: 42 / 127 (33.1%)
推定中: GDP
No description has been provided for this image
No description has been provided for this image

(練習問題2-1) 既知の構造変化点$\tau=1973$についてLR検定¶

In [ ]:
dinput = pd.read_excel('/content/drive/My Drive/Data/data_VAR12.xlsx',
                       sheet_name='data_y', header=0, index_col=0)
y = dinput['GDP']
break_year = 1973
y1 = y[y.index <= break_year]
y2 = y[y.index > break_year]

def estimate_ar1(y_series):
    y_dep = y_series.iloc[1:]
    y_lag = y_series.shift(1).iloc[1:]
    X = sm.add_constant(y_lag)
    model = sm.OLS(y_dep, X).fit()
    SSR = np.sum(model.resid**2)
    return model, SSR

model_full, SSR_full = estimate_ar1(y)
model1, SSR1 = estimate_ar1(y1)
model2, SSR2 = estimate_ar1(y2)
SSR_unrestricted = SSR1 + SSR2
SSR_restricted = SSR_full
n = len(y) - 1
LR = n * np.log(SSR_restricted / SSR_unrestricted)
df = 2
p_value = 1 - chi2.cdf(LR, df)
print("===== LR test for AR(1) structural break =====")
print(f"Break year       : {break_year}")
print(f"LR statistic     : {LR:.4f}")
print(f"Degrees freedom  : {df}")
print(f"P-value          : {p_value:.4f}")

print("\n===== AR(1) coefficients =====")
print("--- 1956-1990 ---")
print(f"Constant : {model1.params.iloc[0]:.4f}")
print(f"AR(1)    : {model1.params.iloc[1]:.4f}")
print("--- 1991-2025 ---")
print(f"Constant : {model2.params.iloc[0]:.4f}")
print(f"AR(1)    : {model2.params.iloc[1]:.4f}")

===== LR test for AR(1) structural break =====
Break year       : 1973
LR statistic     : 28.7219
Degrees freedom  : 2
P-value          : 0.0000

===== AR(1) coefficients =====
--- 1956-1990 ---
Constant : 7.2997
AR(1)    : 0.2202
--- 1991-2025 ---
Constant : 1.1683
AR(1)    : 0.4131

(練習問題2-2) 既知の構造変化点$\tau=1990$についてLR検定¶

In [ ]:
dinput = pd.read_excel('/content/drive/My Drive/Data/data_VAR12.xlsx',
                       sheet_name='data_y', header=0, index_col=0)
y = dinput['GDP']
break_year = 1990
y1 = y[y.index <= break_year]
y2 = y[y.index > break_year]

def estimate_ar1(y_series):
    y_dep = y_series.iloc[1:]
    y_lag = y_series.shift(1).iloc[1:]
    X = sm.add_constant(y_lag)
    model = sm.OLS(y_dep, X).fit()
    SSR = np.sum(model.resid**2)
    return model, SSR

model_full, SSR_full = estimate_ar1(y)
model1, SSR1 = estimate_ar1(y1)
model2, SSR2 = estimate_ar1(y2)
SSR_unrestricted = SSR1 + SSR2
SSR_restricted = SSR_full
n = len(y) - 1
LR = n * np.log(SSR_restricted / SSR_unrestricted)
df = 2
p_value = 1 - chi2.cdf(LR, df)
print("===== LR test for AR(1) structural break =====")
print(f"Break year       : {break_year}")
print(f"LR statistic     : {LR:.4f}")
print(f"Degrees freedom  : {df}")
print(f"P-value          : {p_value:.4f}")

print("\n===== AR(1) coefficients =====")
print("--- 1956-1990 ---")
print(f"Constant : {model1.params.iloc[0]:.4f}")
print(f"AR(1)    : {model1.params.iloc[1]:.4f}")
print("--- 1991-2025 ---")
print(f"Constant : {model2.params.iloc[0]:.4f}")
print(f"AR(1)    : {model2.params.iloc[1]:.4f}")

===== LR test for AR(1) structural break =====
Break year       : 1990
LR statistic     : 23.3670
Degrees freedom  : 2
P-value          : 0.0000

===== AR(1) coefficients =====
--- 1956-1990 ---
Constant : 2.4178
AR(1)    : 0.6320
--- 1991-2025 ---
Constant : 0.8632
AR(1)    : -0.1101

(練習問題2-3) 未知の構造変化点$\tau$について上限LR検定¶

In [ ]:
dinput = pd.read_excel('/content/drive/My Drive/Data/data_VAR12.xlsx',
                       sheet_name='data_y', header=0, index_col=0)
y = dinput['GDP'].squeeze().dropna()
series1  = dinput['1956-1973'].dropna()
series2  = dinput['1974-2025'].dropna()
y = y[(y.index >= 1956) & (y.index <= 2025)]
trim = 0
n = len(y)
start = int(n * trim)
end   = int(n * (1 - trim))
candidate_breaks = y.index[start:end]

def estimate_ar1(y_series):
    y_dep = y_series.iloc[1:]
    y_lag = y_series.shift(1).iloc[1:]
    X = sm.add_constant(y_lag)
    model = sm.OLS(y_dep, X).fit()
    SSR = np.sum(model.resid**2)
    return model, SSR

model_full, SSR_full = estimate_ar1(y)
LR_list = []
break_list = []
for b in candidate_breaks:
    y1 = y[y.index <= b]
    y2 = y[y.index > b]
    if len(y1) < 10 or len(y2) < 10:
        continue
    model1, SSR1 = estimate_ar1(y1)
    model2, SSR2 = estimate_ar1(y2)
    SSR_unrestricted = SSR1 + SSR2
    n_effective = len(y) - 1
    LR = n_effective * np.log(SSR_full / SSR_unrestricted)
    LR_list.append(LR)
    break_list.append(b)
LR_array = np.array(LR_list)
QLR = np.max(LR_array)
best_break = break_list[np.argmax(LR_array)]
df = 2
p_value = 1 - chi2.cdf(QLR, df=df)
print("===== Quandt LR Test for AR(1) =====")
print(f"Sample            : 1956-2025")
print(f"Estimated break   : {best_break}")
print(f"Pre-break sample  : 1956-{best_break}")
print(f"Post-break sample : {best_break + 1}-2025")
print(f"Sup LR statistic  : {QLR:.4f}")
print(f"Approx p-value    : {p_value:.4f}")
y1 = y[y.index <= best_break]
y2 = y[y.index > best_break]
model1, _ = estimate_ar1(y1)
model2, _ = estimate_ar1(y2)
print("\n===== AR(1) estimates =====")
print(f"--- Pre-break (<= {best_break}) ---")
print(f"Constant : {model1.params.iloc[0]:.4f}")
print(f"AR(1)    : {model1.params.iloc[1]:.4f}")
print(f"\n--- Post-break (> {best_break}) ---")
print(f"Constant : {model2.params.iloc[0]:.4f}")
print(f"AR(1)    : {model2.params.iloc[1]:.4f}")

===== Quandt LR Test for AR(1) =====
Sample            : 1956-2025
Estimated break   : 1973
Pre-break sample  : 1956-1973
Post-break sample : 1974-2025
Sup LR statistic  : 28.7219
Approx p-value    : 0.0000

===== AR(1) estimates =====
--- Pre-break (<= 1973) ---
Constant : 7.2997
AR(1)    : 0.2202

--- Post-break (> 1973) ---
Constant : 1.1683
AR(1)    : 0.4131