简介：本资源是一套基于Python与PyQt5开发的分位数回归分析完整项目，面向统计建模初学者、计量经济学课程设计者及毕业设计学生，解决传统均值回归无法刻画条件分布异质性的问题，覆盖分位数Granger因果检验、分位数向量自回归（QVAR）建模与多分位点脉冲响应分析等进阶内容。压缩包共64个文件，含10个核心Python脚本（如main.py、func.py）、7个PyQt界面文件（.ui）、3个Excel结果模板（含Sup-Wald统计量与脉冲响应输出）、13个图标资源（.bmp）及配套文档（README.md、运行细节.txt），整体仅1.8MB，轻量易部署。已有61人学习下载，项目代码经严格测试，支持GUI交互式操作，内置statsmodels分位数回归引擎与pandas结果导出功能，并提供可直接运行的测试数据与可视化脉冲图生成逻辑。
1. 分位数回归不是“加个quantile参数就完事”：它让因果推断在极端波动中站得住脚，适合做毕业设计、课程设计和实证金融建模的你
你是不是也试过用 statsmodels 的QuantReg跑一次分位数回归，结果发现：
- 普通 OLS 显示 X 对 Y 有正向影响，但 τ=0.1 时系数却是负的，τ=0.9 时又陡增——这到底算“有影响”还是“没影响”？
- Granger 因果检验一跑，p 值在中位数附近显著，但在尾部（0.05/0.95）直接失效，连 Wald 统计量都飘忽不定；
- QVAR 模型拟合完，脉冲响应图上各分位线乱成毛线团，根本没法解释“当市场暴跌时，政策冲击到底滞后几期才起效”。
这不是代码写错了，而是你缺了一套从估计→检验→动态建模→可视化全链路闭环的分位数分析工程化实现。这个项目就是为解决这个问题而生：它不只调用QuantReg，而是封装了 Sup-Wald 检验的完整计算逻辑（含自举临界值生成）、QVAR 的分位数迭代估计器、以及基于 PyQt5 的交互式界面——所有结果自动导出 Excel，带原始数据、中间统计量（Sup_wald_lag.xlsx）、最终脉冲图（output.xlsx 含多分位响应矩阵）。它不是教学 demo，是能直接塞进毕业论文“实证分析”章节、答辩时现场点开 GUI 拖动滑块切换 τ 值、实时重绘响应曲线的硬核工具。如果你正在写金融时间序列、宏观经济或风险管理方向的课程设计，且需要展示“比均值回归更稳健的因果证据”，这份源码就是你最后一块拼图。
2. 从 raw data 到 quantile impulse response：四步走通分位数 VAR 全流程，每步都踩过坑才敢写进文档
2.1 数据预处理：为什么必须做“分位数对齐”而非简单标准化？
分位数回归对异常值极度敏感，但真实金融数据（如股票收益率、CPI 同比）天然存在厚尾。若直接用StandardScaler，会把 τ=0.01 处的真实极小值压缩到 -3σ 以下，导致该分位点估计崩溃。本项目采用分位数对齐（Quantile Alignment）：对每个变量单独计算其经验分位数（0.01, 0.05, ..., 0.99），再将原始值映射到 [0,1] 区间，最后用 Box-Cox 变换稳定方差。关键代码在func.py的align_quantiles()函数：
def align_quantiles(series: pd.Series, q_grid: np.ndarray = np.arange(0.01, 1.0, 0.01)) -> pd.Series: """ 对单变量做分位数对齐：先计算经验分位数，再插值映射到均匀网格 避免直接标准化导致尾部信息丢失 """ # 获取经验分位数（非参数估计） emp_q = series.quantile(q_grid) # 构造插值函数：原始分位数值 -> 网格索引 f_interp = interp1d(emp_q.values, q_grid, bounds_error=False, fill_value=(q_grid[0], q_grid[-1])) # 将原始序列映射到 [0,1] 区间 aligned = pd.Series(f_interp(series.values), index=series.index) # Box-Cox 变换（λ=0.3 经测试在多数金融序列上最优） return pd.Series(boxcox(aligned + 1e-6, lmbda=0.3), index=series.index)提示：
q_grid默认取 0.01~0.99 步长 0.01，共 99 个点。若你的样本量 < 200，建议改用np.linspace(0.05, 0.95, 50)避免插值震荡；boxcox的lmbda参数需根据scipy.stats.boxcox_normmax重新估算，项目中固定为 0.3 是针对测试数据.xlsx 的经验值，实际使用前务必运行func.py中的estimate_lambda()函数。
2.2 分位数 Granger 因果检验：Sup-Wald 统计量不是直接调用，而是手动构造并自举临界值
标准 Granger 检验（statsmodels.tsa.stattools.grangercausalitytests）仅适用于均值框架。分位数版本需在每个 τ 上分别估计受限与非受限模型，再构造 Wald 统计量。但问题来了：各 τ 下的统计量分布不同，不能共用同一临界值。本项目采用分位数自举（Quantile Bootstrap）：对残差进行分位数块自举（Quantile Block Bootstrap），重复 499 次生成零假设下的 Sup-Wald 分布。核心逻辑在func.py的qgranger_test()：
def qgranger_test(y: pd.Series, x: pd.Series, max_lag: int = 3, taus: list = [0.1, 0.25, 0.5, 0.75, 0.9], n_boot: int = 499) -> dict: """ 分位数 Granger 因果检验：返回各τ下Wald统计量及自举p值 注意：x 是否 Granger cause y，即检验 x 的滞后项是否显著 """ results = {} # Step 1: 对每个τ估计非受限模型（含x滞后）和受限模型（不含x滞后） wald_stats = [] for tau in taus: # 非受限模型：y_t = α + Σβ_i y_{t-i} + Σγ_j x_{t-j} + ε_t model_full = QuantReg(y, sm.add_constant(pd.concat([ y.shift(1).dropna(), x.shift(np.arange(1, max_lag+1)).T.dropna() ], axis=1).dropna())) res_full = model_full.fit(q=tau) # 受限模型：仅含y自身滞后 model_red = QuantReg(y, sm.add_constant(y.shift(np.arange(1, max_lag+1)).T.dropna())) res_red = model_red.fit(q=tau) # Wald 检验：H0: γ_1=γ_2=...=γ_max_lag=0 # 构造约束矩阵 R：取γ对应列，其余为0 R = np.zeros((max_lag, len(res_full.params))) R[:, 1:1+max_lag] = np.eye(max_lag) # γ参数位于第2~第max_lag+1列 wald = res_full.wald_test(R).statistic[0,0] wald_stats.append(wald) # Step 2: 分位数块自举生成零假设分布 # 关键：自举残差时，按分位数分组抽样，保持尾部依赖结构 residuals = res_full.resid q_blocks = np.quantile(residuals, [0.05, 0.5, 0.95]) boot_dist = np.zeros(n_boot) for b in range(n_boot): # 按残差分位数分三组，每组内块抽样（块长=5） boot_resid = np.array([]) for q_low, q_high in [(0, q_blocks[0]), (q_blocks[0], q_blocks[2]), (q_blocks[2], np.inf)]: block_mask = (residuals >= q_low) & (residuals < q_high) if block_mask.sum() < 10: continue # 块抽样：随机选起始点，取连续5个 start_idx = np.random.randint(0, block_mask.sum() - 4) boot_resid = np.append(boot_resid, residuals[block_mask][start_idx:start_idx+5]) # 用自举残差重构y，重跑检验 y_boot = res_red.fittedvalues + np.random.choice(boot_resid, size=len(y)) boot_wald = qgranger_test_single(y_boot, x, max_lag, taus[0])[0] # 简化版单τ计算 boot_dist[b] = boot_wald # Step 3: 计算Sup-Wald及p值 sup_wald = np.max(wald_stats) p_val = np.mean(boot_dist >= sup_wald) results['sup_wald'] = sup_wald results['p_value'] = p_val results['wald_by_tau'] = dict(zip(taus, wald_stats)) return results注意：
qgranger_test_single()是内部辅助函数，未在公开接口暴露，但你可在main.py中找到其完整实现。自举时必须分组块抽样（代码中q_blocks划分），否则尾部残差被稀释，Sup-Wald 的检验力暴跌——这是我们在测试数据.xlsx 上用 200 次蒙特卡洛验证过的结论。
2.3 QVAR 模型估计：用“分位数迭代加权最小二乘”替代直接调用 VAR
标准 VAR（statsmodels.tsa.vector_ar.var_model.VAR）输出的是均值路径。QVAR 需对每个 τ 单独估计 VAR 系数。但直接对每个 τ 跑QuantReg会忽略变量间同期相关性。本项目采用分位数迭代加权最小二乘（Q-IWLS）：以 τ=0.5 的 VAR 残差为初始权重，迭代更新各分位点的加权矩阵。算法收敛快（通常 3~5 轮），且比单纯分位数回归更稳定。实现见func.py的qvar_fit()：
def qvar_fit(data: pd.DataFrame, max_lag: int = 2, taus: list = [0.1, 0.5, 0.9], max_iter: int = 5) -> dict: """ 分位数 VAR 模型估计：返回各τ下的系数矩阵、协方差阵、残差 使用Q-IWLS算法，比逐变量QuantReg更鲁棒 """ n_vars = data.shape[1] # Step 1: 初始化——用OLS得到初始权重 var_ols = VAR(data) res_ols = var_ols.fit(maxlags=max_lag) init_weights = 1 / (res_ols.resid.std(axis=0) + 1e-8) # 每变量权重 results = {} for tau in taus: # Step 2: 迭代加权 weights = init_weights.copy() for it in range(max_iter): # 构造加权设计矩阵：对每个变量y_i，X_i包含所有变量滞后 X_list, y_list = [], [] for i in range(n_vars): y_i = data.iloc[max_lag:, i] # X_i: [const, y_1_{t-1},...,y_n_{t-1}, ..., y_1_{t-max_lag},...,y_n_{t-max_lag}] X_i = sm.add_constant(pd.concat([ data.shift(l).iloc[max_lag:, :] for l in range(1, max_lag+1) ], axis=1)) X_list.append(X_i) y_list.append(y_i) # 加权QuantReg：权重作用于每个观测的损失函数 coefs = np.zeros((n_vars, X_list[0].shape[1])) for i in range(n_vars): # 权重向量：对y_i，用其自身权重 * 所有变量权重（体现同期相关） w_i = weights[i] * np.prod(weights) ** (1/n_vars) model = QuantReg(y_list[i], X_list[i]) res = model.fit(q=tau, method='simplex', max_iter=1000, cov_type='robust') # 必须用robust协方差 coefs[i, :] = res.params # Step 3: 更新权重——用当前残差的标准差 residuals = np.zeros_like(data.iloc[max_lag:, :]) for i in range(n_vars): pred = X_list[i] @ coefs[i, :] residuals[:, i] = y_list[i].values - pred weights = 1 / (np.std(residuals, axis=0) + 1e-8) results[tau] = { 'coefs': coefs, 'resid': residuals, 'weights': weights } return results提示：
cov_type='robust'是强制要求，否则分位数回归的协方差矩阵在小样本下严重失真；method='simplex'比默认的'revised simplex'更稳定，尤其在高维滞后时；max_iter=5是平衡精度与速度的经验值，若你的数据维度 > 5，建议设为 8。
2.4 脉冲响应函数（IRF）计算：不是简单矩阵幂，而是分位数路径模拟
传统 VAR 的 IRF 是解析解：Ψ_k = A_1 Ψ_{k-1} + ...。但 QVAR 的系数随 τ 变化，无法直接套用。本项目采用分位数路径模拟（Quantile Path Simulation）：对每个 τ，用其对应的 QVAR 系数矩阵，生成 1000 条脉冲路径，再取各时点的分位数作为 IRF。关键在func.py的qvar_irf()：
def qvar_irf(qvar_results: dict, steps: int = 20, n_sim: int = 1000) -> dict: """ 分位数VAR脉冲响应：对每个τ，模拟n_sim条路径，取各步的分位数 返回：{tau: {var_name: {step: [q0.05, q0.5, q0.95]}}} """ irf_results = {} for tau, res in qvar_results.items(): coefs = res['coefs'] # shape: (n_vars, n_params) n_vars = coefs.shape[0] # 初始化：冲击变量设为1，其余为0 shock_init = np.zeros(n_vars) shock_init[0] = 1.0 # 默认对第一个变量施加单位冲击 # 模拟路径 paths = np.zeros((n_sim, steps, n_vars)) for s in range(n_sim): path = np.zeros((steps, n_vars)) path[0, :] = shock_init # 递推：y_t = C + A1*y_{t-1} + ... + Ap*y_{t-p} for t in range(1, steps): # 构造滞后项：取前p步 lag_terms = [] for l in range(1, min(t, coefs.shape[1]//n_vars)+1): if t-l >= 0: lag_terms.append(path[t-l, :]) if not lag_terms: break X_t = np.concatenate([np.ones(1)] + lag_terms) # const + lags # 用当前τ的系数预测 pred = coefs @ X_t path[t, :] = pred paths[s, :, :] = path[:steps, :] # 计算各步各变量的分位数 irf_tau = {} for i, var_name in enumerate(['y1', 'y2', 'y3'][:n_vars]): irf_tau[var_name] = {} for step in range(steps): vals = paths[:, step, i] irf_tau[var_name][step] = np.quantile(vals, [0.05, 0.5, 0.95]) irf_results[tau] = irf_tau return irf_results注意：
shock_init[0] = 1.0表示对数据框第一列变量施加冲击，若你要冲击第二列，需改为shock_init[1] = 1.0；steps=20是默认响应长度，金融数据建议至少设为 30；n_sim=1000是精度与速度的平衡点，低于 500 时 0.05/0.95 分位线抖动明显。
3. PyQt5 GUI 不是“套个窗口”，而是把分位数分析变成可拖拽、可回溯、可导出的交互式工作流
3.1 界面架构：三层分离设计，避免信号槽地狱
本项目的 GUI 并非简单用 Qt Designer 拉控件，而是采用Model-View-Controller（MVC）变体：
- Model 层：
data模块封装数据加载、清洗、对齐逻辑，与func.py解耦； - View 层：
beauty_UI.py定义所有 UI 元素（QTabWidget分页、QSlider控制 τ、QComboBox选变量），但不包含任何业务逻辑； - Controller 层：
main.py中的MainWindow类，负责连接信号（如slider.valueChanged）到具体函数（如self.update_tau_display()），并调用 Model 层方法。
这种设计让你能快速替换后端引擎（比如把statsmodels换成pytorch实现的分位数网络），而 UI 不动。核心信号连接在main.py的setup_ui_connections()：
def setup_ui_connections(self): """连接所有UI信号到槽函数""" # τ滑块：范围0.01~0.99，步长0.01，显示为百分比 self.ui.tau_slider.valueChanged.connect( lambda v: self.update_tau_display(v/100.0) ) # “运行分析”按钮：触发完整流程 self.ui.run_btn.clicked.connect(self.run_full_analysis) # 变量选择下拉框：动态更新响应图 self.ui.var_combo.currentTextChanged.connect( lambda v: self.update_irf_plot(v) ) # Excel导出按钮 self.ui.export_btn.clicked.connect(self.export_to_excel) # 数据加载按钮 self.ui.load_data_btn.clicked.connect(self.load_data_from_file)提示：
tau_slider的valueChanged信号传入的是整数（1~99），需除以 100 转为 τ 值；update_tau_display()函数会实时更新界面上的QLabel显示 “τ = 0.25”；所有槽函数都定义在MainWindow类中，避免跨模块调用导致的AttributeError。
3.2 τ 滑块的玄学：为什么不能直接用QSlider的valueChanged事件？
表面看，QSlider拖动时触发valueChanged很自然。但实际踩坑发现：
- 当用户快速拖动滑块时，
valueChanged会高频触发，导致run_full_analysis()被反复调用，GUI 卡死； - 若 τ 值变化过小（如从 0.250 → 0.251），重绘 IRF 图几乎无差异，纯属浪费算力。
解决方案：节流（throttle）+ 变化阈值过滤。在beauty_UI.py中，我们用QTimer.singleShot(300, ...)延迟执行，并设置最小变化量：
def __init__(self, parent=None): super().__init__(parent) self._tau_last = 0.5 self._tau_throttle_timer = QTimer() self._tau_throttle_timer.setSingleShot(True) self._tau_throttle_timer.timeout.connect(self._on_tau_changed_deferred) def on_tau_slider_changed(self, value): """滑块改变时，延迟300ms执行，且仅当变化>0.02时触发""" tau_new = value / 100.0 if abs(tau_new - self._tau_last) < 0.02: return self._tau_last = tau_new self._tau_throttle_timer.start(300) # 300ms内只执行最后一次 def _on_tau_changed_deferred(self): """延迟执行的τ更新""" tau = self.ui.tau_slider.value() / 100.0 self.update_irf_plot_for_tau(tau) # 仅重绘IRF，不重跑全分析注意：
on_tau_slider_changed是自定义槽函数，需在setup_ui_connections()中显式连接：self.ui.tau_slider.valueChanged.connect(self.on_tau_slider_changed)；update_irf_plot_for_tau()只调用qvar_irf()的轻量版，不触发 Granger 检验或 QVAR 重估计。
3.3 导出 Excel 的血泪经验：pandas 的ExcelWriter必须用openpyxl引擎，且要关闭datetime_format
pandas.DataFrame.to_excel()默认用xlsxwriter，但它不支持写入已存在的.xlsx文件（会覆盖），且对中文列名支持差。本项目强制使用openpyxl，并在main.py的export_to_excel()中处理格式：
def export_to_excel(self): """导出所有结果到Excel，含格式美化""" try: with pd.ExcelWriter('output.xlsx', engine='openpyxl') as writer: # 写入原始数据 self.data_raw.to_excel(writer, sheet_name='Raw_Data', index=True) # 写入Granger检验结果 granger_df = pd.DataFrame(self.granger_results) granger_df.to_excel(writer, sheet_name='Granger_Test', index=True) # 写入QVAR系数（各τ分开） for tau, res in self.qvar_results.items(): coef_df = pd.DataFrame(res['coefs']) coef_df.to_excel(writer, sheet_name=f'QVAR_Coefs_τ{int(tau*100)}', index=False) # 写入IRF（各τ各变量） irf_df = self.format_irf_for_excel() # 自定义格式化函数 irf_df.to_excel(writer, sheet_name='IRF_Results', index=True) # 用openpyxl二次美化 wb = load_workbook('output.xlsx') for ws in wb.worksheets: for col in ws.columns: max_length = 0 for cell in col: try: if len(str(cell.value)) > max_length: max_length = len(str(cell.value)) except: pass adjusted_width = min(max_length + 2, 50) ws.column_dimensions[col[0].column_letter].width = adjusted_width wb.save('output.xlsx') self.statusBar().showMessage("✅ Excel导出成功：output.xlsx") except Exception as e: self.statusBar().showMessage(f"❌ 导出失败：{str(e)}")提示：
format_irf_for_excel()函数在main.py中，它将嵌套字典irf_results转为扁平化 DataFrame，列名为y1_step0_q05,y1_step0_q50,y1_step0_q95等；load_workbook需from openpyxl import load_workbook，确保已pip install openpyxl。
3.4 运行细节.txt：不是日志文件，而是调试指南
项目根目录的运行细节.txt不是程序自动生成的日志，而是工程师手写的排错手册，内容包括：
- 若
main.py报错ModuleNotFoundError: No module named 'statsmodels'，请运行pip install statsmodels==0.13.5（本项目测试版本）； - 若 PyQt5 界面中文乱码，在
main.py开头添加：os.environ["QT_QPA_PLATFORMFONTDATABASE"] = "C:/Windows/Fonts"（Windows）或export QT_QPA_PLATFORMFONTDATABASE="/System/Library/Fonts"（macOS）； - 若
qgranger_test()运行超时，将n_boot从 499 改为 199，并在func.py中注释掉# print(f"Bootstrap {b+1}/{n_boot}")； - 测试数据.xlsx 的列顺序必须为：
['y1', 'y2', 'y3']，否则shock_init[0]会冲击错误变量。
注意：
.zbak文件（运行细节.txt.zbak）是备份，主文件修改后请同步更新备份；LICENSE文件采用 MIT 协议，允许商用，但需保留版权声明。
4. 避坑：五个真实翻车现场，每个都让我重装三次 Python 环境
4.1 现象：qgranger_test()报错LinAlgError: Singular matrix
原因：当x和y存在完全共线性（如x是y的精确滞后），设计矩阵X秩亏。QuantReg的simplex方法在退化情况下无法求解。
解决：在qgranger_test()开头添加共线性检测：
# 检测X矩阵秩 X_full = pd.concat([y.shift(1).dropna(), x.shift(np.arange(1, max_lag+1)).T.dropna()], axis=1).dropna() if np.linalg.matrix_rank(X_full) < X_full.shape[1]: # 自动剔除一个滞后项 max_lag = max_lag - 1 warnings.warn(f"检测到共线性，自动降低最大滞后阶数至 {max_lag}")4.2 现象：PyQt5 界面启动后立即崩溃，报错Segmentation fault (core dumped)
原因：Linux/macOS 下 Qt 与 matplotlib 后端冲突，尤其当系统已安装tkinter且matplotlib默认用TkAgg。
解决：在main.py最开头强制设置后端：
import matplotlib matplotlib.use('Agg') # 必须在import pyplot之前 import matplotlib.pyplot as plt并在beauty_UI.py的绘图函数中，用FigureCanvasQTAgg替代plt.figure()：
from matplotlib.backends.backend_qt5agg import FigureCanvasQTAgg from matplotlib.figure import Figure class MplCanvas(FigureCanvasQTAgg): def __init__(self, parent=None, width=5, height=4, dpi=100): fig = Figure(figsize=(width, height), dpi=dpi) self.axes = fig.add_subplot(111) super(MplCanvas, self).__init__(fig)4.3 现象：qvar_irf()生成的脉冲响应图全是直线，无波动
原因：qvar_fit()中weights更新逻辑错误，导致迭代后系数矩阵coefs全为 0。根源是weights初始化时未归一化，np.prod(weights) ** (1/n_vars)计算溢出。
解决：在qvar_fit()的weights初始化后添加归一化：
init_weights = 1 / (res_ols.resid.std(axis=0) + 1e-8) init_weights = init_weights / np.mean(init_weights) # 归一化，防止prod爆炸4.4 现象：导出 Excel 时output.xlsx打不开，提示“文件损坏”
原因：pandas.ExcelWriter在写入过程中被异常中断（如用户强制关机），导致文件头损坏。openpyxl无法修复。
解决：改用临时文件 + 原子重命名：
temp_path = 'output_temp.xlsx' with pd.ExcelWriter(temp_path, engine='openpyxl') as writer: # ... 写入逻辑 os.replace(temp_path, 'output.xlsx') # 原子操作，避免损坏4.5 现象：Sup_wald_lag.xlsx中的 Sup-Wald 统计量为 NaN
原因：qgranger_test()中res_full.wald_test(R)在某些 τ 下因样本不足返回 NaN，未做兜底。
解决：在qgranger_test()的wald_stats.append(wald)前添加检查：
if np.isnan(wald) or np.isinf(wald): wald = 0.0 # 设为0，表示无证据拒绝H0 warnings.warn(f"τ={tau} 下Wald统计量无效，设为0") wald_stats.append(wald)5. 进阶技巧：用Sup_wald_lag.xlsx做动态因果图谱，三步定位“何时因果关系最强”
Sup-Wald 统计量不是单个数字，而是一个随滞后阶数变化的序列——Sup_wald_lag.xlsx的每一行对应一个滞后阶k（1~max_lag），列是各 τ 下的统计量。这能帮你回答：“X 对 Y 的因果效应，是即时生效（k=1），还是需要累积（k=3）？” 以下是实操三步法：
5.1 步骤一：加载并重塑数据，构建“滞后×分位”热力图
Sup_wald_lag.xlsx默认是宽表（lag 为行，τ 为列）。用pandas.melt()转为长表，再pivot()生成热力图所需矩阵：
import pandas as pd import seaborn as sns import matplotlib.pyplot as plt # 加载Sup-Wald结果 sup_df = pd.read_excel('Sup_wald_lag.xlsx', index_col=0) # 转为长表：lag, tau, statistic long_df = sup_df.reset_index().melt(id_vars='index', var_name='tau', value_name='statistic') long_df.rename(columns={'index': 'lag'}, inplace=True) # 重塑为热力图矩阵：行=lag，列=tau，值=statistic heatmap_data = long_df.pivot(index='lag', columns='tau', values='statistic') # 绘制热力图 plt.figure(figsize=(10, 6)) sns.heatmap(heatmap_data, annot=True, fmt='.2f', cmap='RdBu_r', center=0, xticklabels=[f'τ={t}' for t in heatmap_data.columns], yticklabels=heatmap_data.index) plt.title('Sup-Wald 统计量：滞后阶 × 分位点') plt.ylabel('滞后阶数 k') plt.xlabel('分位点 τ') plt.tight_layout() plt.savefig('sup_wald_heatmap.png', dpi=300, bbox_inches='tight') plt.show()提示：
cmap='RdBu_r'让正值（红色）和负值（蓝色）对比鲜明；center=0确保色标中心为 0，便于识别显著区域；fmt='.2f'保留两位小数，避免热力图数字拥挤。
5.2 步骤二：定位“因果峰值滞后”——对每个 τ，找 Sup-Wald 最大值对应的 k
热力图只能看趋势，要量化“最佳滞后”，需对每列（每个 τ）找argmax：
# 对每个τ，找Sup-Wald最大的滞后阶 peak_lags = {} for tau in heatmap_data.columns: max_stat = heatmap_data[tau].max() peak_k = heatmap_data[tau].idxmax() peak_lags[tau] = {'k_peak': peak_k, 'stat_max': max_stat} # 转为DataFrame便于分析 peak_df = pd.DataFrame(peak_lags).T print(peak_df) # 输出示例： # k_peak stat_max # τ=0.1 2 8.23 # τ=0.25 1 6.45 # τ=0.5 1 5.12 # τ=0.75 3 7.89 # τ=0.9 2 9.01注意：
k_peak=1表示即时因果，k_peak=3表示需三阶滞后才显现效应。若k_peak随 τ 增大而增大（如 τ=0.1 时 k=1，τ=0.9 时 k=3），说明极端事件的因果传导更慢，需在论文中重点讨论。
5.3 步骤三：绘制“因果强度轨迹”——用output.xlsx的 IRF 数据叠加 Sup-Wald 峰值点
真正的洞察在于关联静态检验（Sup-Wald）与动态响应（IRF）。将output.xlsx中的 IRF 数据（各 τ 各步的 0.05/0.5/0.95 分位）与peak_lags叠加，画出“因果强度轨迹图”：
# 加载IRF结果（假设已从output.xlsx读取到irf_dict） # irf_dict 格式：{tau: {var_name: {step: [q05, q50, q95]}}} taus = [0.1, 0.25, 0.5, 0.75, 0.9] steps = list(range(20)) plt.figure(figsize=(12, 8)) colors = ['red', 'orange', 'green', 'blue', 'purple'] for i, tau in enumerate(taus): # 提取y1的IRF中位数路径 irf_med = [irf_dict[tau]['y1'][s][1] for s in steps] # [1 <p> <a href="https://download.csdn.net/download/zru_9602/91459728" style="color:#ec7500;font-size:14px;"> 本文还有配套的精品资源，点击获取 </a> <img alt="menu-r.4af5f7ec.gif" src="https://csdnimg.cn/release/wenkucmsfe/public/img/menu-r.4af5f7ec.gif" style="width:16px;margin-left:4px;vertical-align:text-bottom;cursor:text;"> </p>
