拓冰网站定制 
简介面向金融时间序列分析与量化研究者这份MATLAB实现包聚焦马尔科夫转换GARCHMS-GARCH模型覆盖数据预处理、状态定义、GARCH参数设定、最大似然估计与模型诊断等完整建模流程适合已有GARCH基础、希望掌握状态切换波动率建模的读者。压缩包共14个文件大小1.34MB其中12个m脚本为主体包含主程序、似然函数、协方差矩阵等核心代码模块并附1份Markov-switching EGARCH的PDF说明文档及1个txt文本文件配合PDF可快速理解算法推导、参数初始化与转换概率矩阵估计等关键步骤。代码模块划分清晰涵盖普通GARCH、EGARCH及QML估计等多种变体可直接运行或改造用于股市指数、汇率等收益率序列的状态识别与波动预测。已有490人学习下载对需要实证复现MS-GARCH/马尔科夫转换模型的科研与实务人员具有直接参考价值。1. MS-GARCH 在解决什么问题一段收益率序列里可能住着好几个波动率做波动率建模的人早晚会遇到一个尴尬拿 GARCH(1,1) 拟合沪深300 的日收益平稳年份效果不错一到 2015 年、2020 年这种极端行情模型残差直接爆表预测区间看起来像笑话。问题往往不在 GARCH 本身而是波动率的“脾气”会换档。MS-GARCH马尔可夫转换 GARCH也叫 Markov-Switching GARCH中文常写作“马尔科夫转换”或“马尔可夫机制转换”就是让模型自己在低波动和高波动两种或多种状态之间切换每个状态配一套独立的 GARCH 参数同时估计状态转移概率。它解决的痛点是传统 GARCH 用一套参数描述所有时期而实际市场波动并不稳结构会变化。适合读这篇文章的人是做波动率预测、VaR 和压力测试、期权隐含波动率研究的从业者也包括用量价序列做状态择时的量化研究员。下面从模型结构、R 语言实操、Python 替代路线一直讲到高频踩坑现场最后收在样本外预测的验证方法上。2. 马尔可夫链和 GARCH 方程怎么耦合先把模型结构看明白2.1 为什么用马尔可夫转换而不是人肉断点或者滚动窗口很多人面对波动率结构变化的第一个想法是找断点。常见做法是用 Chow 检验或者 Bai-Perron 检验把样本切成几段每段单独拟合 GARCH。这个方法在学术上成熟但落地有一个硬伤断点是外生给定的你得先知道 2015 年股灾和 2020 年疫情是断点样本外再出现新断点时模型完全来不及反应。滚动窗口则是另一个思路窗口宽度本身就是个玄学窗口短了参数噪声大长了又对结构变化反应迟钝。马尔可夫转换的做法不一样。它把状态变量 s_t 设成一个不可观测的马尔可夫链状态之间的转移概率是模型内生的也就是说模型自己从数据里学“现在更像哪个波动状态”而不是由人事先指定。这个特性对预测特别重要你不需要提前知道危机哪天来只需要知道当前处于高波动状态的概率在上升就可以提前调整风险敞口。相比之下断点法的“后悔药”属性很强只适合事后归因不适合实时决策。2.2 两状态 MS-GARCH 的基本方程与参数MS-GARCH 的完整形式可以写成一个带状态条件均值和状态条件方差的方程组。以最常见的两状态 MS-GARCH(1,1) 为例观测方程y_t μ_{s_t} ε_tε_t σ_{t,s_t} · z_tz_t ~ N(0, 1)方差方程每个状态各一套σ_t² ω_{s_t} α_{s_t} · ε_{t-1}² β_{s_t} · σ_{t-1}²其中 s_t ∈ {1, 2}状态转移概率矩阵是P [[p, 1-p], [1-q, q]]p 表示状态 1 下一期仍停留在状态 1 的概率q 表示状态 2 下期仍停留在状态 2 的概率。状态持续期的期望是 1/(1-p) 和 1/(1-q)。这里有一个新手最容易绕进去的点σ_{t-1}² 应该取哪个状态的方差因为在 t-1 时刻我们并不知道当时处于哪个状态模型实际上需要维护两条方差路径分别对应“当前是状态 1”和“当前是状态 2”。在似然函数里每个时刻的状态概率通过 Hamilton 滤波递归更新然后把两个状态的密度按概率加权。这个机制学术上叫 Gray 折叠Gray, 1996或 Haas 的独立路径法。落地时我们一般用 Gray 折叠后面第四章会给出 Python 代码骨架。参数说明一下每个状态有一组 (μ, ω, α, β)加上转移概率 p、q一共是 2×4 2 10 个待估参数。样本量少于 500 时10 个参数很容易过拟合所以实际使用中常用 μ_1 μ_2 的约束版本只让方差参数跨状态变化。2.3 MS-GARCH 的识别难点为什么它比普通 GARCH 更难估计普通 GARCH 的似然函数通常是单峰的用优化器随便给个初值也能收敛。MS-GARCH 不一样似然函数是多峰的而且存在标签切换问题label switching把状态 1 和状态 2 的参数整体互换似然值完全一样优化器可能收敛到两套等价解上让你分不清哪个是低波动哪个是高波动。另一个识别难点是状态退化如果某个状态只出现过几次比如 2020 年 2 月那种持续十几天的高波动模型会把转移概率拉到接近 1状态持续时间估计方差极大参数置信区间宽得没有实用价值。还有一个数学上必须注意的地方状态条件方差方程里ε_{t-1}² 与状态 s_t 之间的关联不是标准的 GARCH 平滑结构所以在估计时不能直接套用普通 GARCH 的方差目标化variance targeting技巧。很多人在 R 包MSGARCH里选了 variance.target TRUE结果发现状态参数估计非常不稳定。正确的做法是先分别以全样本方差作为每个状态的未条件方差初值再让优化器自由搜索。3. R 语言实战用 MSGARCH 包跑通一次完整建模3.1 拟合前的数据处理与平稳性预处理用 MS-GARCH 之前数据必须是平稳序列。常规做法是取对数收益率r_t ln(P_t/P_{t-1}) × 100乘以 100 是为了让参数量级更友好。对 A 股日线数据我一般要求至少 1000 个交易日也就是四年左右。少于这个量级状态转移概率的估计误差会非常大高波动状态可能只被观测到十几个样本点参数不可信。数据预处理有一个常见坑不要对收益率再做差分也不要做去均值。MS-GARCH 的条件均值方程会自己估计 μ你提前减掉均值会导致似然函数里的 μ 失去识别力状态概率估计会偏移。如果原始价格序列有明显趋势对数收益已经基本消除了。异常值处理方面不建议用 winsorize 截尾因为极端收益正是高波动状态的重要信息如果存在明显数据错误比如价格跌为 0 或负值直接删除对应行即可。预处理代码用 R 可以这样写library(quantmod) # 以某指数日线为例实际使用时可以换成自己的行情源 getSymbols(000300.SS, from 2016-01-01, to 2024-12-31, auto.assign TRUE) px - Cl(000300.SS) # 不渲染价格序列 ret - as.numeric(diff(log(px)) * 100) ret - ret[is.finite(ret)] # 去掉 NA 和异常值 plot(ret, type l, main 日对数收益率(%))这段代码的逻辑是先用diff(log(px))计算对数收益率并乘 100然后过滤掉非有限值。很多人忽略is.finite这一步导致后面优化器在遇到 NA 时报错或者收敛到奇怪的位置。另外要留意getSymbols获取国内指数可能需要额外设置数据源正式落地时我一般直接用本地数据库或 Wind 导出的 CSV代码只是展示数据处理流程。3.2 最小可复现拟合代码与输出解读R 生态里最常用的包是MSGARCH它把参数估计、状态概率滤波和预测都封装好了。拟合一个两状态、正态分布、GARCH(1,1) 的 MS-GARCH 模型最小代码是library(MSGARCH) # 创建模型规格两状态、标准 GARCH(1,1)、正态扰动 spec - CreateSpec( variance.spec list(model sGARCH), distribution.spec list(distribution norm), switch.spec list(K 2) ) # 极大似然估计 fit - FitML(spec, series ret) summary(fit)CreateSpec的参数含义variance.spec list(model sGARCH)指定每个状态内部用的是标准 GARCH(1,1) 方差方程distribution.spec控制扰动项分布可选norm、std学生 t、sstd偏 t等switch.spec list(K 2)就是马尔可夫转换的状态个数。对于日频收益率我一般直接用std而不是norm因为收益尾部更厚但厚尾分布会显著增加优化时间如果只是先跑通流程用norm更合适。FitML返回的 summary 里重点看三类内容。第一是每个状态的参数估计包括 ω_1、ω_2、α_1、α_2、β_1、β_2注意对比两个状态的 αβ 之和它反映波动率衰减速度高波动状态的 αβ 通常更接近 1说明冲击持续更久。第二是转移概率矩阵 p、q如果 p 估计为 0.99 以上说明状态 1 非常持久。第三是持续期期望1/(1-p) 和 1/(1-q) 直接用笔算就行。3.3 提取状态概率与绘制平滑概率图参数估计完下一步就是把每个时刻属于状态 1 和状态 2 的概率提取出来。MSGARCH包提供了StateProb和PosteriorProb两个函数前者返回滤波概率后者返回平滑概率。平滑概率用了全样本信息事后看更准实时决策只能用滤波概率。# 平滑概率每个时点位于各状态的后验概率 sp - PosteriorProb(fit) matplot(sp, type l, lty 1, col c(steelblue, firebrick), xlab t, ylab Probability, main Smoothed State Probabilities) legend(topright, legend c(State 1 (Low Vol), State 2 (High Vol)), col c(steelblue, firebrick), lty 1)画图后通常能明显看到高波动状态概率在 2015 年、2020 年、2022 年出现尖峰状抬升。如果画出来两个状态概率一直在 0.5 上下纠缠不清说明模型没有把两状态区分开后续使用要谨慎。这是模型失败的最常见信号之一原因可能出在初始值或者状态数选择上。处理办法放在第五章专门讲。4. Python 里的替代路线自写似然函数复现两状态 MS-GARCH4.1 为什么选择从零写似然而不是直接找现成包装很多人希望像 R 一样pip install ms_garch一把梭。现实是 Python 生态里没有一个统一维护、开箱即用的 MS-GARCH 包。有的库叫markovregime或者hmmlearn但前者只支持均值转换不支持方差方程里的 GARCH 结构后者是通用隐马尔可夫模型观测分布必须是同一族分布不能天然和 GARCH 递归耦合。因此实际工作中Python 落地 MS-GARCH 有两条路一是通过rpy2直接把 R 的MSGARCH包拉过来跑适合不想重复造轮子的场景二是自己写负对数似然函数配合scipy.optimize或者statsmodels的优化器做估计。本节给出第二种方案的最小骨架它还能帮你把第二章没想清楚的 Hamilton 滤波和 Gray 折叠落到实处。4.2 负对数似然与状态概率滤波的代码骨架两状态 MS-GARCH 的似然计算核心是 Hamilton 滤波。每个时点 t 迭代四步计算每个状态的观测密度、计算混合密度、更新滤波概率、预测下一期状态概率。方差递归使用 Gray 折叠即在每个时点用滤波概率加权的混合方差作为下一期方差递归的输入。这样处理虽然会损失一部分理论上的一致性但工程上稳定不容易发散是实证中最常见的做法。import numpy as np from scipy.stats import norm from scipy.optimize import minimize def ms_garch_nll(params, r): # 参数顺序: [mu1, mu2, omega1, omega2, alpha1, alpha2, beta1, beta2, p_logit, q_logit] mu params[0:2] omega np.exp(params[2:4]) # 保证方差为正 alpha np.exp(params[4:6]) # 保证非负 beta np.exp(params[6:8]) # 保证非负 p 1.0 / (1.0 np.exp(-params[8])) # 反logit到(0,1) q 1.0 / (1.0 np.exp(-params[9])) T len(r) xi np.array([0.5, 0.5]) # 初始状态概率取均匀 sigma2 np.array([np.var(r), np.var(r)]) # 两条状态方差路径 loglik 0.0 for t in range(T): # 1. 每个状态在当前值下的正态密度 eta norm.pdf(r[t], locmu, scalenp.sqrt(sigma2)) # 2. 混合密度用于累积似然 dens np.dot(xi, eta) if dens 1e-12: dens 1e-12 loglik np.log(dens) # 3. 滤波后验概率 xi_f xi * eta / dens # 4. 预测下一期状态概率 xi np.array([p * xi_f[0] (1 - q) * xi_f[1], (1 - p) * xi_f[0] q * xi_f[1]]) # 5. Gray折叠更新方差 mixed_var np.dot(xi_f, sigma2) sigma2 omega alpha * (r[t] ** 2) beta * mixed_var return -loglik这段代码有几个必须说明的参数细节。第一ω、α、β 用exp参数化是为了让优化器在无约束空间里工作同时保证估计值非负如果你直接用原始参数scipy.optimize很容易在迭代过程中把方差推成负数报错或者收敛到 NaN。第二p 和 q 用反 logit 变换映射到 (0,1)避免转移概率越界。第三初始状态概率取 {0.5, 0.5} 并且不加约束对结果影响不大但如果样本里两个状态占比极不均衡用 {0.9, 0.1} 之类的先验值有时能帮助优化器更快找到好区域。第四方差初值取全样本方差比取 0 或者取第一个观测本身的残差平方更平稳能显著减少迭代初期的数值振荡。4.3 用 scipy.optimize 做数值优化的参数陷阱有了目标函数后调用优化器的写法如下# 初始化均值取样本均值附近方差参数取小正数p/q 取 0.9 init [0.0, 0.0, np.log(np.var(ret)/2), np.log(np.var(ret)*2), np.log(0.05), np.log(0.1), np.log(0.9), np.log(0.85), 2.0, 2.0] # logit(0.9)≈2.197 res minimize(ms_garch_nll, init, args(ret,), methodL-BFGS-B, options{maxiter: 2000}) print(res.x, res.fun)这里最容易翻车的点是scipy.optimize默认的迭代次数远不够 MS-GARCH 这种高维非凸问题收敛建议至少给maxiter2000。即便如此一次优化也经常落进局部最优。我常用的策略是并行做 20 到 50 组随机初值取负对数似然最小的那一组作为最终解。具体操作时用numpy.random.uniform给均值、方差参数和转移概率各生成一组初值然后把每组初值代入优化器最后对比res.fun。这个多起点搜索在 Python 里跑起来很快因为每个序列长度也就几百到几千单次优化耗时通常不到一秒到几秒。另一个容易忽略的问题是收敛标志。res.success为True只代表优化器行程结束不代表找到全局最优。更加可本文还有配套的精品资源点击获取 MS-GARCH解析：马尔可夫转换下波动率状态切换建模与实战
简介面向金融时间序列分析与量化研究者这份MATLAB实现包聚焦马尔科夫转换GARCHMS-GARCH模型覆盖数据预处理、状态定义、GARCH参数设定、最大似然估计与模型诊断等完整建模流程适合已有GARCH基础、希望掌握状态切换波动率建模的读者。压缩包共14个文件大小1.34MB其中12个m脚本为主体包含主程序、似然函数、协方差矩阵等核心代码模块并附1份Markov-switching EGARCH的PDF说明文档及1个txt文本文件配合PDF可快速理解算法推导、参数初始化与转换概率矩阵估计等关键步骤。代码模块划分清晰涵盖普通GARCH、EGARCH及QML估计等多种变体可直接运行或改造用于股市指数、汇率等收益率序列的状态识别与波动预测。已有490人学习下载对需要实证复现MS-GARCH/马尔科夫转换模型的科研与实务人员具有直接参考价值。1. MS-GARCH 在解决什么问题一段收益率序列里可能住着好几个波动率做波动率建模的人早晚会遇到一个尴尬拿 GARCH(1,1) 拟合沪深300 的日收益平稳年份效果不错一到 2015 年、2020 年这种极端行情模型残差直接爆表预测区间看起来像笑话。问题往往不在 GARCH 本身而是波动率的“脾气”会换档。MS-GARCH马尔可夫转换 GARCH也叫 Markov-Switching GARCH中文常写作“马尔科夫转换”或“马尔可夫机制转换”就是让模型自己在低波动和高波动两种或多种状态之间切换每个状态配一套独立的 GARCH 参数同时估计状态转移概率。它解决的痛点是传统 GARCH 用一套参数描述所有时期而实际市场波动并不稳结构会变化。适合读这篇文章的人是做波动率预测、VaR 和压力测试、期权隐含波动率研究的从业者也包括用量价序列做状态择时的量化研究员。下面从模型结构、R 语言实操、Python 替代路线一直讲到高频踩坑现场最后收在样本外预测的验证方法上。2. 马尔可夫链和 GARCH 方程怎么耦合先把模型结构看明白2.1 为什么用马尔可夫转换而不是人肉断点或者滚动窗口很多人面对波动率结构变化的第一个想法是找断点。常见做法是用 Chow 检验或者 Bai-Perron 检验把样本切成几段每段单独拟合 GARCH。这个方法在学术上成熟但落地有一个硬伤断点是外生给定的你得先知道 2015 年股灾和 2020 年疫情是断点样本外再出现新断点时模型完全来不及反应。滚动窗口则是另一个思路窗口宽度本身就是个玄学窗口短了参数噪声大长了又对结构变化反应迟钝。马尔可夫转换的做法不一样。它把状态变量 s_t 设成一个不可观测的马尔可夫链状态之间的转移概率是模型内生的也就是说模型自己从数据里学“现在更像哪个波动状态”而不是由人事先指定。这个特性对预测特别重要你不需要提前知道危机哪天来只需要知道当前处于高波动状态的概率在上升就可以提前调整风险敞口。相比之下断点法的“后悔药”属性很强只适合事后归因不适合实时决策。2.2 两状态 MS-GARCH 的基本方程与参数MS-GARCH 的完整形式可以写成一个带状态条件均值和状态条件方差的方程组。以最常见的两状态 MS-GARCH(1,1) 为例观测方程y_t μ_{s_t} ε_tε_t σ_{t,s_t} · z_tz_t ~ N(0, 1)方差方程每个状态各一套σ_t² ω_{s_t} α_{s_t} · ε_{t-1}² β_{s_t} · σ_{t-1}²其中 s_t ∈ {1, 2}状态转移概率矩阵是P [[p, 1-p], [1-q, q]]p 表示状态 1 下一期仍停留在状态 1 的概率q 表示状态 2 下期仍停留在状态 2 的概率。状态持续期的期望是 1/(1-p) 和 1/(1-q)。这里有一个新手最容易绕进去的点σ_{t-1}² 应该取哪个状态的方差因为在 t-1 时刻我们并不知道当时处于哪个状态模型实际上需要维护两条方差路径分别对应“当前是状态 1”和“当前是状态 2”。在似然函数里每个时刻的状态概率通过 Hamilton 滤波递归更新然后把两个状态的密度按概率加权。这个机制学术上叫 Gray 折叠Gray, 1996或 Haas 的独立路径法。落地时我们一般用 Gray 折叠后面第四章会给出 Python 代码骨架。参数说明一下每个状态有一组 (μ, ω, α, β)加上转移概率 p、q一共是 2×4 2 10 个待估参数。样本量少于 500 时10 个参数很容易过拟合所以实际使用中常用 μ_1 μ_2 的约束版本只让方差参数跨状态变化。2.3 MS-GARCH 的识别难点为什么它比普通 GARCH 更难估计普通 GARCH 的似然函数通常是单峰的用优化器随便给个初值也能收敛。MS-GARCH 不一样似然函数是多峰的而且存在标签切换问题label switching把状态 1 和状态 2 的参数整体互换似然值完全一样优化器可能收敛到两套等价解上让你分不清哪个是低波动哪个是高波动。另一个识别难点是状态退化如果某个状态只出现过几次比如 2020 年 2 月那种持续十几天的高波动模型会把转移概率拉到接近 1状态持续时间估计方差极大参数置信区间宽得没有实用价值。还有一个数学上必须注意的地方状态条件方差方程里ε_{t-1}² 与状态 s_t 之间的关联不是标准的 GARCH 平滑结构所以在估计时不能直接套用普通 GARCH 的方差目标化variance targeting技巧。很多人在 R 包MSGARCH里选了 variance.target TRUE结果发现状态参数估计非常不稳定。正确的做法是先分别以全样本方差作为每个状态的未条件方差初值再让优化器自由搜索。3. R 语言实战用 MSGARCH 包跑通一次完整建模3.1 拟合前的数据处理与平稳性预处理用 MS-GARCH 之前数据必须是平稳序列。常规做法是取对数收益率r_t ln(P_t/P_{t-1}) × 100乘以 100 是为了让参数量级更友好。对 A 股日线数据我一般要求至少 1000 个交易日也就是四年左右。少于这个量级状态转移概率的估计误差会非常大高波动状态可能只被观测到十几个样本点参数不可信。数据预处理有一个常见坑不要对收益率再做差分也不要做去均值。MS-GARCH 的条件均值方程会自己估计 μ你提前减掉均值会导致似然函数里的 μ 失去识别力状态概率估计会偏移。如果原始价格序列有明显趋势对数收益已经基本消除了。异常值处理方面不建议用 winsorize 截尾因为极端收益正是高波动状态的重要信息如果存在明显数据错误比如价格跌为 0 或负值直接删除对应行即可。预处理代码用 R 可以这样写library(quantmod) # 以某指数日线为例实际使用时可以换成自己的行情源 getSymbols(000300.SS, from 2016-01-01, to 2024-12-31, auto.assign TRUE) px - Cl(000300.SS) # 不渲染价格序列 ret - as.numeric(diff(log(px)) * 100) ret - ret[is.finite(ret)] # 去掉 NA 和异常值 plot(ret, type l, main 日对数收益率(%))这段代码的逻辑是先用diff(log(px))计算对数收益率并乘 100然后过滤掉非有限值。很多人忽略is.finite这一步导致后面优化器在遇到 NA 时报错或者收敛到奇怪的位置。另外要留意getSymbols获取国内指数可能需要额外设置数据源正式落地时我一般直接用本地数据库或 Wind 导出的 CSV代码只是展示数据处理流程。3.2 最小可复现拟合代码与输出解读R 生态里最常用的包是MSGARCH它把参数估计、状态概率滤波和预测都封装好了。拟合一个两状态、正态分布、GARCH(1,1) 的 MS-GARCH 模型最小代码是library(MSGARCH) # 创建模型规格两状态、标准 GARCH(1,1)、正态扰动 spec - CreateSpec( variance.spec list(model sGARCH), distribution.spec list(distribution norm), switch.spec list(K 2) ) # 极大似然估计 fit - FitML(spec, series ret) summary(fit)CreateSpec的参数含义variance.spec list(model sGARCH)指定每个状态内部用的是标准 GARCH(1,1) 方差方程distribution.spec控制扰动项分布可选norm、std学生 t、sstd偏 t等switch.spec list(K 2)就是马尔可夫转换的状态个数。对于日频收益率我一般直接用std而不是norm因为收益尾部更厚但厚尾分布会显著增加优化时间如果只是先跑通流程用norm更合适。FitML返回的 summary 里重点看三类内容。第一是每个状态的参数估计包括 ω_1、ω_2、α_1、α_2、β_1、β_2注意对比两个状态的 αβ 之和它反映波动率衰减速度高波动状态的 αβ 通常更接近 1说明冲击持续更久。第二是转移概率矩阵 p、q如果 p 估计为 0.99 以上说明状态 1 非常持久。第三是持续期期望1/(1-p) 和 1/(1-q) 直接用笔算就行。3.3 提取状态概率与绘制平滑概率图参数估计完下一步就是把每个时刻属于状态 1 和状态 2 的概率提取出来。MSGARCH包提供了StateProb和PosteriorProb两个函数前者返回滤波概率后者返回平滑概率。平滑概率用了全样本信息事后看更准实时决策只能用滤波概率。# 平滑概率每个时点位于各状态的后验概率 sp - PosteriorProb(fit) matplot(sp, type l, lty 1, col c(steelblue, firebrick), xlab t, ylab Probability, main Smoothed State Probabilities) legend(topright, legend c(State 1 (Low Vol), State 2 (High Vol)), col c(steelblue, firebrick), lty 1)画图后通常能明显看到高波动状态概率在 2015 年、2020 年、2022 年出现尖峰状抬升。如果画出来两个状态概率一直在 0.5 上下纠缠不清说明模型没有把两状态区分开后续使用要谨慎。这是模型失败的最常见信号之一原因可能出在初始值或者状态数选择上。处理办法放在第五章专门讲。4. Python 里的替代路线自写似然函数复现两状态 MS-GARCH4.1 为什么选择从零写似然而不是直接找现成包装很多人希望像 R 一样pip install ms_garch一把梭。现实是 Python 生态里没有一个统一维护、开箱即用的 MS-GARCH 包。有的库叫markovregime或者hmmlearn但前者只支持均值转换不支持方差方程里的 GARCH 结构后者是通用隐马尔可夫模型观测分布必须是同一族分布不能天然和 GARCH 递归耦合。因此实际工作中Python 落地 MS-GARCH 有两条路一是通过rpy2直接把 R 的MSGARCH包拉过来跑适合不想重复造轮子的场景二是自己写负对数似然函数配合scipy.optimize或者statsmodels的优化器做估计。本节给出第二种方案的最小骨架它还能帮你把第二章没想清楚的 Hamilton 滤波和 Gray 折叠落到实处。4.2 负对数似然与状态概率滤波的代码骨架两状态 MS-GARCH 的似然计算核心是 Hamilton 滤波。每个时点 t 迭代四步计算每个状态的观测密度、计算混合密度、更新滤波概率、预测下一期状态概率。方差递归使用 Gray 折叠即在每个时点用滤波概率加权的混合方差作为下一期方差递归的输入。这样处理虽然会损失一部分理论上的一致性但工程上稳定不容易发散是实证中最常见的做法。import numpy as np from scipy.stats import norm from scipy.optimize import minimize def ms_garch_nll(params, r): # 参数顺序: [mu1, mu2, omega1, omega2, alpha1, alpha2, beta1, beta2, p_logit, q_logit] mu params[0:2] omega np.exp(params[2:4]) # 保证方差为正 alpha np.exp(params[4:6]) # 保证非负 beta np.exp(params[6:8]) # 保证非负 p 1.0 / (1.0 np.exp(-params[8])) # 反logit到(0,1) q 1.0 / (1.0 np.exp(-params[9])) T len(r) xi np.array([0.5, 0.5]) # 初始状态概率取均匀 sigma2 np.array([np.var(r), np.var(r)]) # 两条状态方差路径 loglik 0.0 for t in range(T): # 1. 每个状态在当前值下的正态密度 eta norm.pdf(r[t], locmu, scalenp.sqrt(sigma2)) # 2. 混合密度用于累积似然 dens np.dot(xi, eta) if dens 1e-12: dens 1e-12 loglik np.log(dens) # 3. 滤波后验概率 xi_f xi * eta / dens # 4. 预测下一期状态概率 xi np.array([p * xi_f[0] (1 - q) * xi_f[1], (1 - p) * xi_f[0] q * xi_f[1]]) # 5. Gray折叠更新方差 mixed_var np.dot(xi_f, sigma2) sigma2 omega alpha * (r[t] ** 2) beta * mixed_var return -loglik这段代码有几个必须说明的参数细节。第一ω、α、β 用exp参数化是为了让优化器在无约束空间里工作同时保证估计值非负如果你直接用原始参数scipy.optimize很容易在迭代过程中把方差推成负数报错或者收敛到 NaN。第二p 和 q 用反 logit 变换映射到 (0,1)避免转移概率越界。第三初始状态概率取 {0.5, 0.5} 并且不加约束对结果影响不大但如果样本里两个状态占比极不均衡用 {0.9, 0.1} 之类的先验值有时能帮助优化器更快找到好区域。第四方差初值取全样本方差比取 0 或者取第一个观测本身的残差平方更平稳能显著减少迭代初期的数值振荡。4.3 用 scipy.optimize 做数值优化的参数陷阱有了目标函数后调用优化器的写法如下# 初始化均值取样本均值附近方差参数取小正数p/q 取 0.9 init [0.0, 0.0, np.log(np.var(ret)/2), np.log(np.var(ret)*2), np.log(0.05), np.log(0.1), np.log(0.9), np.log(0.85), 2.0, 2.0] # logit(0.9)≈2.197 res minimize(ms_garch_nll, init, args(ret,), methodL-BFGS-B, options{maxiter: 2000}) print(res.x, res.fun)这里最容易翻车的点是scipy.optimize默认的迭代次数远不够 MS-GARCH 这种高维非凸问题收敛建议至少给maxiter2000。即便如此一次优化也经常落进局部最优。我常用的策略是并行做 20 到 50 组随机初值取负对数似然最小的那一组作为最终解。具体操作时用numpy.random.uniform给均值、方差参数和转移概率各生成一组初值然后把每组初值代入优化器最后对比res.fun。这个多起点搜索在 Python 里跑起来很快因为每个序列长度也就几百到几千单次优化耗时通常不到一秒到几秒。另一个容易忽略的问题是收敛标志。res.success为True只代表优化器行程结束不代表找到全局最优。更加可本文还有配套的精品资源点击获取