简介面向金融时间序列分析的一份MATLAB实现资源主要围绕马尔科夫转换GARCH模型即MS-GARCH用于刻画经济扩张或收缩、市场情绪高低等状态下波动率的非线性切换。该资源适合金融计量经济学研究者、量化分析师及研究生作为算法学习与复现参考。压缩包共14个文件以12个MATLAB脚本为主分别用于似然函数计算、最大似然估计、协方差矩阵求解、起始值设定等环节另含1篇Markov-switching EGARCH的Matlab说明文档和1个文本格式的似然函数文件帮助对照理解。包体仅1.34MB内容紧凑便于快速上手。目前已有490人学习下载受到一定关注。从代码结构看内容覆盖MS-GARCH与MS-EGARCH两类实现包含QGML估计等关键步骤并附有教学型PDF。读者可借助脚本快速搭建波动率状态转换模型进行估计、模拟与预测也能在此基础上修改设定扩展自身研究整体实用性强。1. 先从“波动率突然放大”说起马尔科夫转换 MS-GARCH 解决了什么问题2015年夏天A股、2020年3月全球市场波动率不是缓慢爬升的是“跳”起来的。用单一GARCH(1,1)拟合这类序列αβ经常跑到0.98以上意思是任何一次冲击都需要上百个交易日才能消化这与实际感受完全不符。马尔科夫转换模型Markov Switching的思路很简单让数据存在两种或多种状态——低波动状态和高波动状态状态之间的切换由一个随时间演化的概率矩阵控制在每个状态内部再用GARCH刻画短期的波动聚集。把这个结构拼起来就是标题里的 MS-GARCHMarkov-Switching GARCH。这个模型不新鲜但在量化风控和波动率预测场景里一直很能打它能给出“当前处在哪个波动区制”的后验概率也能做比单一GARCH更贴近现实的多步方差预测。适合谁做波动率预测、期权定价、目标波动率调仓、以及需要区分“常态波动”和“危机波动”的风控从业者。2. MS-GARCH 的模型结构为什么是“马尔科夫转换 GARCH”而不是别家组合2.1 单一 GARCH 的 αβ 聚合问题波动聚集与结构突变被混在一起先说一个我经常拿来当开场白的现象。对沪深300日收益率直接拟合GARCH(1,1)在2015到2022这个区间上α估计值一般在0.06到0.1之间β在0.9到0.94附近αβ≈0.98甚至更高。这意味着条件方差的衰减半衰期大约在ln(0.5)/ln(0.98)≈35天一次大波动要五周才能衰减一半。但如果你去看波动率曲线的实际形状2015年6月那轮跳升之后波动率在两个月内就明显收敛了2016年初熔断前后又来了一波。这种反复“跳起来、落下去”的形态单一GARCH只能用一个高持久性参数去平均化拟合结果就是平时高估波动惯性危机后高估波动持续。Markov-Switching的想法不是修改GARCH方程而是把样本分成几个“波动区制”每个区制有自己的GARCH参数。低波动区制里αβ可能只有0.85冲击衰减快高波动区制里αβ可能接近0.99表示危机中波动聚集更强烈。这样就把“结构性变化”和“短期的波动聚集”分离开来。2.2 马尔科夫转换的数学形式S_t、转移概率矩阵 P、平稳分布MS-GARCH 的核心是一个不可观测的状态变量 S_t它只能取有限个离散值。两区制情况下 S_t1 表示低波动S_t2 表示高波动。状态的变化由一个齐次马尔科夫链控制P [p11, 1-p11; 1-p22, p22]其中 p11 Pr(S_t1 | S_{t-1}1)p22 Pr(S_t2 | S_{t-1}2)。对角元素衡量状态持续性p11 和 p22 越高系统越“恋战”一旦进入高波动区制就不容易出来。稳态概率 π1 (1-p22)/(2-p11-p22)π2 (1-p11)/(2-p11-p22)这个公式经常用来做初值设定。观测方程每个区制内仍是标准GARCH(1,1)σ_t²(s) ω_s α_s * ε_{t-1}² β_s * σ_{t-1}²(s)注意这里条件方差是带上标(s)的。问题在于t 时刻的 σ_t² 理论上依赖从 t0 到 t 的完整状态路径路径数量是 2^t这就是后面避坑章节要单独说的路径依赖问题。实现上必须做近似我在第3章会给出常用的近似估计方法。2.3 区制数K2 还是 K3大多数实证用两个区制就够了一个代表“平静期”一个代表“危机/高波动期”。两区制的优势是参数的局部识别相对容易转移概率和后验状态概率都比较好解释。三个区制也不是不能用它通常会把低波动进一步拆成“温和波动”和“极低波动”或者把高波动拆成“高波动”和“极端行情”。但K3的代价是GARCH参数从6个变成9个转移矩阵从2个参数变成6个样本量不足的时候估计结果对初值高度敏感而且容易出现两个高波动区制互相抢数据的翻车现场。我的经验是先用K2跑通全流程再拿AIC/BIC去看K3是否有实质改善如果BIC改善不到10个点就继续用K2。2.4 均值要不要也切换MS-GARCH 与 MS-AR-GARCH 的边界标题是MS-GARCH但落在“马尔科夫转换”这个大框架下均值方程的处理是个绕不开的选型点。最纯粹的形式是均值不切换μ 对所有区制相同区分完全靠波动率结构。好处是干净利落也是R语言MSGARCH包默认的做法。但实际收益率经常出现“高波动区制里均值更负”的现象比如2015年下半年、2018年波动率上行往往伴随负收益。这种情况下如果你只切波动率低波动区制的残差会被负收益拉偏导致方差估计失准。所以很多从业方案包括我在实际项目里的做法会这么选先跑纯 MS-GARCH看低波动区制残差均值是否显著不为零如果显著再升级到区制均值也切换的 MS-AR-GARCH。要注意的是均值切换会显著增加似然面的复杂度也让标签问题更严重所以不是越复杂越好。另外一个常被忽略的坑是条件分布。GARCH在险价值建模里正态分布假设会让高波动区制的尾部被低估。建议至少用 t 分布或偏 t 分布sstd做敏感性测试R的MSGARCH里可以直接指定distributionstd。我用正态只是为了让代码和讲解更直观。3. 用 Python 手写两区制 MS-GARCH(1,1)模拟数据、Hamilton 滤波、最大似然估计3.1 先造一份“答案已知”的数据模拟两区制 GARCH 序列模型没跑通之前先不要碰真实行情。我的习惯是先用已知参数模拟数据这样能确认估计器能不能恢复真值也方便后面排查。模拟步骤是从状态序列开始每个时点的条件方差由当期的状态决定然后按正态分布抽收益。import numpy as np import pandas as pd np.random.seed(42) n 2000 # 真实参数区制1低波动区制2高波动 mu_true 0.05 omega np.array([0.02, 0.50]) alpha np.array([0.05, 0.20]) beta np.array([0.90, 0.70]) p11_true, p22_true 0.98, 0.95 P_true np.array([[p11_true, 1 - p11_true], [1 - p22_true, p22_true]]) # 平稳分布采样初始状态 pi0 np.array([(1 - p22_true) / (2 - p11_true - p22_true), (1 - p11_true) / (2 - p11_true - p22_true)]) state np.zeros(n, dtypeint) state[0] np.random.choice([0, 1], ppi0) r np.zeros(n) h np.zeros(n) h[0] omega[state[0]] # 用一个非负初值 for t in range(n): if t 0: # 状态切换 state[t] np.random.choice([0, 1], pP_true[state[t-1]]) # 当前状态下的GARCH(1,1) h[t] omega[state[t]] alpha[state[t]] * r[t-1]**2 beta[state[t]] * h[t-1] r[t] mu_true np.sqrt(h[t]) * np.random.randn() ret pd.Series(r, nameret) print(ret.describe())逻辑说明这里用 r[t-1]² 而不是 (r[t-1]-mu)² 做模拟是为了让代码短一些实际估计时我用中心化残差。α_true(区制1)0.05 表示低波动区制内的杠杆效应很弱区制2的α0.20意味着高波动区制内“当日冲击”对次日方差的贡献大得多。β在不同区制的差异体现的就是波动衰减速度不同。实际项目里很少有人真的关心“真实参数”这个模拟环节主要用来检验估计器是否可识别如果模拟数据都估不准真实数据上更没法信。3.2 写似然函数用 Hamilton 滤波做状态概率递推用 Klaassen 近似处理路径依赖MS-GARCH 的似然难点不在马尔科夫链而在条件方差依赖状态历史。严格求和需要对 2^T 条路径做平均这在长度 2000 的序列上是天文数字。工程上的通行近似是 Klaassen2002的 collapsion 做法计算条件方差时用上一期每个状态下的滤波概率对 h[t-1] 做加权平均权重要乘上状态转移概率。from scipy.stats import norm def ms_garch_neg_loglik(params, r, n_states2): mu params[0] omega params[1:1n_states] alpha params[1n_states:12*n_states] beta params[12*n_states:13*n_states] p11, p22 params[-2], params[-1] P np.array([[p11, 1-p11], [1-p22, p22]]) T len(r) eps2 (r - mu) ** 2 # 初始状态概率P 的平稳分布 denom_ss 2 - p11 - p22 pi0 np.array([(1-p22)/denom_ss, (1-p11)/denom_ss]) filt np.zeros((T, n_states)) h np.zeros((T, n_states)) h[0] np.var(r) # 无条件样本方差作为共同起点 ll 0.0 for t in range(T): if t 0: dens norm.pdf(r[0], mu, np.sqrt(h[0])) num pi0 * dens ll np.log(num.sum()) filt[0] num / num.sum() else: # 一步预测状态概率Pr(S_t | F_{t-1}) pred filt[t-1] P # Klaassen 近似给定当前状态 s对上一期 h 做概率加权 for s in range(n_states): p_prev_given_s P[:, s] * filt[t-1] e_h (p_prev_given_s * h[t-1]).sum() / p_prev_given_s.sum() h[t, s] omega[s] alpha[s] * eps2[t-1] beta[s] * e_h # 当前状态下残差的分布密度 dens norm.pdf(r[t], mu, np.sqrt(h[t])) num pred * dens ll np.log(num.sum()) filt[t] num / num.sum() return -ll这段代码有几点要仔细说。第一为什么不是 h[t,s] ω_s α_s * ε_{t-1}² β_s * h[t-1,s]因为 t-1 时刻的真实方差并不是“某一个状态下的方差”而是对所有状态的不确定性加权平均。如果上一期滤波概率显示有30%概率处于高波动状态那么这一期低波动区制的方差起点也应该带上那30%的惯性。Klaassen 近似的本质就是把路径依赖用一个条件期望替代代价是似然不再是精确的但优化稳定性大幅提高。第二filt[t] 的更新是标准的 Hamilton 滤波先做一步预测 pred[t] filt[t-1] × P再乘以当前观测密度归一化后得到后验状态概率。这个概率序列本身就是重要的输出后面做区制识别和风控信号就是用它。第三初始状态概率用平稳分布。注意当 p11 或 p22 接近1时denom_ss 接近0pi0 会不稳定。我一般会加一个下限约束p11、p22 都落在 [0.5, 0.999] 区间内。3.3 用 SLSQP 做极大似然估计参数边界与 αβ1 约束优化时GARCH 参数必须满足 ω0、α≥0、β≥0、αβ1转移概率落在开区间(0,1)。这里 αβ 小于1是线性约束scipy.optimize.minimize 里用 SLSQP 同时支持 bounds 和 constraints。from scipy.optimize import minimize def fit_ms_garch(r, init_paramsNone): var0 np.var(r) if init_params is None: init_params np.array([ 0.0, # mu 0.05 * var0, 0.20 * var0, # omega1, omega2 0.05, 0.20, # alpha1, alpha2 0.90, 0.70, # beta1, beta2 0.90, 0.90 # p11, p22 ]) bounds [ (None, None), # mu (1e-6, None), (1e-6, None), # omega (1e-6, None), (1e-6, None), # alpha (0.0, None), (0.0, None), # beta (0.5, 0.999), (0.5, 0.999) # p ] constraints [ {type: ineq, fun: lambda p: 1 - (p[3] p[5])}, # alpha1 beta1 1 {type: ineq, fun: lambda p: 1 - (p[4] p[6])} # alpha2 beta2 1 ] res minimize( ms_garch_neg_loglik, init_params, args(r,), methodSLSQP, boundsbounds, constraintsconstraints, options{maxiter: 1000, ftol: 1e-10} ) return res res fit_ms_garch(ret.values) print(res.x)参数说明init_params 里 mu 从0开始ω 用样本方差的5%和20%起步——这个比例不对也不要紧SLSQP 对 ω 相对宽容。真正难的是转移概率和 α/β 之间的耦合所以我特意把 p11、p22 初值都放到0.9意思是“先假定状态很持久”再让优化器去调整。SLSQP 的 ftol 我习惯放到 1e-10因为优化目标是一个累积似然尺度大默认 1e-8 经常在离最优几步的地方就停住。maxiter 给1000威慑一下那些需要长时间跑的坏初值。跑完先看 res.success再看 res.x 与模拟真值是否接近。从实践角度看手写最大似然估计适合样本内拟合、快速原型、以及理解模型内部机制。如果你要的是贝叶斯后验区间或者样本量很大、状态数更多我建议调 R 的 MSGARCH 包——它内置了贝叶斯 MCMC还能输出平滑概率和后验区间。3.4 用 R 的 MSGARCH 包交叉验证结果R 语言这个包是当前 MS-GARCH 工程落地最省事的一条路径。代码短到几乎不用解释library(MSGARCH) ret - read.csv(ret.csv)$ret spec - CreateSpec(variance.spec list(model c(sGARCH, sGARCH)), distribution.spec list(distribution norm), switch.spec list(K 2)) set.seed(2024) fit - FitML(spec, data ret) summary(fit)CreateSpec 里 model 写两个 sGARCH表示两个区制都用标准 GARCHK2 指定区制数distributionnorm 和 Python 侧保持一致。FitML 是最大似然估计对应我上面的 SLSQP 过程。交叉验证不是拿两个人写的包互相对数值——最大似然本来就允许在浮点层面差一点。我的做法是看三件事方差参数的量级是否一致、转移概率是否同时指向同区间、平滑概率的区制划分时点是否基本重合。如果 R 包的平滑概率在某个时间点说“进入高波动”Python 手写版本说“还在低波动”那就说明滤波实现里某个细节有问题优先检查 Klaassen 近似和初始条件。4. 参数设定与诊断初值、约束、区制数选择与收敛检查4.1 一个能直接用的初始值模板真实数据上最怕的不是 garch 参数差一点而是转移概率初值给得离谱导致状态根本不切换。下面这组初值是我在日度金融收益率上最常用的默认值可以直接抄参数默认初值依据μ0.0日度收益率均值接近0ω1低波动0.05 * Var(r)低波动区制占样本方差的小头ω2高波动0.20 * Var(r)高波动区制承担主要方差α1、α20.05、0.20高波动区制对当日冲击更敏感β1、β20.90、0.70低波动区制衰减更慢但水平低p11、p220.90、0.90状态切换是低频事件这个模板对应“低波动区制承载大部分时间、高波动区制短暂但猛烈”的典型状态划分。如果你的数据是债券收益率或汇率ω2 和 α2 的比例要调低债券日度波动率变化比股票缓和得多。4.2 约束的必要性与实现细节αβ1 是 GARCH 弱平稳条件在 MS-GARCH 框架里应当按区制分别满足。注意我前面 Python 代码的 constraints 列表里下标没有写错p[3] 是 alpha1p[4] 是 alpha2p[5] 是 beta1p[6] 是 beta2因为参数向量是 [mu, omega1, omega2, alpha1, alpha2, beta1, beta2, p11, p22]。p11、p22 我约束在 [0.5, 0.999]。下界0.5是防止优化器把状态识别成“每天抛硬币切换”那种结果没有经济学含义上界0.999是防止 p 贴上1导致过滤概率变成0/1的退化序列。如果你用的是 R 的 FitML默认也内置了类似约束但默认上界是0.999999我建议调低到0.999数值上更稳。4.3 区制数选择先看 BIC再看经济含义区制数 K 的模型选择常见做法是拟合 K2 和 K3比较 BIC。注意这里有个 MS-GARCH 特有的现象因为每个区制自带一组 GARCH 参数K3 在样本上的似然改进往往显著但 BIC 的惩罚不一定拉得住结果出现第三区制只是把某段极端行情单独切出来换个样本区间就消失。我的选择顺序第一K2 和 K3 的 BIC 差第二K3 的平滑概率里是否有一个区制的平均持续期短于20个交易日——如果是说明那个区制不是“状态”只是几个离群值第三样本外预测方差是否真的优于 K2。如果 BIC 偏好 K3 但样本外预测没有改善我会放弃 K3。这个“预测验证优先于信息准则”的习惯帮我避开了好几次过度建模。4.4 收敛与识别诊断不要只盯着 log-likelihoodML 估计器报 success True 不等于模型可信。我每次拟合完都会做三个固定动作第一看看平滑概率filt[t]的分布。好的估计应该大部分时点接近0或1也就是说模型能清晰地区分状态。如果大量时点的概率在0.4到0.6之间徘徊说明两个区制的 GARCH 参数太接近或者数据根本不支持区制划分此时参数不稳定是必然的。第二重跑三组初值。把 ω2 初值改成0.1、0.3、0.5倍的样本方差分别拟合看是否落在同一个最优区域。这个操作能暴露似然面的多峰性。对于 MS-GARCH多峰几乎是常态原因就是后面要讲的标签切换。两峰之间隔着一个状态的互换log-likelihood 完全相同。第三检查转移概率的标准误。如果 p11 的标准误大于0.02或者估计值顶着0.999的上界说明这个状态的持续性没有被样本信息充分识别。原因可能是样本里高波动事件太少或者数据长度不够。5. 避坑排查MS-GARCH 的 5 个高发踩坑点5.1 标签切换同样的似然参数却“互换”了现象同一份数据换一个随机种子或初值估计结果里 ω1 和 ω2 互换α、β 也跟着互换log-likelihood 却几乎一样。我看似收敛了但区制1到底是低波动还是高波动变得不确定。原因两区制模型的似然函数关于状态编号是对称的。把“状态1”和“状态2”整体交换转移概率矩阵也对称交换似然不变。这是马尔科夫转换模型与生俱来的识别问题不是估计器写错了。解决最常见做法是在估计后做排序按 ω 的大小重命名区制ω 小的叫低波动、ω 大的叫高波动平滑概率跟着重排。如果不在意区制的具体编号这只影响展示不影响预测。但要注意如果你要做跨模型的对比比如两个市场分别估计后比较高波动区制参数排序必须统一。我在代码里加约束 ω1 ω2 也可以但 SLSQP 在边界处容易卡住所以我更倾向事后排序。5.2 路径依赖我们的似然本质上是一个近似现象Python 手写版本和 R 包对同一条序列的平滑概率高度一致但对数似然值差出一截或者你在某种极端参数下发现 h[t] 突然飙到离谱。原因严格似然需要对 2^T 条状态路径求和我们用的 Klaassen 近似是把条件方差对状态作期望等价于把“精确似然”换成了“近似似然”。在大多数参数下近似质量很好但当转移概率接近1、或者两个区制方差差异极大时近似的误差会被放大。解决承认近似然后做两件事。一是把 Klaassen 近似的结果和 R 包的 FitML 做对照R 包内部用的是它自己的数值近似偏差方向可以帮你判断二是如果模型用于高频风控参数估计后务必做滚动样本外验证近似导致的偏差最终会体现在预测方差的可预测性上。不要试图在论文层面和精确似然较劲工程场景下 Klaassen 是成本最低的后悔药。5.3 初始方差 h[0] 选择前几十个样本点经常“发疯”现象平滑概率图的前30到50个点出现一段不稳定的快速摆动之后才恢复正常。原因滤波从 h[0]Var(r) 出发这个初值对两个区制都不一定合理。高波动区制的真实方差可能远大于样本均值低波动区制又远小于它。滤波需要一段时间来“消化”错误的初始条件这段消化期的状态概率就会异常。解决最简单的做法是 burn-in从第 t100 个点开始累积似然前100个点只用于让滤波进入稳态。或者更讲究一点把 h[0] 当作额外参数估计但新增参数会加剧识别问题。我实际用的还是 burn-in代码里在 ll 累加时直接跳过前50或100个点。这也是 R 包不显式暴露这个问题的原因——它的内部实现已经处理了。5.4 正态分布假设低波动区制被负收益残差“拖垮”现象估计结果里低波动区制的 ω 偏大平滑概率却在下跌市里频繁跳到高波动但高波动区制的 α 和 β 看起来并不可靠。原因日度收益率有偏度和厚尾。用正态分布估计时所有负的大残差都要被高波动区制解释导致区制划分被“跌出来的波动”主导。特别是那些“阴跌”行情方差并不大但单边负收益让正态似然觉得很不舒服。解决把条件分布换成偏 tsstd或用稳健标准误重估。R 包里分布换成 sstdPython 手写代码则把 norm.pdf 换成 scipy.stats 的偏 t 密度。很多量化团队实际生产用的是偏 t你把它当默认也不过分。我的习惯是先纯正态跑通再切 sstd 看关键结论是否动摇——如果区制划分完全变了说明原来的风险判断不可靠。5.5 转移概率“顶到”边界状态持久性被高估现象p22 估计到0.999的上界平滑概率显示高波动区制一旦进入就长期存在但样本外预测方差偏大。原因高波动事件在样本中次数太少样本信息不足以精确估计退出概率。ML 估计趋向于把 p22 推到边界因为那会让训练似然稍微好一点。这个现象在只有一次极端段落的短样本里尤其明显。解决改用收缩先验或者在约束上手动封顶。贝叶斯思路是在 p22 上放 Beta 先验把均值拉到0.95附近比如 Beta(19, 1)等价于给5%的退出率先验。ML 思路则是在 SLSQP 的边界里把 p 的上界直接设为0.99当作工程约束。无论哪种都要认识到p22 点到上界附近时高波动区制的平均持续期是 1/(1-p22)实际只会越长越悲观风控上要把这个值当悲观假设用。6. 进阶用法把区制平滑概率变成能上线的风控信号前面五章解决的是“估出来”最后一章讲“怎么用”。我见过太多项目估计完 MS-GARCH 画个概率图就结束了其实这个模型最有价值的产出是区制平滑概率和一步向前预测方差。先说区制概率。滤波概率 filt[t] 给出的是 t 时刻基于全部历史信息的状态后验概率。工程上我会对它做平滑Kalman 后向递推得到 smooth probability。两者在样本内差别不大但平滑概率更稳。用高波动区制平滑概率超过0.5作为风险开关可以做一个很粗糙但有效的仓位控制器概率超过0.5就减半仓回到0.3以下再加回来。回测里这比固定阈值调节 GARCH 预测波动率更早反应因为状态切换是离散的。再做一步向前预测。t 时刻预测 t1 的方差需要对当前状态求期望# 假设 params 是已估计的参数filt_last 是终期滤波概率 def forecast_variance(params, r, filt_last, h_last, mu): omega params[1:3] alpha params[3:5] beta params[5:7] p11, p22 params[7], params[8] P np.array([[p11, 1-p11], [1-p22, p22]]) pred filt_last P var_forecast 0.0 for s in range(2): h_next_s omega[s] alpha[s] * (r[-1] - mu)**2 beta[s] * h_last[s] var_forecast pred[s] * h_next_s return var_forecast这个函数把“区制概率”和“区制内 GARCH 预测”做了加权。相比单一 GARCH这个预测在波动率跳升的第二天就能给更高的方差因为状态概率已经向高波动区制倾斜。最后说验证。我不信任任何没有对照的回归测试。每次接新的数据集都会让单一 GARCH 和 MS-GARCH 各出一份预测方差序列用 MSE 和分位数覆盖率做对比。MSE 当然要低但更关键的是分位数覆盖用正态分布的5%分位看实际收益率落在5%之外的频率是否接近5%。MS-GARCH 高波动区制的方差估计经常偏保守覆盖率偏离方向会直接反映模型假设是否该换成 sstd。做完这些结论才敢放进风控系统。我的习惯是每个季度用滚动窗口重估一次参数避免模型对陈旧危机产生条件反射。模型代码不复杂麻烦的永远是分布假设、初值、和状态编号这些细节。希望这些踩过的坑能够帮到你。本文还有配套的精品资源点击获取 马尔科夫转换MS-GARCH模型：从原理到Python实现与风控应用
简介面向金融时间序列分析的一份MATLAB实现资源主要围绕马尔科夫转换GARCH模型即MS-GARCH用于刻画经济扩张或收缩、市场情绪高低等状态下波动率的非线性切换。该资源适合金融计量经济学研究者、量化分析师及研究生作为算法学习与复现参考。压缩包共14个文件以12个MATLAB脚本为主分别用于似然函数计算、最大似然估计、协方差矩阵求解、起始值设定等环节另含1篇Markov-switching EGARCH的Matlab说明文档和1个文本格式的似然函数文件帮助对照理解。包体仅1.34MB内容紧凑便于快速上手。目前已有490人学习下载受到一定关注。从代码结构看内容覆盖MS-GARCH与MS-EGARCH两类实现包含QGML估计等关键步骤并附有教学型PDF。读者可借助脚本快速搭建波动率状态转换模型进行估计、模拟与预测也能在此基础上修改设定扩展自身研究整体实用性强。1. 先从“波动率突然放大”说起马尔科夫转换 MS-GARCH 解决了什么问题2015年夏天A股、2020年3月全球市场波动率不是缓慢爬升的是“跳”起来的。用单一GARCH(1,1)拟合这类序列αβ经常跑到0.98以上意思是任何一次冲击都需要上百个交易日才能消化这与实际感受完全不符。马尔科夫转换模型Markov Switching的思路很简单让数据存在两种或多种状态——低波动状态和高波动状态状态之间的切换由一个随时间演化的概率矩阵控制在每个状态内部再用GARCH刻画短期的波动聚集。把这个结构拼起来就是标题里的 MS-GARCHMarkov-Switching GARCH。这个模型不新鲜但在量化风控和波动率预测场景里一直很能打它能给出“当前处在哪个波动区制”的后验概率也能做比单一GARCH更贴近现实的多步方差预测。适合谁做波动率预测、期权定价、目标波动率调仓、以及需要区分“常态波动”和“危机波动”的风控从业者。2. MS-GARCH 的模型结构为什么是“马尔科夫转换 GARCH”而不是别家组合2.1 单一 GARCH 的 αβ 聚合问题波动聚集与结构突变被混在一起先说一个我经常拿来当开场白的现象。对沪深300日收益率直接拟合GARCH(1,1)在2015到2022这个区间上α估计值一般在0.06到0.1之间β在0.9到0.94附近αβ≈0.98甚至更高。这意味着条件方差的衰减半衰期大约在ln(0.5)/ln(0.98)≈35天一次大波动要五周才能衰减一半。但如果你去看波动率曲线的实际形状2015年6月那轮跳升之后波动率在两个月内就明显收敛了2016年初熔断前后又来了一波。这种反复“跳起来、落下去”的形态单一GARCH只能用一个高持久性参数去平均化拟合结果就是平时高估波动惯性危机后高估波动持续。Markov-Switching的想法不是修改GARCH方程而是把样本分成几个“波动区制”每个区制有自己的GARCH参数。低波动区制里αβ可能只有0.85冲击衰减快高波动区制里αβ可能接近0.99表示危机中波动聚集更强烈。这样就把“结构性变化”和“短期的波动聚集”分离开来。2.2 马尔科夫转换的数学形式S_t、转移概率矩阵 P、平稳分布MS-GARCH 的核心是一个不可观测的状态变量 S_t它只能取有限个离散值。两区制情况下 S_t1 表示低波动S_t2 表示高波动。状态的变化由一个齐次马尔科夫链控制P [p11, 1-p11; 1-p22, p22]其中 p11 Pr(S_t1 | S_{t-1}1)p22 Pr(S_t2 | S_{t-1}2)。对角元素衡量状态持续性p11 和 p22 越高系统越“恋战”一旦进入高波动区制就不容易出来。稳态概率 π1 (1-p22)/(2-p11-p22)π2 (1-p11)/(2-p11-p22)这个公式经常用来做初值设定。观测方程每个区制内仍是标准GARCH(1,1)σ_t²(s) ω_s α_s * ε_{t-1}² β_s * σ_{t-1}²(s)注意这里条件方差是带上标(s)的。问题在于t 时刻的 σ_t² 理论上依赖从 t0 到 t 的完整状态路径路径数量是 2^t这就是后面避坑章节要单独说的路径依赖问题。实现上必须做近似我在第3章会给出常用的近似估计方法。2.3 区制数K2 还是 K3大多数实证用两个区制就够了一个代表“平静期”一个代表“危机/高波动期”。两区制的优势是参数的局部识别相对容易转移概率和后验状态概率都比较好解释。三个区制也不是不能用它通常会把低波动进一步拆成“温和波动”和“极低波动”或者把高波动拆成“高波动”和“极端行情”。但K3的代价是GARCH参数从6个变成9个转移矩阵从2个参数变成6个样本量不足的时候估计结果对初值高度敏感而且容易出现两个高波动区制互相抢数据的翻车现场。我的经验是先用K2跑通全流程再拿AIC/BIC去看K3是否有实质改善如果BIC改善不到10个点就继续用K2。2.4 均值要不要也切换MS-GARCH 与 MS-AR-GARCH 的边界标题是MS-GARCH但落在“马尔科夫转换”这个大框架下均值方程的处理是个绕不开的选型点。最纯粹的形式是均值不切换μ 对所有区制相同区分完全靠波动率结构。好处是干净利落也是R语言MSGARCH包默认的做法。但实际收益率经常出现“高波动区制里均值更负”的现象比如2015年下半年、2018年波动率上行往往伴随负收益。这种情况下如果你只切波动率低波动区制的残差会被负收益拉偏导致方差估计失准。所以很多从业方案包括我在实际项目里的做法会这么选先跑纯 MS-GARCH看低波动区制残差均值是否显著不为零如果显著再升级到区制均值也切换的 MS-AR-GARCH。要注意的是均值切换会显著增加似然面的复杂度也让标签问题更严重所以不是越复杂越好。另外一个常被忽略的坑是条件分布。GARCH在险价值建模里正态分布假设会让高波动区制的尾部被低估。建议至少用 t 分布或偏 t 分布sstd做敏感性测试R的MSGARCH里可以直接指定distributionstd。我用正态只是为了让代码和讲解更直观。3. 用 Python 手写两区制 MS-GARCH(1,1)模拟数据、Hamilton 滤波、最大似然估计3.1 先造一份“答案已知”的数据模拟两区制 GARCH 序列模型没跑通之前先不要碰真实行情。我的习惯是先用已知参数模拟数据这样能确认估计器能不能恢复真值也方便后面排查。模拟步骤是从状态序列开始每个时点的条件方差由当期的状态决定然后按正态分布抽收益。import numpy as np import pandas as pd np.random.seed(42) n 2000 # 真实参数区制1低波动区制2高波动 mu_true 0.05 omega np.array([0.02, 0.50]) alpha np.array([0.05, 0.20]) beta np.array([0.90, 0.70]) p11_true, p22_true 0.98, 0.95 P_true np.array([[p11_true, 1 - p11_true], [1 - p22_true, p22_true]]) # 平稳分布采样初始状态 pi0 np.array([(1 - p22_true) / (2 - p11_true - p22_true), (1 - p11_true) / (2 - p11_true - p22_true)]) state np.zeros(n, dtypeint) state[0] np.random.choice([0, 1], ppi0) r np.zeros(n) h np.zeros(n) h[0] omega[state[0]] # 用一个非负初值 for t in range(n): if t 0: # 状态切换 state[t] np.random.choice([0, 1], pP_true[state[t-1]]) # 当前状态下的GARCH(1,1) h[t] omega[state[t]] alpha[state[t]] * r[t-1]**2 beta[state[t]] * h[t-1] r[t] mu_true np.sqrt(h[t]) * np.random.randn() ret pd.Series(r, nameret) print(ret.describe())逻辑说明这里用 r[t-1]² 而不是 (r[t-1]-mu)² 做模拟是为了让代码短一些实际估计时我用中心化残差。α_true(区制1)0.05 表示低波动区制内的杠杆效应很弱区制2的α0.20意味着高波动区制内“当日冲击”对次日方差的贡献大得多。β在不同区制的差异体现的就是波动衰减速度不同。实际项目里很少有人真的关心“真实参数”这个模拟环节主要用来检验估计器是否可识别如果模拟数据都估不准真实数据上更没法信。3.2 写似然函数用 Hamilton 滤波做状态概率递推用 Klaassen 近似处理路径依赖MS-GARCH 的似然难点不在马尔科夫链而在条件方差依赖状态历史。严格求和需要对 2^T 条路径做平均这在长度 2000 的序列上是天文数字。工程上的通行近似是 Klaassen2002的 collapsion 做法计算条件方差时用上一期每个状态下的滤波概率对 h[t-1] 做加权平均权重要乘上状态转移概率。from scipy.stats import norm def ms_garch_neg_loglik(params, r, n_states2): mu params[0] omega params[1:1n_states] alpha params[1n_states:12*n_states] beta params[12*n_states:13*n_states] p11, p22 params[-2], params[-1] P np.array([[p11, 1-p11], [1-p22, p22]]) T len(r) eps2 (r - mu) ** 2 # 初始状态概率P 的平稳分布 denom_ss 2 - p11 - p22 pi0 np.array([(1-p22)/denom_ss, (1-p11)/denom_ss]) filt np.zeros((T, n_states)) h np.zeros((T, n_states)) h[0] np.var(r) # 无条件样本方差作为共同起点 ll 0.0 for t in range(T): if t 0: dens norm.pdf(r[0], mu, np.sqrt(h[0])) num pi0 * dens ll np.log(num.sum()) filt[0] num / num.sum() else: # 一步预测状态概率Pr(S_t | F_{t-1}) pred filt[t-1] P # Klaassen 近似给定当前状态 s对上一期 h 做概率加权 for s in range(n_states): p_prev_given_s P[:, s] * filt[t-1] e_h (p_prev_given_s * h[t-1]).sum() / p_prev_given_s.sum() h[t, s] omega[s] alpha[s] * eps2[t-1] beta[s] * e_h # 当前状态下残差的分布密度 dens norm.pdf(r[t], mu, np.sqrt(h[t])) num pred * dens ll np.log(num.sum()) filt[t] num / num.sum() return -ll这段代码有几点要仔细说。第一为什么不是 h[t,s] ω_s α_s * ε_{t-1}² β_s * h[t-1,s]因为 t-1 时刻的真实方差并不是“某一个状态下的方差”而是对所有状态的不确定性加权平均。如果上一期滤波概率显示有30%概率处于高波动状态那么这一期低波动区制的方差起点也应该带上那30%的惯性。Klaassen 近似的本质就是把路径依赖用一个条件期望替代代价是似然不再是精确的但优化稳定性大幅提高。第二filt[t] 的更新是标准的 Hamilton 滤波先做一步预测 pred[t] filt[t-1] × P再乘以当前观测密度归一化后得到后验状态概率。这个概率序列本身就是重要的输出后面做区制识别和风控信号就是用它。第三初始状态概率用平稳分布。注意当 p11 或 p22 接近1时denom_ss 接近0pi0 会不稳定。我一般会加一个下限约束p11、p22 都落在 [0.5, 0.999] 区间内。3.3 用 SLSQP 做极大似然估计参数边界与 αβ1 约束优化时GARCH 参数必须满足 ω0、α≥0、β≥0、αβ1转移概率落在开区间(0,1)。这里 αβ 小于1是线性约束scipy.optimize.minimize 里用 SLSQP 同时支持 bounds 和 constraints。from scipy.optimize import minimize def fit_ms_garch(r, init_paramsNone): var0 np.var(r) if init_params is None: init_params np.array([ 0.0, # mu 0.05 * var0, 0.20 * var0, # omega1, omega2 0.05, 0.20, # alpha1, alpha2 0.90, 0.70, # beta1, beta2 0.90, 0.90 # p11, p22 ]) bounds [ (None, None), # mu (1e-6, None), (1e-6, None), # omega (1e-6, None), (1e-6, None), # alpha (0.0, None), (0.0, None), # beta (0.5, 0.999), (0.5, 0.999) # p ] constraints [ {type: ineq, fun: lambda p: 1 - (p[3] p[5])}, # alpha1 beta1 1 {type: ineq, fun: lambda p: 1 - (p[4] p[6])} # alpha2 beta2 1 ] res minimize( ms_garch_neg_loglik, init_params, args(r,), methodSLSQP, boundsbounds, constraintsconstraints, options{maxiter: 1000, ftol: 1e-10} ) return res res fit_ms_garch(ret.values) print(res.x)参数说明init_params 里 mu 从0开始ω 用样本方差的5%和20%起步——这个比例不对也不要紧SLSQP 对 ω 相对宽容。真正难的是转移概率和 α/β 之间的耦合所以我特意把 p11、p22 初值都放到0.9意思是“先假定状态很持久”再让优化器去调整。SLSQP 的 ftol 我习惯放到 1e-10因为优化目标是一个累积似然尺度大默认 1e-8 经常在离最优几步的地方就停住。maxiter 给1000威慑一下那些需要长时间跑的坏初值。跑完先看 res.success再看 res.x 与模拟真值是否接近。从实践角度看手写最大似然估计适合样本内拟合、快速原型、以及理解模型内部机制。如果你要的是贝叶斯后验区间或者样本量很大、状态数更多我建议调 R 的 MSGARCH 包——它内置了贝叶斯 MCMC还能输出平滑概率和后验区间。3.4 用 R 的 MSGARCH 包交叉验证结果R 语言这个包是当前 MS-GARCH 工程落地最省事的一条路径。代码短到几乎不用解释library(MSGARCH) ret - read.csv(ret.csv)$ret spec - CreateSpec(variance.spec list(model c(sGARCH, sGARCH)), distribution.spec list(distribution norm), switch.spec list(K 2)) set.seed(2024) fit - FitML(spec, data ret) summary(fit)CreateSpec 里 model 写两个 sGARCH表示两个区制都用标准 GARCHK2 指定区制数distributionnorm 和 Python 侧保持一致。FitML 是最大似然估计对应我上面的 SLSQP 过程。交叉验证不是拿两个人写的包互相对数值——最大似然本来就允许在浮点层面差一点。我的做法是看三件事方差参数的量级是否一致、转移概率是否同时指向同区间、平滑概率的区制划分时点是否基本重合。如果 R 包的平滑概率在某个时间点说“进入高波动”Python 手写版本说“还在低波动”那就说明滤波实现里某个细节有问题优先检查 Klaassen 近似和初始条件。4. 参数设定与诊断初值、约束、区制数选择与收敛检查4.1 一个能直接用的初始值模板真实数据上最怕的不是 garch 参数差一点而是转移概率初值给得离谱导致状态根本不切换。下面这组初值是我在日度金融收益率上最常用的默认值可以直接抄参数默认初值依据μ0.0日度收益率均值接近0ω1低波动0.05 * Var(r)低波动区制占样本方差的小头ω2高波动0.20 * Var(r)高波动区制承担主要方差α1、α20.05、0.20高波动区制对当日冲击更敏感β1、β20.90、0.70低波动区制衰减更慢但水平低p11、p220.90、0.90状态切换是低频事件这个模板对应“低波动区制承载大部分时间、高波动区制短暂但猛烈”的典型状态划分。如果你的数据是债券收益率或汇率ω2 和 α2 的比例要调低债券日度波动率变化比股票缓和得多。4.2 约束的必要性与实现细节αβ1 是 GARCH 弱平稳条件在 MS-GARCH 框架里应当按区制分别满足。注意我前面 Python 代码的 constraints 列表里下标没有写错p[3] 是 alpha1p[4] 是 alpha2p[5] 是 beta1p[6] 是 beta2因为参数向量是 [mu, omega1, omega2, alpha1, alpha2, beta1, beta2, p11, p22]。p11、p22 我约束在 [0.5, 0.999]。下界0.5是防止优化器把状态识别成“每天抛硬币切换”那种结果没有经济学含义上界0.999是防止 p 贴上1导致过滤概率变成0/1的退化序列。如果你用的是 R 的 FitML默认也内置了类似约束但默认上界是0.999999我建议调低到0.999数值上更稳。4.3 区制数选择先看 BIC再看经济含义区制数 K 的模型选择常见做法是拟合 K2 和 K3比较 BIC。注意这里有个 MS-GARCH 特有的现象因为每个区制自带一组 GARCH 参数K3 在样本上的似然改进往往显著但 BIC 的惩罚不一定拉得住结果出现第三区制只是把某段极端行情单独切出来换个样本区间就消失。我的选择顺序第一K2 和 K3 的 BIC 差第二K3 的平滑概率里是否有一个区制的平均持续期短于20个交易日——如果是说明那个区制不是“状态”只是几个离群值第三样本外预测方差是否真的优于 K2。如果 BIC 偏好 K3 但样本外预测没有改善我会放弃 K3。这个“预测验证优先于信息准则”的习惯帮我避开了好几次过度建模。4.4 收敛与识别诊断不要只盯着 log-likelihoodML 估计器报 success True 不等于模型可信。我每次拟合完都会做三个固定动作第一看看平滑概率filt[t]的分布。好的估计应该大部分时点接近0或1也就是说模型能清晰地区分状态。如果大量时点的概率在0.4到0.6之间徘徊说明两个区制的 GARCH 参数太接近或者数据根本不支持区制划分此时参数不稳定是必然的。第二重跑三组初值。把 ω2 初值改成0.1、0.3、0.5倍的样本方差分别拟合看是否落在同一个最优区域。这个操作能暴露似然面的多峰性。对于 MS-GARCH多峰几乎是常态原因就是后面要讲的标签切换。两峰之间隔着一个状态的互换log-likelihood 完全相同。第三检查转移概率的标准误。如果 p11 的标准误大于0.02或者估计值顶着0.999的上界说明这个状态的持续性没有被样本信息充分识别。原因可能是样本里高波动事件太少或者数据长度不够。5. 避坑排查MS-GARCH 的 5 个高发踩坑点5.1 标签切换同样的似然参数却“互换”了现象同一份数据换一个随机种子或初值估计结果里 ω1 和 ω2 互换α、β 也跟着互换log-likelihood 却几乎一样。我看似收敛了但区制1到底是低波动还是高波动变得不确定。原因两区制模型的似然函数关于状态编号是对称的。把“状态1”和“状态2”整体交换转移概率矩阵也对称交换似然不变。这是马尔科夫转换模型与生俱来的识别问题不是估计器写错了。解决最常见做法是在估计后做排序按 ω 的大小重命名区制ω 小的叫低波动、ω 大的叫高波动平滑概率跟着重排。如果不在意区制的具体编号这只影响展示不影响预测。但要注意如果你要做跨模型的对比比如两个市场分别估计后比较高波动区制参数排序必须统一。我在代码里加约束 ω1 ω2 也可以但 SLSQP 在边界处容易卡住所以我更倾向事后排序。5.2 路径依赖我们的似然本质上是一个近似现象Python 手写版本和 R 包对同一条序列的平滑概率高度一致但对数似然值差出一截或者你在某种极端参数下发现 h[t] 突然飙到离谱。原因严格似然需要对 2^T 条状态路径求和我们用的 Klaassen 近似是把条件方差对状态作期望等价于把“精确似然”换成了“近似似然”。在大多数参数下近似质量很好但当转移概率接近1、或者两个区制方差差异极大时近似的误差会被放大。解决承认近似然后做两件事。一是把 Klaassen 近似的结果和 R 包的 FitML 做对照R 包内部用的是它自己的数值近似偏差方向可以帮你判断二是如果模型用于高频风控参数估计后务必做滚动样本外验证近似导致的偏差最终会体现在预测方差的可预测性上。不要试图在论文层面和精确似然较劲工程场景下 Klaassen 是成本最低的后悔药。5.3 初始方差 h[0] 选择前几十个样本点经常“发疯”现象平滑概率图的前30到50个点出现一段不稳定的快速摆动之后才恢复正常。原因滤波从 h[0]Var(r) 出发这个初值对两个区制都不一定合理。高波动区制的真实方差可能远大于样本均值低波动区制又远小于它。滤波需要一段时间来“消化”错误的初始条件这段消化期的状态概率就会异常。解决最简单的做法是 burn-in从第 t100 个点开始累积似然前100个点只用于让滤波进入稳态。或者更讲究一点把 h[0] 当作额外参数估计但新增参数会加剧识别问题。我实际用的还是 burn-in代码里在 ll 累加时直接跳过前50或100个点。这也是 R 包不显式暴露这个问题的原因——它的内部实现已经处理了。5.4 正态分布假设低波动区制被负收益残差“拖垮”现象估计结果里低波动区制的 ω 偏大平滑概率却在下跌市里频繁跳到高波动但高波动区制的 α 和 β 看起来并不可靠。原因日度收益率有偏度和厚尾。用正态分布估计时所有负的大残差都要被高波动区制解释导致区制划分被“跌出来的波动”主导。特别是那些“阴跌”行情方差并不大但单边负收益让正态似然觉得很不舒服。解决把条件分布换成偏 tsstd或用稳健标准误重估。R 包里分布换成 sstdPython 手写代码则把 norm.pdf 换成 scipy.stats 的偏 t 密度。很多量化团队实际生产用的是偏 t你把它当默认也不过分。我的习惯是先纯正态跑通再切 sstd 看关键结论是否动摇——如果区制划分完全变了说明原来的风险判断不可靠。5.5 转移概率“顶到”边界状态持久性被高估现象p22 估计到0.999的上界平滑概率显示高波动区制一旦进入就长期存在但样本外预测方差偏大。原因高波动事件在样本中次数太少样本信息不足以精确估计退出概率。ML 估计趋向于把 p22 推到边界因为那会让训练似然稍微好一点。这个现象在只有一次极端段落的短样本里尤其明显。解决改用收缩先验或者在约束上手动封顶。贝叶斯思路是在 p22 上放 Beta 先验把均值拉到0.95附近比如 Beta(19, 1)等价于给5%的退出率先验。ML 思路则是在 SLSQP 的边界里把 p 的上界直接设为0.99当作工程约束。无论哪种都要认识到p22 点到上界附近时高波动区制的平均持续期是 1/(1-p22)实际只会越长越悲观风控上要把这个值当悲观假设用。6. 进阶用法把区制平滑概率变成能上线的风控信号前面五章解决的是“估出来”最后一章讲“怎么用”。我见过太多项目估计完 MS-GARCH 画个概率图就结束了其实这个模型最有价值的产出是区制平滑概率和一步向前预测方差。先说区制概率。滤波概率 filt[t] 给出的是 t 时刻基于全部历史信息的状态后验概率。工程上我会对它做平滑Kalman 后向递推得到 smooth probability。两者在样本内差别不大但平滑概率更稳。用高波动区制平滑概率超过0.5作为风险开关可以做一个很粗糙但有效的仓位控制器概率超过0.5就减半仓回到0.3以下再加回来。回测里这比固定阈值调节 GARCH 预测波动率更早反应因为状态切换是离散的。再做一步向前预测。t 时刻预测 t1 的方差需要对当前状态求期望# 假设 params 是已估计的参数filt_last 是终期滤波概率 def forecast_variance(params, r, filt_last, h_last, mu): omega params[1:3] alpha params[3:5] beta params[5:7] p11, p22 params[7], params[8] P np.array([[p11, 1-p11], [1-p22, p22]]) pred filt_last P var_forecast 0.0 for s in range(2): h_next_s omega[s] alpha[s] * (r[-1] - mu)**2 beta[s] * h_last[s] var_forecast pred[s] * h_next_s return var_forecast这个函数把“区制概率”和“区制内 GARCH 预测”做了加权。相比单一 GARCH这个预测在波动率跳升的第二天就能给更高的方差因为状态概率已经向高波动区制倾斜。最后说验证。我不信任任何没有对照的回归测试。每次接新的数据集都会让单一 GARCH 和 MS-GARCH 各出一份预测方差序列用 MSE 和分位数覆盖率做对比。MSE 当然要低但更关键的是分位数覆盖用正态分布的5%分位看实际收益率落在5%之外的频率是否接近5%。MS-GARCH 高波动区制的方差估计经常偏保守覆盖率偏离方向会直接反映模型假设是否该换成 sstd。做完这些结论才敢放进风控系统。我的习惯是每个季度用滚动窗口重估一次参数避免模型对陈旧危机产生条件反射。模型代码不复杂麻烦的永远是分布假设、初值、和状态编号这些细节。希望这些踩过的坑能够帮到你。本文还有配套的精品资源点击获取