简介马尔科夫转换GARCHMS-GARCH模型是融合GARCH与马尔科夫链的金融波动率分析工具这套MATLAB实现包恰好为研究者与量化学习者提供了可直接运行的代码基础。压缩包共14个文件、约1.34MB其中12个M文件覆盖主程序、似然估计、协方差矩阵与起始值设定等核心模块辅以1份Markov-switching EGARCH的PDF说明文档和1个文本文件便于对照公式理解实现细节。已有490人学习浏览适合具备一定GARCH基础、希望进一步掌握状态转换建模的读者。通过运行main.m并调用相关似然函数可完成MS-GARCH/MS-EGARCH的构建、最大似然估计与波动率模拟也可基于协方差矩阵模块扩展研究是复现经典模型、开展实证分析的高效起点。1. 金融时间序列的状态切换难题这套 MS-GARCH 代码包解决的是什么做金融时间序列实证的人多半会遇到同一个尴尬用单机制 GARCH 去拟合一段跨越牛熊的收益率数据波动率估计往往“不上不下”——既不像恐慌期那么高也不像平静期那么低。马尔科夫转换 GARCHMS-GARCH就是为这种场景准备的它把波动过程拆成若干个状态状态之间按转移矩阵切换每个状态内部再套 GARCH 或 EGARCH 结构。这套名为 markov garch 的 MATLAB 代码包包含主程序、似然函数变体、初值设定、协方差矩阵估计和一份 Markov-switching EGARCH 的 PDF 说明文档基本覆盖了从数据预处理到参数估计再到诊断的完整链路。适合正在写实证论文、做波动率预测或风险度量的人尤其是已经会用 GARCH但想往状态转换方向扩展的 MATLAB 用户。下面直接从代码结构入手讲清楚这套资源怎么用、坑在哪。2. 把波动拆成两套机制MS-GARCH 的状态设定与三种似然函数选型2.1 状态数与转移矩阵为什么两状态是默认起点MS-GARCH 的核心假设是收益率序列在每一个时点处于某个不可观测的状态状态之间的转移服从一阶马尔科夫链。最常见的设定是两个状态一个代表低波动比如平稳期一个代表高波动比如危机期转移矩阵是 2×2% 转移矩阵 P(i,j)从状态 i 转移到状态 j 的概率 P [0.98 0.02; % 状态1低波动留在原地概率0.98跳到状态2概率0.02 0.05 0.95]; % 状态2高波动跳到状态1概率0.05留在原地概率0.95逻辑说明每一行加起来等于 1对角元越大状态越持久。金融数据里波动聚集的典型特征就是“高波动期会持续一段时间”所以对角元一般设成 0.9 以上。参数说明如果跑出来的转移概率对角元小于 0.5说明状态划分大概率有问题——两个状态在频繁互跳等于没分类。这套代码里状态数的设定不是直接写一个k 2那么简单的。main.m 里会初始化状态转移矩阵然后把它作为参数传给似然函数。改三状态的时候转移矩阵变成 3×3同时每个状态的 GARCH 参数要多复制一份GARCH 和 EGARCH 变体都要重新检查维度。选择两状态作为起点不是因为它最准确而是因为它最容易估计。三状态模型的参数个数几乎翻倍对样本长度和初值质量的要求高很多经常出现某个状态拟合到样本尾部的一两笔极端值上导致参数估计毫无经济含义。从两状态跑通再逐渐加状态是实操中最稳的路径。2.2 likelihoodfreiqml.m 与 likelihoodgarchgarchqml.m非对称性从哪来文件列表里的likelihoodfreiqml.m和likelihoodgarchgarchqml.m是两个不同的似然函数变体。从命名上看frei在德文里是“自由”的意思对应英文的free所以likelihoodfrei大概率是指“自由参数版本”的似然——也就是不约束某些系数让优化器自己去估计。而likelihoodgarchgarch对应的是两个状态都设定为 GARCH(1,1) 的版本。qml后缀是 quasi-maximum likelihood 的缩写即拟最大似然估计。三个变体的核心计算框架是一致的function [logL, lik] likelihoodgarchgarchqml(theta, data, k) % theta: 参数向量包含 GARCH 参数和转移矩阵元素 % data: 收益率序列列向量 % k: 状态数量 % 拆出 GARCH 参数omega(1:k), alpha(1:k), beta(1:k) % 拆出转移矩阵 Pk×k % 初始化滤波概率和条件方差 % 逐点递推先算条件方差再算正态密度然后更新状态概率 % logL sum(log(lik)); end逻辑说明这就是 Hamilton 滤波的 GARCH 版本。每一个时间点模型对每个状态算一个条件方差然后按状态概率加权得到混合正态分布的密度值最后在样本上累加得到似然值。参数说明theta的排列顺序必须和startvalues.m里定义的顺序完全一致否则估计结果对不上数据必须是收益率列向量不是价格序列这是初学者最容易犯的一个错误。likelihoodfreiqml.m和likelihoodgarchgarchqml.m的区别点我一般理解为前者在某些系数上做了自由化处理。比较典型的做法是允许两个状态的 GARCH 截距项各自独立地估计而后者可能是共享部分参数的约束版本。具体要看代码里的注释和参数向量长度。两套版本跑同样的数据如果结果差异很大多半是约束条件没有正确设置而不是模型本身的问题。2.3 EGARCH 变体的似然函数为什么偏好非对称设定likelihoodegarchegarch.m和likelihoodegarchgarch.m这两个文件把 EGARCH 引了进来。EGARCH 和 GARCH 的关键区别在于GARCH 对方差方程直接建模EGARCH 对对数方差建模并且加入了一个非对称项。这个非对称项非常贴合金融数据的一个经典特征利空消息对波动率的冲击比利好消息更大。% EGARCH(1,1) 方差方程log(sigma_t^2) % omega alpha * (|eps_{t-1}|/sigma_{t-1}) gamma * (eps_{t-1}/sigma_{t-1}) beta * log(sigma_{t-1}^2) % 其中 gamma 项就是非对称项eps 为负时gamma*eps/sigma 为负对 log 方差产生额外冲击逻辑说明这个公式里gamma的符号和大小决定了非对称性的方向与强度。如果估计出的gamma显著为负说明负残差对波动率的推高作用强于正残差符合“坏消息引起更大波动”的典型事实。参数说明EGARCH 里因为log(sigma^2)是可正可负的优化时不需要像普通 GARCH 那样强制约束方差为正数值稳定性更好这也是很多实证研究偏好 EGARCH 的原因。likelihoodegarchgarch.m的意思应该是状态 1 用 EGARCH、状态 2 用 GARCH这种混合设定在实践中是有意义的一个状态里可以存在非对称效应另一个状态里可以不存在。这种灵活度是 MS-GARCH 相对单一 GARCH 模型的优势所在。跑这类变体的时候注意状态编号和模型形式的对应关系别把 EGARCH 的参数套到 GARCH 的状态上。3. 跑通全套代码从 startvalues.m 到 main.m 的完整复现步骤3.1 startvalues.m初值怎么给才不崩MS-GARCH 的似然函数是非凸的局部最优解很多初值给不好优化器可能直接发散到 NaN。startvalues.m做的事就是把参数向量变成一个合理的初始值。常见的做法是先跑一个单一机制的 GARCH(1,1)或 EGARCH把估计结果作为状态参数的初值然后设定两个状态的截距项往上和往下偏离转移矩阵给一个“高持久”的设定function theta0 startvalues(data, k, modeltype) % 先用单一机制 GARCH 估计得到基准参数 % 以基准参数的 1.2 倍作为状态1初值、0.8 倍作为状态2初值 % 转移矩阵初始化为对角元 0.95、非对角元 0.05 % theta0 [omega_1, alpha_1, beta_1, omega_2, alpha_2, beta_2, P_11, P_22]; end逻辑说明这种“从单机制扩散到多机制”的做法能保证初值落在合理的参数区域内。盲目随机初始化会让优化器把参数推到似然面的边界上。参数说明初值只影响优化路径不影响最终估计结果的理论性质但实际计算中差初值经常导致估计卡在局部最优上所以这一步值得多花时间。几种常被验证有效的初值策略一是用样本方差作为两个状态共同的截距起点二是把估计出的单机制参数加上正负 20% 的扰动作为双状态初值三是固定转移矩阵跑几步再放开估计避免转移概率在初始阶段就粘在边界上。这套代码里 startvalues.m 就是类似逻辑改状态数的时候注意参数向量的长度要同步调整。3.2 maximierung_likelihood.m 的优化器设置fmincon 与 fminunc 的选择maximierung_likelihood.m是优化入口负责调用 MATLAB 的优化器最大化似然函数。文件名的 “maximierung” 是德文的“最大化”。默认的做法是用 fmincon 做约束优化因为 GARCH 参数需要满足约束条件omega 0、alpha 0、beta 0、alpha beta 1或者 EGARCH 里对参数符号的约束。% 设置优化选项 options optimoptions(fmincon, ... Algorithm, interior-point, ... Display, iter, ... MaxIterations, 1000, ... MaxFunctionEvaluations, 5000); % 参数下界和上界omega 和 alpha 的下界可以设一个很小的正值如 1e-6避免贴零 lb [1e-6 * ones(1, 2*k), 0, 0, -1 * ones(1, k)]; % 下界 ub [10 * ones(1, 2*k), 1 - 1e-4, 1 - 1e-4, 1 * ones(1, k)]; % 上界逻辑说明MaxFunctionEvaluations设 5000 是常规做法如果数据量大或状态数多可能不够需要调高到 1e5 甚至更高。用Displayiter可以看到优化迭代过程发现似然值长期不降或者卡在原地说明步长设置有问题或初值不理想。参数说明fmincon 的 interior-point 算法对非线性约束的处理较稳妥但计算量比 fminunc 大如果不涉及复杂约束也可以试试 fminunc速度会快不少。一个重要的优化细节MS-GARCH 似然函数里经常出现概率密度接近零的情况导致log(lik)变成负无穷。很多实现会在似然函数里加一个极小的下界保护比如lik max(lik, 1e-300)防止优化器在似然值为零的区域反复试探。这套代码里如果没有这个保护遇到极端数据时会崩可以自己补上。3.3 covariancematrix.m标准误、Hessian 与参数显著性covariancematrix.m负责在估计完成后计算参数的标准误。标准做法是用数值 Hessian 矩阵的逆配合 QML 框架下的 sandwich 形式即稳健标准误。在 QML 设定下数据不要求严格服从正态分布标准误需要在普通 Hessian 矩阵基础上做稳健化修正% 在最优参数处计算数值 Hessian H hessian((theta) likelihoodgarchgarchqml(theta, data, k), theta_hat); % 稳健标准误 sqrt(diag(inv(H) * J * inv(H))) % J 是得分向量的外积和逻辑说明如果不做 sandwich 修正标准误会偏小参数显著性检验容易被误导。很多论文里报告显著系数实际上用的就是 QML 标准误。参数说明hessian函数需要自己写或用工具箱里的fminunc输出网格步长的选择对结果有影响一般取参数值乘以 1e-4 作为最小步长。标准误估计是整个代码包里的隐藏要点。不同似然函数变体GARCH vs EGARCH对应的 Hessian 结构不同covariancematrix.m在读取参数向量时要注意当前用的是哪个变体避免参数索引错位导致标准误全错。3.4 main.m 的输入输出规范收益率序列怎么预处理main.m是整套流程的入口。它负责读数据、调 startvalues、调优化器、调 covariancematrix最后把结果拼成一份可读的输出。实际上手时数据处理这一步最容易翻车% 读取价格数据计算简单收益率 prices readtable(your_data.csv); returns diff(log(prices.Close)); % 对数收益率 returns returns(~isnan(returns)); % 去掉缺失值 returns returns * 100; % 转成百分数数值上更稳定逻辑说明用对数收益率的理由是它在时间上加总等于总收益且分布性质比简单收益率更接近正态假设。乘以 100 是因为收益率数值太小0.001 级别会让优化器的梯度在数值精度边缘徘徊乘完以后参数量级更友好。参数说明如果你用的是日频数据可以尝试乘以 100如果是低频月度数据可以不乘。一个常见的预处理误区直接用原始价格序列跑 MS-GARCH。价格序列是非平稳的方差随时间增长估计出的波动率没有任何意义。另一个高频错误是忘记去均值——收益率序列最好先减去样本均值或拟合一个常数均值否则均值会吸收一部分波动率信息导致 GARCH 参数偏向 0。main.m 跑完后会输出各状态的条件方差、状态滤波概率和参数估计。我的习惯是先把滤波概率画出来看每个状态对应的样本分段是否与实际市场周期大致对应。如果两个状态完全随机穿插说明模型设定有问题——状态数太多、参数有冗余或数据不适合用 MS-GARCH。4. 避坑MS-GARCH 估计的五个高频翻车现场4.1 初值太差导致似然值 NaN 或虚数现象模型跑起来第一轮迭代目标函数直接返回 NaN优化器报错退出。原因GARCH 方差方程的递归更新中某个时点的条件方差变成了负数或零取对数后得到 NaN。常见触发点包括初始方差设成 0、alpha 和 beta 为负、数据里存在极端离群值。解决初始方差用样本方差而非 0不等式约束里给 alpha、beta 加下界 1e-6如果数据里有正负 10 倍标准差的极端值做截尾处理或换 EGARCH天然避免方差为负。另外在似然函数内加if ~isfinite(logL), logL 1e12; return; end的保护性分支让优化器跳过无效区域。4.2 状态标签翻转估计出的状态 1 和状态 2 定义反了现象代码里设定状态 1 是低波动跑出来状态 1 的均值方差反而更大低频波动率高于高频状态。原因似然函数的标签不具可识别性——把状态 1 和状态 2 互换得到的似然值完全一样。优化器可能收敛到任一标号方向和你的先验设定无关。解决估计完成后加一步标签重排比较两个状态的 GARCH 截距或无条件方差如果状态 1 的方差大于状态 2就对换参数向量的对应区块和转移矩阵的行列。画图的时候按“高波动 vs 低波动”来命名状态而不是按代码里的“状态 1 vs 状态 2”。这也是为什么看结果时不能只看状态编号要看参数数值。4.3 QML 标准误与普通标准误混用现象参数估计收敛但报告的标准误小得离谱t 值动辄几十。原因直接用数值逆 Hessian 当协方差矩阵忽略了 QML 情境下的 sandwich 修正。金融收益率通常有厚尾正态拟似然设定下的误差项并非真实正态分布普通标准误低估了不确定性。解决统一按inv(H) * J * inv(H)计算。检查covariancematrix.m的输出是不是 sandwich 形式如果不是在main.m后处理阶段补上。注意很多情况下 J 矩阵的数值波动很大建议用外积梯度在移一个样本点后的量级做平滑而不是直接逐点累加。4.4 EGARCH 参数的指数运算溢出现象EGARCH 变体跑着跑着条件方差变成 Infcovariancematrix 计算崩溃。原因EGARCH 的方差方程里exp(omega ...)的指数部分在参数迭代中变得很大exp 直接溢出到 Inf。这种情况在赫赛矩阵计算时尤为常见因为数值差分会让参数暂时漂移到极端区域。解决给指数运算上限sigma2 min(sigma2, 1e8)或者在目标函数外部包一层if any(~isfinite(sigma2))的检查。另一种方案是修改参数化方式把 omega 写成log(omega)让优化器在无约束空间发挥作用估计完再转化回来。4.5 转移矩阵收敛到边界导致状态退化现象估计出的 P_11 或 P_22 等于 0.9999状态几乎不切换或者某个状态的对角元掉到 0.5 以下状态变成高频噪音。原因数据里如果某个状态持续时间极短转移概率会推到边界上。样本长度不足是主因其次是状态之间差异不明显模型无法可靠区分。解决一是增加约束把转移概率限制在 [0.85, 0.995] 区间内二是减少状态数从 2 到 2 再试其他结构三是检查每个状态的平均持续期1/(1-P_ii)——小于 5 天的状态基本没有实证意义。另一种做法是用 t 分布替代正态分布作为条件分布厚尾可以吸收部分极端值避免模型为了拟合尾部极端值而产生退化状态。5. 把代码改造成通用模块从单数据集到批量实证研究5.1 统一接口函数用结构体收纳全部参数原始代码的传参方式是散装的参数向量、数据、模型类型各传各的跑多了容易乱。建议包装一个统一的调用入口把所有选项装进一个结构体function results run_msgarch(data, opts) % opts 字段 % .k 状态数量默认 2 % .model garch / egarch / mix % .start_method singleregime / grid % .optimizer fmincon / fminunc % .max_iter 最大迭代次数 % .print_level 输出级别 if ~isfield(opts, k), opts.k 2; end if ~isfield(opts, model), opts.model garch; end % 后续按 opts.model 选择对应的似然函数和约束条件 % 估计完成后统一生成标准化输出结构体 results end逻辑说明把原来的 main.m 改成函数式调用后批量跑多组数据就只需要一个 for 循环。参数说明start_method设为grid时代码会用多组初值并行起跑选出似然值最高的作为最终估计——这个策略对 MS-GARCH 这类多峰问题非常有效虽然耗时翻倍但能显著降低局部最优概率。5.2 多状态与多滞后扩展两状态到三状态的改动点切换三状态时不是简单把k改成 3 就完事有四处必须检查第一处转移矩阵从 2×2 变成 3×3参数向量里多出 3 个转移概率。第二处每个状态的 GARCH 参数块复制成三份omega、alpha、beta长度都变成 k。第三处滤波递推矩阵的维度全部跟着改Hamilton 滤波里alpha_t向量长度变 3P_ij的求和逻辑要重写。第四处状态初始概率pi不再是二元素向量而是三元素向量。全部改完后用模拟数据验证状态识别效果再上真实数据。实践中三状态最常见的设定是“低波动 / 中波动 / 高波动”转移矩阵的秩结构不再是简单对角而是允许跳跃式转换。比如低波动可以直接跳到高波动危机冲击而不一定经过中波动。这种跳跃设定在代码层面就是在转移矩阵中让某些非对角元素也可能显著不为零。5.3 平滑概率与状态划分从估计结果到可解释的结论估计完成后main.m输出的滤波概率是P(state_t | data_1:t)也就是只用了过去信息的状态推断。实务上更常用的是平滑概率P(state_t | data_1:T)它用完整样本的信息重新推断每一个时点的状态识别结果更稳定。平滑概率的实现并不复杂% 平滑概率递推从 T 到 1 反向计算 % sm(:,T) filter(:,T) % for t T-1:-1:1 % sm(:,t) filter(:,t) .* (P * (sm(:,t1) ./ (P * filter(:,t1)))); % end逻辑说明这个反向递推是 Kim (1994) 的平滑算法计算量不大但能显著改善状态划分质量。参数说明平滑概率里如果出现除零说明某个时刻滤波概率所有元素都接近零需要加上下界保护。很多实证论文里的状态图严格来说应该用平滑概率而非滤波概率。把平滑概率按时间序列可视化可以看出状态切换的时点是否与经济事件吻合。我的经验是如果平滑概率在 0.5 附近反复徘徊说明模型对这个时点属于哪个状态没有把握与其强行划分不如报告这个不确定性。6. 参数识别的一个硬检查重估扰动法评估状态划分置信度MS-GARCH 最容易被审稿人或合作者质疑的一点是“你划分出来的那两个状态是真的存在还是模型自己编出来的”一个可行的验证方法是重估扰动法也叫 multistart 检查法。做法是固定数据不动生成 50 组不同的初值分别跑完整估计流程记录每一组的最终似然值、参数和状态划分结果。然后看两个东西第一50 组估计的最大似然值是否集中在同一水平——如果分散得厉害说明优化没有收敛到同一区域第二各组估计出的状态序列是否高度相似——如果状态时点错开说明模型对该数据的状态划分缺乏唯一性。% 伪代码多初值重估 results_list cell(1, 50); for i 1:50 theta0 random_initialization(data, k); results_list{i} estimate_msgarch(data, theta0); end % 比较所有估计的 logL找到最高值及对应参数 % 检查处于最高 logL 附近比如相差 1 以内的参数是否一致逻辑说明这个检查的本质是验证估计结果不是某一个局部最优的偶然产物。参数说明50 组初值是比较平衡的选择样本量小或状态数多时可以加到 100 组虽然耗时但值得。做这个检查时不要只看似然值要同时看状态划分的交叉表——如果两组估计的 logL 很接近但状态完全错位说明识别性存疑。这个方法也适用于另一个场景给模型选状态数。如果 k2 在重估时稳定复现k3 的估计结果在 50 组初值中频繁翻车且不同初值收敛到不同局部最优那就是数据支撑不了三状态的有力信号。这套 markov garch 代码包里把所有似然变体和工具都配齐了做这种 robust check 不需要改太多代码只要把 main.m 包进一个循环就行。从那以后我每次拿到一套新数据跑 MS-GARCH都会强制走一遍重估扰动这个步骤看到 50 组结果的高度集中才敢写结论。状态划分这种事模型给出漂亮的分界线是一回事能不能经得起初值扰动还稳定是另一回事。希望这个检查方法和前面这些参数设置与避坑经验能帮到你省一点反复调试的时间。本文还有配套的精品资源点击获取