简介这份资源是面向统计学、数据分析与科研人员的MATLAB因果分析源码合集聚焦线性结构因果模型、收敛交叉映射与向量自回归模型三类方法适合具备一定MATLAB基础、希望深入探究变量间因果而非仅相关性的研究者使用。压缩包共169个文件、约14.34MB以90个m脚本为核心配合22个gif动态演示、8个pdf理论文档、7个py脚本、6个txt说明及5个xlsx数据表另含mat数据、html与png结果文件覆盖建模、检验与推断的完整流程。目前已有357人学习下载。读者可据此复现线性结构因果模型的构建、用收敛交叉映射检测时间序列间的潜在因果并通过VAR模型分析多变量动态响应与反馈效应gif直观呈现计算过程pdf补充方法论与案例Python脚本则提供跨语言实现思路便于对照学习与二次开发。1. 因果分析在 MATLAB 里到底怎么落地从相关到因果的那一步很多人第一次做因果分析都是从一张相关性热力图开始的。两路信号相关系数 0.87看着挺高于是下结论说 A 影响 B。结果把数据滞后几拍再算相关系数掉到 0.1之前的结论直接翻车。问题出在相关是对称的而因果是有方向的还带时间延迟。MATLAB 里做因果分析设计源码核心就是把「方向」和「延迟」这两件事用可计算的指标固定下来。线性结构方法用自回归系数和偏相关判断谁先谁后收敛交叉映射处理非线性耦合向量自回归模型则给出多变量之间的动态响应。这三套东西不是互相替代而是覆盖不同数据假设线性高斯、非线性弱耦合、多变量系统。适合手里已经有时间序列、想从「一起动」推到「谁驱动谁」的工程师尤其是做故障溯源、传感器因果排序、生理信号耦合分析这类场景。下面按选型、实现、参数、避坑的顺序拆开讲。2. 线性结构方法先验图与偏相关怎么定方向2.1 为什么先做线性结构而不是直接上非线性线性结构方法的价值在于它便宜、可解释、能当基线。它的假设是变量之间关系近似线性用向量自回归拟合之后看系数矩阵的显著性来判断因果方向。常见做法是格兰杰因果检验但格兰杰因果只告诉你「X 的过去能帮助预测 Y」不等于物理因果。所以在线性结构这一步我一般会先构建一个先验因果图再用偏相关做条件独立性检验把明显不成立的边剪掉。偏相关的意义是在控制其他变量之后X 和 Y 还剩多少相关。如果控制 Z 之后 X 和 Y 的偏相关接近零那 X 到 Y 的边很可能是通过 Z 传导的不该直接连。这一步能避免后面 VAR 模型里出现一堆虚假显著系数。MATLAB 里算偏相关没有现成一行函数但可以用协方差矩阵求逆来推。下面这段代码是通用的偏相关计算输入是数据矩阵每列一个变量。function pcorr partial_corr(X) % X: n x p 数据矩阵每列一个变量 % 返回 p x p 偏相关矩阵 C cov(X); % 协方差矩阵 p size(X, 2); invC inv(C); % 协方差逆矩阵 D diag(1 ./ sqrt(diag(invC))); pcorr -D * invC * D; % 偏相关公式 pcorr(1:p1:end) 1; % 对角线置 1 end逻辑说明偏相关矩阵由协方差逆矩阵标准化得到负号是因为精度矩阵的非对角元与偏相关符号相反。参数上X 的列顺序就是变量编号后面画图要对应。如果变量量纲差异大先做 zscore 标准化否则协方差矩阵条件数会很差inv 出来的结果不稳定。拿到偏相关矩阵后用阈值筛边。阈值不是拍脑袋我一般用 Fisher z 变换做显著性检验样本量 n、变量数 p 已知时偏相关系数 r 的 z 值约等于 sqrt(n - p - 3) * 0.5 * log((1r)/(1-r))绝对值大于 1.96 就保留。这一步做完因果图基本就瘦下来了。2.2 用 VAR 系数矩阵做格兰杰检验的实操线性结构定完边接下来用向量自回归把动态关系拟合出来。VAR 模型形式是 Y_t c A1Y_{t-1} ... ApY_{t-p} e_t其中 A 矩阵的非对角元就是跨变量的影响系数。MATLAB 里可以用 varm 和 estimate也可以自己写最小二乘。自己写的好处是能控制滞后阶数和约束。function [A, resid] fit_var(Y, p) % Y: T x k 时间序列p: 滞后阶数 T size(Y, 1); k size(Y, 2); Ylag []; for i 1:p Ylag [Ylag, Y(p-i1:end-i, :)]; end Ytarget Y(p1:end, :); X [ones(size(Ylag,1),1), Ylag]; % 加常数项 B X \ Ytarget; % 最小二乘 A B(2:end, :); % 去掉常数项 resid Ytarget - X * B; end逻辑说明Ylag 按滞后从近到远拼接X 第一列是常数。B 是回归系数A 按变量块排列。参数 p 的选择用 AIC 或 BICMATLAB 里可以循环 p1 到 10算 AIC log(det(Sigma)) 2k^2p/T取最小。注意 k 是变量数p 太大时参数个数 k^2*p 会爆炸T 不够就过拟合残差白噪声检验会失败。格兰杰检验看的是 A 里对应块的联合显著性。比如检验变量 j 是否格兰杰引起变量 i就看 A 中所有滞后项里 i 行 j 列那组系数是否整体为零用 F 检验。这一步做完线性结构能给出的方向就齐了。3. 收敛交叉映射 CCM非线性耦合的检测与参数调优3.1 CCM 的核心思想与适用边界收敛交叉映射Convergent Cross Mapping是给非线性、弱耦合、不能做实验干预的系统准备的。它的逻辑是如果 X 引起 Y那么 Y 的历史里应该包含 X 的信息用 Y 的延迟嵌入去重构 X 的流形预测精度会随样本长度增加而收敛。反过来如果只是相关没有因果预测精度不会随长度收敛。CCM 的适用边界很明确系统要是确定性的、有低维吸引子、耦合不能太强也不能太弱。太强会同步太弱重构不出来。我一般先用互信息或者距离相关性判断有没有非线性依赖再上 CCM。如果数据是纯随机或者高噪声CCM 的收敛曲线会一直平着这时候别硬解释。MATLAB 里实现 CCM 主要分四步延迟嵌入、找近邻、交叉映射预测、算相关系数随长度的收敛。下面给核心函数。function rho ccm_score(X, Y, E, tau, L) % X 引起 Y 的检验用 Y 的嵌入预测 X % E: 嵌入维数, tau: 延迟, L: 使用的样本长度 X X(:); Y Y(:); N length(Y); % 构造 Y 的延迟嵌入 M N - (E-1)*tau; Yemb zeros(M, E); for i 1:E Yemb(:, i) Y((1:M) (i-1)*tau); end % 用前 L 个点做预测 L min(L, M); pred zeros(L, 1); for t 1:L dist sqrt(sum((Yemb - Yemb(t, :)).^2, 2)); dist(t) inf; % 排除自身 [~, idx] sort(dist); nn idx(1:E1); % E1 个近邻 w exp(-dist(nn) / (dist(nn(1)) eps)); w w / sum(w); pred(t) sum(w .* X(nn)); end rho corr(pred, X(1:L)); end逻辑说明Yemb 是 Y 的延迟坐标近邻用欧氏距离权重用指数核预测目标是 X。参数 E 用假近邻法或 FNN 确定tau 用自相关第一次过零或互信息第一个极小值。L 从短到长循环看 rho 是否随 L 增大而上升并趋于稳定。如果 rho 一直低或者不收敛说明 Y 里没有 X 的信息X 到 Y 的因果证据弱。3.2 参数 E、tau、L 怎么设才不玄学E 太小流形折叠近邻找错E 太大近邻稀疏预测方差大。我一般用 FNN 从 1 试到 10取假近邻比例低于 5% 的最小 E。tau 用互信息第一个局部极小值对应的延迟就是合适值通常在 1 到 10 之间。L 的循环范围从 20 到 N 的 80%画收敛曲线。有个血泪经验CCM 对噪声很敏感如果数据信噪比低先做滑动平均或者小波去噪但别过度否则会把非线性结构抹平。另外CCM 的 rho 绝对值不重要重要的是趋势。X 引起 Y 时用 Y 预测 X 的 rho 随 L 上升同时反向用 X 预测 Y 的 rho 可能也上升但更平。两个方向都算比较收敛斜率。4. 向量自回归模型多变量动态与脉冲响应4.1 VAR 的阶数选择与稳定性检验VAR 是线性结构方法的延伸但它的重点从「有没有因果」转到「因果有多强、持续多久」。做 VAR 之前必须检验平稳性单位根检验用 adftest 或 kpsstest。不平稳就差分差分阶数用 pp 检验定。平稳之后选滞后阶数AIC、BIC、HQ 三个准则一起看多数情况取 BIC 最小的阶数因为 BIC 惩罚更重模型更简约。maxP 10; aic zeros(maxP,1); bic zeros(maxP,1); for p 1:maxP [~, resid] fit_var(Y, p); Sigma cov(resid); k size(Y,2); T size(resid,1); aic(p) log(det(Sigma)) 2*k^2*p/T; bic(p) log(det(Sigma)) log(T)*k^2*p/T; end [~, p_opt] min(bic);逻辑说明det(Sigma) 是残差协方差行列式越小拟合越好但参数多了会过拟合所以加惩罚项。p_opt 就是选定的阶数。选完阶数后做稳定性检验看伴随矩阵的特征根是否都在单位圆内。如果有根接近 1说明模型不稳定脉冲响应会发散这时候要么降阶要么检查数据里有没有结构性断点。4.2 脉冲响应与方差分解的解读脉冲响应函数IRF描述一个变量受到冲击后其他变量怎么随时间响应。VAR 的 IRF 依赖变量排序排序不同结果不同这是最容易被误读的地方。我一般用 Cholesky 分解排序按「外生性」从强到弱或者用广义脉冲响应避免排序依赖。% 假设已估计 VARA 为系数Sigma 为残差协方差 horizon 20; k size(Sigma, 2); IRF zeros(horizon, k, k); P chol(Sigma, lower); IRF(1, :, :) P; for h 2:horizon for j 1:k acc zeros(k,1); for i 1:min(h-1, size(A,1)/k) acc acc A((i-1)*k1:i*k, :) * squeeze(IRF(h-i, :, j)); end IRF(h, :, j) acc; end end逻辑说明P 是 Cholesky 下三角代表当期冲击。IRF(h,:,j) 是第 j 个变量冲击后第 h 期各变量的响应。参数 horizon 一般取 20 到 50看响应是否衰减到零。如果 IRF 不收敛说明模型不稳定或者阶数不够。方差分解则是看每个变量预测误差方差被各冲击解释的比例通常和 IRF 一起看互相印证。5. 避坑与排查因果分析里最容易翻车的五件事5.1 现象格兰杰显著但物理上说不通原因格兰杰因果只反映预测增益如果两个变量都被第三个未观测变量驱动会出现虚假显著。解决把可疑的共同驱动变量加进 VAR或者用偏相关先筛条件独立。如果加进去之后显著性消失原来的边就是伪的。5.2 现象CCM 收敛曲线两个方向都上升原因耦合太强导致同步或者数据长度不够近邻统计不可靠。解决检查耦合强度如果两个方向 rho 都接近 1说明系统同步了CCM 失效。可以降采样或者用更长的数据但根本办法是换方法比如用转移熵。5.3 现象VAR 残差有自相关原因滞后阶数不够或者数据有未建模的周期性。解决先看残差自相关函数如果滞后 p 之后还有显著尖峰加阶数。如果加阶数后 AIC 还在降但模型不稳定检查数据有没有季节性做季节差分。5.4 现象脉冲响应置信区间宽到没法解读原因样本量小或者参数太多。解决用 bootstrap 算置信区间如果区间跨零说明响应不显著。减少变量数或者用贝叶斯 VAR 加先验收缩。别硬解释一个跨零的 IRF。5.5 现象偏相关矩阵求逆报奇异原因变量高度共线协方差矩阵条件数太大。解决先算条件数大于 1e6 就说明共线严重。删掉冗余变量或者用伪逆 pinv 代替 inv但伪逆会改变偏相关的解释慎用。更稳的做法是先做 PCA 降维但降维后变量物理意义会模糊。6. 把三套方法串成一条流水线验证与进阶技巧三套方法单独用都有盲区串起来才稳。我一般的流水线是先做平稳性检验和去趋势然后偏相关筛边得到初始因果图接着用 VAR 拟合线性部分并做格兰杰检验对线性解释不了的残差再做 CCM 看非线性耦合最后用脉冲响应和方差分解量化动态强度。这条流水线的好处是每一步都有验证不会一上来就陷进非线性方法里出不来。验证环节我固定做三件事。第一置换检验把 X 的时间顺序随机打乱重算因果指标重复 500 次看真实指标在零分布里的分位数。如果 p 值大于 0.05不管指标多好看都不认。第二分窗稳定性把数据切成前后两半分别算因果指标如果符号相反或者差异巨大说明关系不稳定可能是非平稳或者有断点。第三反向检验X 到 Y 显著时一定算 Y 到 X如果两个方向都显著且强度接近要警惕共同驱动或者同步。% 置换检验示例检验 X 到 Y 的格兰杰因果 nPerm 500; F_real granger_F(Y, X, p_opt); % 自定义函数返回 F 统计量 F_perm zeros(nPerm, 1); for i 1:nPerm X_shuffle X(randperm(length(X))); F_perm(i) granger_F(Y, X_shuffle, p_opt); end p_value mean(F_perm F_real);逻辑说明granger_F 内部做两次 VAR 回归一次含 X 滞后项一次不含比较残差平方和。置换打乱 X 的时间顺序破坏因果但保留边缘分布。p_value 小于 0.05 才认为显著。参数 nPerm 至少 500太少分位数不准。进阶技巧方面如果变量数超过 10 个VAR 参数会爆炸这时候用稀疏 VAR加 L1 正则MATLAB 里可以用 lasso 逐变量做。CCM 的多变量版本可以用多元嵌入但维数灾难严重我一般先做变量聚类每组选一个代表。另外所有因果指标都对时间延迟敏感建议在 0 到 20 拍范围内扫一遍延迟取指标最强的延迟作为主延迟但别只报最强的那一个把延迟曲线一起放出来读者才能判断稳健性。最后说个我自己的习惯每次跑完因果分析先把所有中间结果存成 mat 文件包括偏相关矩阵、VAR 系数、CCM 收敛曲线、IRF。过一周再回来看如果换个预处理或者换个阶数结论就变那这个结论本来就不该信。因果分析最怕的就是只报一个 p 值把过程藏起来。把过程摊开自己和自己对账比任何花哨方法都管用。希望帮到你。本文还有配套的精品资源点击获取