from __future__ import annotations import json import math from pathlib import Path import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt import numpy as np import pandas as pd import statsmodels.api as sm from matplotlib import font_manager from statsmodels.stats.multitest import multipletests from statsmodels.tsa.stattools import adfuller ROOT = Path(__file__).resolve().parent DATA = ROOT / "data" OUTPUTS = ROOT / "outputs" FIGURES = OUTPUTS / "figures" OUTPUTS.mkdir(exist_ok=True) FIGURES.mkdir(exist_ok=True) FONT = ROOT.parents[1] / "public" / "fonts" / "NotoSansTC-Regular.ttf" if FONT.exists(): font_manager.fontManager.addfont(str(FONT)) plt.rcParams["font.family"] = "Noto Sans TC" plt.rcParams["axes.unicode_minus"] = False plt.rcParams["figure.dpi"] = 180 plt.rcParams["savefig.facecolor"] = "white" COLORS = { "brown": "#6F3519", "orange": "#C56E33", "gold": "#B88A3B", "green": "#3F6B57", "blue": "#3C6682", "red": "#A54842", "purple": "#6B5B7A", "cream": "#F5EEE6", "ink": "#2D211B", "muted": "#756B64", } OUTCOME_SPECS = { "fao_cereals_index": { "label": "FAO國際穀物", "change": "dlog_fao_cereals_index", "channels": ["dlog_brent_usd_bbl"], }, "import_cereals_usd": { "label": "進口穀類（美元）", "change": "dlog_import_cereals_usd", "channels": ["dlog_fao_cereals_index", "dlog_brent_usd_bbl"], }, "import_cereals_twd": { "label": "進口穀類（新臺幣）", "change": "dlog_import_cereals_twd", "channels": [ "dlog_fao_cereals_index", "dlog_brent_usd_bbl", "fx_depreciation", ], }, "cpi_food": { "label": "臺灣食品CPI", "change": "dlog_cpi_food", "channels": [ "dlog_fao_food_index", "dlog_import_cereals_twd", "dlog_brent_usd_bbl", "fx_depreciation", ], }, "cpi_grain_products": { "label": "臺灣穀類製品CPI", "change": "dlog_cpi_grain_products", "channels": [ "dlog_fao_cereals_index", "dlog_import_cereals_twd", "dlog_brent_usd_bbl", "fx_depreciation", ], }, "cpi_rice_products": { "label": "臺灣米類製品CPI", "change": "dlog_cpi_rice_products", "channels": [ "dlog_fao_cereals_index", "dlog_import_cereals_twd", "dlog_brent_usd_bbl", "fx_depreciation", ], }, } def write_csv(frame: pd.DataFrame, name: str) -> None: frame.to_csv(OUTPUTS / name, index=False, encoding="utf-8-sig") def add_month_dummies(frame: pd.DataFrame, design: pd.DataFrame) -> pd.DataFrame: dummies = pd.get_dummies( frame.loc[design.index, "month"].astype(int), prefix="month", drop_first=True, dtype=float, ) return pd.concat([design, dummies], axis=1) def fit_hac(frame: pd.DataFrame, target: str, regressors: list[str]): selected = frame[[target, *regressors, "month"]].replace([np.inf, -np.inf], np.nan).dropna() design = selected[regressors].astype(float) design = add_month_dummies(selected, design) design = sm.add_constant(design, has_constant="add") model = sm.OLS(selected[target].astype(float), design).fit( cov_type="HAC", cov_kwds={"maxlags": 12, "use_correction": True}, ) return model, selected def linear_combination(model, names: list[str], scale: float = 1.0) -> dict[str, float]: vector = np.zeros(len(model.params)) positions = {name: index for index, name in enumerate(model.params.index)} for name in names: vector[positions[name]] = scale test = model.t_test(vector) interval = np.asarray(test.conf_int()).reshape(-1) return { "estimate": float(np.asarray(test.effect).reshape(-1)[0]), "standard_error": float(np.asarray(test.sd).reshape(-1)[0]), "p_value": float(np.asarray(test.pvalue).reshape(-1)[0]), "ci_low": float(interval[0]), "ci_high": float(interval[1]), } def event_cluster_starts(frame: pd.DataFrame) -> pd.Series: prior = ( frame["stock_tail_loss"] .shift(1) .rolling(3, min_periods=1) .max() .fillna(0) ) return ((frame["stock_tail_loss"] == 1) & (prior == 0)).astype(int) def classify_event(row: pd.Series) -> str: fx = row["fx_depreciation"] oil = row["dlog_brent_usd_bbl"] grain = row["dlog_fao_cereals_index"] if fx >= 0.02: return "新臺幣貶值達2%以上" if oil > 0 and grain > 0: return "Brent與FAO穀物同漲" if oil < 0 and grain < 0: return "Brent與FAO穀物同跌" return "混合型" def descriptive_event_table(frame: pd.DataFrame) -> pd.DataFrame: event_frame = frame.loc[frame["tail_event_start"] == 1].copy() event_frame["event_type"] = event_frame.apply(classify_event, axis=1) columns = [ "date", "event_type", "taiex_log_return", "taiex_drawdown_24m", "fx_depreciation", "dlog_brent_usd_bbl", "dlog_fao_cereals_index", "dlog_import_cereals_twd", "dlog_cpi_food", "dlog_cpi_grain_products", "dlog_cpi_rice_products", ] output = event_frame[columns].copy() for column in columns[2:]: output[f"{column}_pct"] = output[column] * 100 return output def build_event_study(frame: pd.DataFrame) -> pd.DataFrame: event_indices = frame.index[frame["tail_event_start"] == 1].tolist() rows: list[dict[str, object]] = [] for level, spec in OUTCOME_SPECS.items(): log_level = np.log(frame[level]) for horizon in range(-6, 13): values: list[float] = [] for event_index in event_indices: baseline_index = event_index - 1 target_index = event_index + horizon if baseline_index < 0 or target_index < 0 or target_index >= len(frame): continue values.append(float(log_level.iloc[target_index] - log_level.iloc[baseline_index])) array = np.asarray(values, dtype=float) mean = float(np.mean(array)) median = float(np.median(array)) se = float(np.std(array, ddof=1) / np.sqrt(len(array))) if len(array) > 1 else np.nan rows.append( { "outcome": level, "outcome_label": spec["label"], "horizon_month": horizon, "events": len(array), "mean_cumulative_log_change": mean, "median_cumulative_log_change": median, "standard_error": se, "ci_low": mean - 1.96 * se if np.isfinite(se) else np.nan, "ci_high": mean + 1.96 * se if np.isfinite(se) else np.nan, } ) return pd.DataFrame(rows) def build_local_projections(frame: pd.DataFrame) -> pd.DataFrame: rows: list[dict[str, object]] = [] work = frame.copy() work["taiex_log_return_lag1"] = work["taiex_log_return"].shift(1) channel_columns = sorted( { channel for spec in OUTCOME_SPECS.values() for channel in spec["channels"] } ) for channel in channel_columns: work[f"{channel}_lag1"] = work[channel].shift(1) for level, spec in OUTCOME_SPECS.items(): change = spec["change"] work[f"{change}_lag1"] = work[change].shift(1) log_level = np.log(work[level]) for horizon in range(0, 13): target = f"lp_{level}_h{horizon}" work[target] = log_level.shift(-horizon) - log_level.shift(1) base_regressors = [ "stock_tail_loss", "taiex_log_return_lag1", f"{change}_lag1", ] for model_name, channels in [ ("base", []), ("channel_conditioned", spec["channels"]), ]: regressors = list(base_regressors) for channel in channels: regressors.extend([channel, f"{channel}_lag1"]) model, selected = fit_hac(work, target, regressors) rows.append( { "outcome": level, "outcome_label": spec["label"], "horizon_month": horizon, "model": model_name, "tail_loss_coefficient": float(model.params["stock_tail_loss"]), "standard_error": float(model.bse["stock_tail_loss"]), "p_value": float(model.pvalues["stock_tail_loss"]), "ci_low": float(model.conf_int().loc["stock_tail_loss", 0]), "ci_high": float(model.conf_int().loc["stock_tail_loss", 1]), "observations": int(len(selected)), "controls": "；".join(channels) if channels else "落後股市報酬、落後依變數與月份固定效果", } ) return pd.DataFrame(rows) def build_event_definition_robustness(frame: pd.DataFrame) -> pd.DataFrame: """Compare conclusions across transparent crash definitions. These regressions are sensitivity checks, not additional causal identification. The leave-GFC-out definition sets tail-loss months from 2008-06 through 2008-11 to zero because that cluster otherwise contributes five of the fifteen primary tail observations. """ work = frame.copy() work["fixed_loss_10pct"] = (work["taiex_log_return"] <= -0.10).astype(int) work["tail_event_start"] = event_cluster_starts(work) work["tail_loss_leave_gfc_out"] = work["stock_tail_loss"].copy() gfc_mask = work["date"].between(pd.Timestamp("2008-06-01"), pd.Timestamp("2008-11-01")) work.loc[gfc_mask, "tail_loss_leave_gfc_out"] = 0 work["taiex_log_return_lag1"] = work["taiex_log_return"].shift(1) definitions = { "empirical_bottom_5pct": "stock_tail_loss", "clustered_event_start": "tail_event_start", "fixed_monthly_loss_10pct": "fixed_loss_10pct", "rolling_24m_drawdown_onset_20pct": "drawdown_onset_20pct", "empirical_bottom_5pct_leave_gfc_out": "tail_loss_leave_gfc_out", } rows: list[dict[str, object]] = [] for level, spec in OUTCOME_SPECS.items(): change = spec["change"] work[f"{change}_lag1"] = work[change].shift(1) for channel in spec["channels"]: if f"{channel}_lag1" not in work: work[f"{channel}_lag1"] = work[channel].shift(1) log_level = np.log(work[level]) for horizon in [0, 3, 6, 12]: target = f"robust_{level}_h{horizon}" work[target] = log_level.shift(-horizon) - log_level.shift(1) for definition, event_variable in definitions.items(): regressors = [ event_variable, "taiex_log_return_lag1", f"{change}_lag1", ] for channel in spec["channels"]: regressors.extend([channel, f"{channel}_lag1"]) model, selected = fit_hac(work, target, regressors) rows.append( { "outcome": level, "outcome_label": spec["label"], "horizon_month": horizon, "definition": definition, "event_variable": event_variable, "event_observations_full_sample": int(work[event_variable].sum()), "coefficient": float(model.params[event_variable]), "standard_error": float(model.bse[event_variable]), "p_value": float(model.pvalues[event_variable]), "ci_low": float(model.conf_int().loc[event_variable, 0]), "ci_high": float(model.conf_int().loc[event_variable, 1]), "observations": int(len(selected)), } ) return pd.DataFrame(rows) def build_asymmetric_distributed_lags(frame: pd.DataFrame) -> pd.DataFrame: rows: list[dict[str, object]] = [] work = frame.copy() for lag in range(0, 7): work[f"stock_loss_lag{lag}"] = work["stock_loss_magnitude"].shift(lag) work[f"stock_gain_lag{lag}"] = work["stock_gain_magnitude"].shift(lag) for level, spec in OUTCOME_SPECS.items(): change = spec["change"] work[f"{change}_lag1"] = work[change].shift(1) regressors = [ *[f"stock_loss_lag{lag}" for lag in range(0, 7)], *[f"stock_gain_lag{lag}" for lag in range(0, 7)], f"{change}_lag1", ] for channel in spec["channels"]: if f"{channel}_lag1" not in work: work[f"{channel}_lag1"] = work[channel].shift(1) regressors.extend([channel, f"{channel}_lag1"]) model, selected = fit_hac(work, change, regressors) loss_names = [f"stock_loss_lag{lag}" for lag in range(0, 7)] gain_names = [f"stock_gain_lag{lag}" for lag in range(0, 7)] loss = linear_combination(model, loss_names, scale=0.10) gain = linear_combination(model, gain_names, scale=0.10) # loss_magnitude = max(-r, 0), while gain_magnitude = max(r, 0). # Under sign-symmetric price responses, the two cumulative # coefficients have equal magnitude and opposite sign. Therefore the # correct null is loss + gain = 0, not loss - gain = 0. difference_vector_names = loss_names + gain_names difference_weights = [0.10] * len(loss_names) + [0.10] * len(gain_names) vector = np.zeros(len(model.params)) positions = {name: index for index, name in enumerate(model.params.index)} for name, weight in zip(difference_vector_names, difference_weights): vector[positions[name]] = weight difference_test = model.t_test(vector) difference_ci = np.asarray(difference_test.conf_int()).reshape(-1) for direction, result in [("10%股市下跌", loss), ("10%股市上漲", gain)]: rows.append( { "outcome": level, "outcome_label": spec["label"], "direction": direction, "cumulative_0_6m_log_effect": result["estimate"], "cumulative_0_6m_pct_approx": result["estimate"] * 100, "standard_error": result["standard_error"], "p_value": result["p_value"], "ci_low": result["ci_low"], "ci_high": result["ci_high"], "observations": int(len(selected)), "loss_plus_gain_log_effect": float( np.asarray(difference_test.effect).reshape(-1)[0] ), "sign_symmetry_p_value": float( np.asarray(difference_test.pvalue).reshape(-1)[0] ), "loss_plus_gain_ci_low": float(difference_ci[0]), "loss_plus_gain_ci_high": float(difference_ci[1]), } ) return pd.DataFrame(rows) def stationarity_table(frame: pd.DataFrame) -> pd.DataFrame: levels = { "TAIEX月平均": "taiex_monthly_average", "NTD/USD": "ntd_per_usd", "FAO食品": "fao_food_index", "FAO穀物": "fao_cereals_index", "進口穀類（美元）": "import_cereals_usd", "進口穀類（新臺幣）": "import_cereals_twd", "臺灣食品CPI": "cpi_food", "臺灣穀類製品CPI": "cpi_grain_products", "臺灣米類製品CPI": "cpi_rice_products", "Brent": "brent_usd_bbl", } rows: list[dict[str, object]] = [] for label, column in levels.items(): log_level = np.log(frame[column]).dropna() for transform, series in [ ("對數水準", log_level), ("一階對數差分", log_level.diff().dropna()), ]: statistic, p_value, used_lag, observations, critical, _ = adfuller( series, regression="c", autolag="AIC", ) rows.append( { "series": label, "transform": transform, "adf_statistic": statistic, "p_value": p_value, "used_lag": used_lag, "observations": observations, "critical_5pct": critical["5%"], "reject_unit_root_at_5pct": p_value < 0.05, } ) return pd.DataFrame(rows) def descriptive_statistics(frame: pd.DataFrame) -> pd.DataFrame: changes = { "TAIEX月報酬": "taiex_log_return", "NTD/USD變動": "fx_depreciation", "Brent變動": "dlog_brent_usd_bbl", "FAO食品變動": "dlog_fao_food_index", "FAO穀物變動": "dlog_fao_cereals_index", "進口穀類美元變動": "dlog_import_cereals_usd", "進口穀類新臺幣變動": "dlog_import_cereals_twd", "臺灣食品CPI變動": "dlog_cpi_food", "臺灣穀類製品CPI變動": "dlog_cpi_grain_products", "臺灣米類製品CPI變動": "dlog_cpi_rice_products", } rows: list[dict[str, object]] = [] for label, column in changes.items(): values = frame[column].dropna() * 100 rows.append( { "series": label, "observations": len(values), "mean_pct": values.mean(), "std_pct": values.std(), "minimum_pct": values.min(), "median_pct": values.median(), "maximum_pct": values.max(), } ) return pd.DataFrame(rows) def correlation_table(frame: pd.DataFrame) -> pd.DataFrame: mapping = { "TAIEX": "taiex_log_return", "NTD/USD": "fx_depreciation", "Brent": "dlog_brent_usd_bbl", "FAO穀物": "dlog_fao_cereals_index", "進口穀物USD": "dlog_import_cereals_usd", "進口穀物TWD": "dlog_import_cereals_twd", "食品CPI": "dlog_cpi_food", "穀類CPI": "dlog_cpi_grain_products", "米類CPI": "dlog_cpi_rice_products", } correlation = frame[list(mapping.values())].corr() correlation.index = mapping correlation.columns = mapping return correlation.reset_index(names="series") def plot_indexed_series(frame: pd.DataFrame) -> None: fig, axes = plt.subplots( 2, 1, figsize=(12.4, 7.8), sharex=True, gridspec_kw={"height_ratios": [1.45, 1]}, ) upper_series = [ ("taiex_monthly_average", "TAIEX月平均", COLORS["brown"]), ("fao_cereals_index", "FAO國際穀物", COLORS["orange"]), ("import_cereals_twd", "進口穀類（新臺幣）", COLORS["blue"]), ] lower_series = [ ("cpi_food", "臺灣食品CPI", COLORS["green"]), ("cpi_grain_products", "臺灣穀類製品CPI", COLORS["gold"]), ("cpi_rice_products", "臺灣米類製品CPI", COLORS["purple"]), ] for ax, series in zip(axes, [upper_series, lower_series]): for column, label, color in series: indexed = frame[column] / frame[column].iloc[0] * 100 ax.plot(frame["date"], indexed, label=label, color=color, linewidth=1.7) for date in frame.loc[frame["tail_event_start"] == 1, "date"]: ax.axvline(date, color=COLORS["red"], alpha=0.12, linewidth=1) ax.legend(ncol=3, frameon=False, loc="upper left") ax.grid(axis="y", alpha=0.2) ax.spines[["top", "right"]].set_visible(False) ax.set_ylabel("2002-01=100") axes[0].set_title("圖1　股市、國際／進口穀價與臺灣零售物價的分層路徑") axes[0].text( 0.995, 0.06, "垂直線為10個尾端事件起點", transform=axes[0].transAxes, ha="right", color=COLORS["muted"], fontsize=8, ) fig.tight_layout() fig.savefig(FIGURES / "figure-1-indexed-price-paths.png", bbox_inches="tight") plt.close(fig) def plot_tail_events(frame: pd.DataFrame, events: pd.DataFrame) -> None: fig, ax = plt.subplots(figsize=(12.4, 6.2)) colors = np.where(frame["stock_tail_loss"] == 1, COLORS["red"], COLORS["blue"]) ax.bar(frame["date"], frame["taiex_log_return"] * 100, color=colors, width=25, alpha=0.82) for event_number, (_, event) in enumerate(events.iterrows(), start=1): date = pd.Timestamp(event["date"]) y = event["taiex_log_return_pct"] ax.annotate( str(event_number), xy=(date, y), xytext=(0, -15 if y < 0 else 8), textcoords="offset points", ha="center", va="top" if y < 0 else "bottom", fontsize=8, color=COLORS["ink"], ) ax.axhline(0, color=COLORS["ink"], linewidth=0.8) ax.set_title("圖2　TAIEX月報酬與最差5%尾端月份") ax.set_ylabel("月對數報酬（%）") ax.spines[["top", "right"]].set_visible(False) ax.grid(axis="y", alpha=0.18) fig.tight_layout() fig.savefig(FIGURES / "figure-2-tail-loss-events.png", bbox_inches="tight") plt.close(fig) def plot_event_study(event_study: pd.DataFrame) -> None: fig, ax = plt.subplots(figsize=(11.2, 6.8)) selected = { "fao_cereals_index": COLORS["orange"], "import_cereals_twd": COLORS["blue"], "cpi_food": COLORS["green"], "cpi_rice_products": COLORS["purple"], } for outcome, color in selected.items(): subset = event_study.loc[event_study["outcome"] == outcome] x = subset["horizon_month"].to_numpy(dtype=float) mean = subset["mean_cumulative_log_change"].to_numpy(dtype=float) * 100 low = subset["ci_low"].to_numpy(dtype=float) * 100 high = subset["ci_high"].to_numpy(dtype=float) * 100 label = subset["outcome_label"].iloc[0] ax.plot(x, mean, marker="o", markersize=3.5, linewidth=1.7, color=color, label=label) ax.fill_between(x, low, high, color=color, alpha=0.10) ax.axvline(-1, color=COLORS["ink"], linestyle="--", linewidth=1) ax.axhline(0, color=COLORS["ink"], linewidth=0.8) ax.set_title("圖3　尾端股市事件前後的平均累積價格變化（描述性）") ax.set_xlabel("事件月相對月份（-1為共同基準）") ax.set_ylabel("相對事件前一月的累積對數變化（%）") ax.legend(frameon=False, ncol=2) ax.grid(alpha=0.18) ax.spines[["top", "right"]].set_visible(False) fig.tight_layout() fig.savefig(FIGURES / "figure-3-event-study-descriptive.png", bbox_inches="tight") plt.close(fig) def plot_local_projections(local: pd.DataFrame) -> None: fig, axes = plt.subplots(1, 3, figsize=(14.4, 5.3), sharey=False) selected = [ ("fao_cereals_index", "FAO國際穀物", COLORS["orange"]), ("import_cereals_twd", "進口穀類（新臺幣）", COLORS["blue"]), ("cpi_food", "臺灣食品CPI", COLORS["green"]), ] for ax, (outcome, title, color) in zip(axes, selected): subset = local.loc[ (local["outcome"] == outcome) & (local["model"] == "channel_conditioned") ] x = subset["horizon_month"].to_numpy(dtype=float) estimate = subset["tail_loss_coefficient"].to_numpy(dtype=float) * 100 low = subset["ci_low"].to_numpy(dtype=float) * 100 high = subset["ci_high"].to_numpy(dtype=float) * 100 ax.plot(x, estimate, color=color, marker="o", markersize=3.5) ax.fill_between(x, low, high, color=color, alpha=0.16) ax.axhline(0, color=COLORS["ink"], linewidth=0.8) ax.set_title(title) ax.set_xlabel("事件後月數") ax.grid(alpha=0.18) ax.spines[["top", "right"]].set_visible(False) axes[0].set_ylabel("尾端股市月份的條件關聯（近似%）") fig.suptitle("圖4　加入匯率、能源與國際價格後的局部投影關聯", y=1.02) fig.tight_layout() fig.savefig(FIGURES / "figure-4-local-projections.png", bbox_inches="tight") plt.close(fig) def plot_asymmetry(asymmetric: pd.DataFrame) -> None: labels = [spec["label"] for spec in OUTCOME_SPECS.values()] losses = [] gains = [] loss_errors_low = [] loss_errors_high = [] gain_errors_low = [] gain_errors_high = [] for outcome in OUTCOME_SPECS: subset = asymmetric.loc[asymmetric["outcome"] == outcome] loss_row = subset.loc[subset["direction"] == "10%股市下跌"].iloc[0] gain_row = subset.loc[subset["direction"] == "10%股市上漲"].iloc[0] losses.append(loss_row["cumulative_0_6m_pct_approx"]) gains.append(gain_row["cumulative_0_6m_pct_approx"]) loss_errors_low.append((loss_row["cumulative_0_6m_log_effect"] - loss_row["ci_low"]) * 100) loss_errors_high.append((loss_row["ci_high"] - loss_row["cumulative_0_6m_log_effect"]) * 100) gain_errors_low.append((gain_row["cumulative_0_6m_log_effect"] - gain_row["ci_low"]) * 100) gain_errors_high.append((gain_row["ci_high"] - gain_row["cumulative_0_6m_log_effect"]) * 100) y = np.arange(len(labels)) fig, ax = plt.subplots(figsize=(10.8, 6.4)) width = 0.36 ax.barh( y - width / 2, losses, height=width, xerr=np.asarray([loss_errors_low, loss_errors_high]), color=COLORS["red"], alpha=0.85, capsize=3, label="10%股市下跌", ) ax.barh( y + width / 2, gains, height=width, xerr=np.asarray([gain_errors_low, gain_errors_high]), color=COLORS["blue"], alpha=0.85, capsize=3, label="10%股市上漲", ) ax.axvline(0, color=COLORS["ink"], linewidth=0.8) ax.set_yticks(y, labels) ax.invert_yaxis() ax.set_xlabel("0-6個月累積關聯（近似%）") ax.set_title("圖5　股市下跌與上漲的非對稱分散落後估計") ax.legend(frameon=False) ax.grid(axis="x", alpha=0.18) ax.spines[["top", "right"]].set_visible(False) fig.tight_layout() fig.savefig(FIGURES / "figure-5-asymmetric-distributed-lags.png", bbox_inches="tight") plt.close(fig) def plot_event_heatmap(events: pd.DataFrame) -> None: variables = [ ("taiex_log_return_pct", "TAIEX"), ("fx_depreciation_pct", "NTD/USD"), ("dlog_brent_usd_bbl_pct", "Brent"), ("dlog_fao_cereals_index_pct", "FAO穀物"), ("dlog_import_cereals_twd_pct", "進口穀物TWD"), ("dlog_cpi_food_pct", "食品CPI"), ] matrix = events[[column for column, _ in variables]].to_numpy(dtype=float) scale = np.nanstd(matrix, axis=0, ddof=1) standardized = np.divide(matrix, scale, out=np.zeros_like(matrix), where=scale > 0) standardized = np.clip(standardized, -3, 3) fig, ax = plt.subplots(figsize=(10.8, 6.8)) image = ax.imshow(standardized, cmap="RdBu_r", vmin=-3, vmax=3, aspect="auto") ax.set_xticks(range(len(variables)), [label for _, label in variables], rotation=20, ha="right") ax.set_yticks( range(len(events)), [ f"{pd.Timestamp(row.date).strftime('%Y-%m')}｜{row.event_type}" for row in events.itertuples() ], ) for i in range(matrix.shape[0]): for j in range(matrix.shape[1]): ax.text(j, i, f"{matrix[i, j]:.1f}", ha="center", va="center", fontsize=8) ax.set_title("圖6　尾端事件月的跨市場變化異質性（格內為月變動%，顏色為標準化值）") fig.colorbar(image, ax=ax, fraction=0.035, pad=0.03, label="欄內標準差") fig.tight_layout() fig.savefig(FIGURES / "figure-6-event-channel-heatmap.png", bbox_inches="tight") plt.close(fig) def finite_or_none(value: float) -> float | None: return float(value) if math.isfinite(float(value)) else None def main() -> None: frame = pd.read_csv(DATA / "monthly-analysis-panel.csv") frame["date"] = pd.to_datetime(frame["date"], format="%Y-%m") frame["tail_event_start"] = event_cluster_starts(frame) descriptive = descriptive_statistics(frame) stationarity = stationarity_table(frame) correlations = correlation_table(frame) events = descriptive_event_table(frame) event_study = build_event_study(frame) local = build_local_projections(frame) local["fdr_q_value"] = np.nan for model_name, indices in local.groupby("model").groups.items(): local.loc[indices, "fdr_q_value"] = multipletests( local.loc[indices, "p_value"].to_numpy(dtype=float), method="fdr_bh", )[1] robustness = build_event_definition_robustness(frame) asymmetric = build_asymmetric_distributed_lags(frame) write_csv(descriptive, "descriptive-statistics.csv") write_csv(stationarity, "stationarity-tests.csv") write_csv(correlations, "correlation-matrix.csv") write_csv(events, "tail-event-months.csv") write_csv(event_study, "event-study-descriptive.csv") write_csv(local, "local-projection-results.csv") write_csv(robustness, "event-definition-robustness.csv") write_csv(asymmetric, "asymmetric-distributed-lag-results.csv") attenuation = ( local.loc[local["horizon_month"] == 6] .pivot( index=["outcome", "outcome_label"], columns="model", values="tail_loss_coefficient", ) .reset_index() ) attenuation["coefficient_change_after_channels"] = ( attenuation["channel_conditioned"] - attenuation["base"] ) write_csv(attenuation, "channel-conditioning-comparison-h6.csv") plot_indexed_series(frame) plot_tail_events(frame, events) plot_event_study(event_study) plot_local_projections(local) plot_asymmetry(asymmetric) plot_event_heatmap(events) h6 = local.loc[ (local["horizon_month"] == 6) & (local["model"] == "channel_conditioned") ].set_index("outcome") asym_loss = asymmetric.loc[ asymmetric["direction"] == "10%股市下跌" ].set_index("outcome") event_counts = events["event_type"].value_counts().to_dict() results = { "report_number": "SHWRP-2026-028", "analysis_date": "2026-07-30", "sample_start": frame["date"].min().strftime("%Y-%m"), "sample_end": frame["date"].max().strftime("%Y-%m"), "months": int(len(frame)), "tail_threshold_log_return": float(frame["taiex_log_return"].quantile(0.05)), "tail_loss_months": int(frame["stock_tail_loss"].sum()), "clustered_tail_events": int(frame["tail_event_start"].sum()), "rolling_24m_drawdown_onsets": int(frame["drawdown_onset_20pct"].sum()), "event_definition_counts": { str(row["definition"]): int(row["event_observations_full_sample"]) for _, row in ( robustness.loc[robustness["outcome"] == "cpi_food"] .drop_duplicates("definition") .iterrows() ) }, "event_type_counts": {str(key): int(value) for key, value in event_counts.items()}, "h6_channel_conditioned_tail_coefficients": { outcome: { "estimate_log_points": finite_or_none(row["tail_loss_coefficient"]), "p_value": finite_or_none(row["p_value"]), "ci_low": finite_or_none(row["ci_low"]), "ci_high": finite_or_none(row["ci_high"]), } for outcome, row in h6.iterrows() }, "asymmetric_0_6m_effect_for_10pct_stock_loss": { outcome: { "estimate_log_points": finite_or_none(row["cumulative_0_6m_log_effect"]), "approx_percent": finite_or_none(row["cumulative_0_6m_pct_approx"]), "p_value": finite_or_none(row["p_value"]), "sign_symmetry_p_value": finite_or_none(row["sign_symmetry_p_value"]), } for outcome, row in asym_loss.iterrows() }, "figure_count": 6, "channel_conditioned_lp_tests": int( (local["model"] == "channel_conditioned").sum() ), "channel_conditioned_lp_nominal_p_below_0_05": int( ( (local["model"] == "channel_conditioned") & (local["p_value"] < 0.05) ).sum() ), "channel_conditioned_lp_fdr_q_below_0_05": int( ( (local["model"] == "channel_conditioned") & (local["fdr_q_value"] < 0.05) ).sum() ), "causal_identification": False, "interpretation": ( "All event-study, local-projection, and distributed-lag estimates are " "time-series associations. The design does not identify a structural " "causal effect of stock-market crashes on food prices." ), "main_claim_boundary": ( "A stock-market tail month has no fixed theoretical sign for food prices. " "Direction and timing must be read jointly with exchange-rate, energy, " "international food-price, and domestic retail-price layers." ), } (ROOT / "analysis-results.json").write_text( json.dumps(results, ensure_ascii=False, indent=2), encoding="utf-8", ) print(json.dumps(results, ensure_ascii=False, indent=2)) if __name__ == "__main__": main()