简介MS_Regress_FEX是一套面向Matlab用户的马尔科夫状态转换Regime Switching模型估计工具箱适合经济学、金融学等领域研究人员与数据分析师处理具有状态切换特征的时序数据。压缩包共38个文件以32个m函数文件为主辅以cpp/h源码接口、txt说明与PDF文档整体约337KB集成了模型设定、参数估计、模拟、预测和诊断等完整功能。目前已有164人浏览学习。使用者可通过附带的多元/多状态模拟与估计、t分布/GED分布、约束系数设定、MSVAR等多种示例脚本快速上手还可借助MEX滤波接口提升计算效率从而深入理解马尔科夫转换模型在宏观经济、金融市场和政策评估中的实际应用是一份兼具理论参考与代码实践的紧凑型工具箱。1. 为什么时间序列需要Markov Regime Switching从GNP数据到MS_Regress_FEX做宏观经济和金融时间序列最难受的一点是同一个回归系数对区间A和数据B都不适用。我早年处理美国实际GNP环比数据时衰退期和扩张期的截距、方差差得非常大直接用单一AR模型会让残差在1973和1982年附近出现明显的“系统性异常”。马尔科夫区制转换模型Markov Regime Switching也叫状态切换模型把这种参数变化直接建模成一个受马尔科夫链控制的状态变量也就是说系统在扩张、衰退、高波动、低波动这些状态之间转移状态的持续性和转移概率一起被估计出来。MS_Regress_FEX 是 MATLAB 生态里相当完整的实现压缩包里不仅有核心的 MS_Regress_Fit、MS_Regress_Sim 和滤波、似然函数还带了 Hamilton 原始数据的复现脚本、多变量扩展、Student-t/GED 分布假设甚至提供了一个mex_MS_Filter.cpp用来加速滤波循环。对已经懂一点状态空间或隐马尔科夫模型、但不想自己写滤波和优化循环的人这个包能直接把“模型设定”到“参数估计”之间的路缩短一大半。2. 模型设定与工具箱的核心函数MS_Regress_Fit、spec2param 与参数向量转换2.1 状态变量、转移概率与观测方程的整体结构马尔科夫区制转换模型的核心是假设一个不可直接观测的状态变量 (S_t \in {1,2,\dots,k}) 按一阶马尔科夫链变化即 (P(S_tj|S_{t-1}i,S_{t-2}\dots)P(S_tj|S_{t-1}i)p_{ij})。观测方程通常写成y_t X_t * beta_{S_t} sigma_{S_t} * epsilon_t也就是说回归系数和噪声方差完全由当前所处的状态决定。MS_Regress_FEX在实现时没有把状态序列当作可观测数据而是在似然函数里对所有可能的状态路径做滤波递推。滤波的核心是用上一期的预测概率和当期观测值更新状态概率再通过 Kim 平滑算法得到整个样本上的状态概率。这类模型对初始值、转移概率矩阵约束和误差分布形态都非常敏感工具箱里的checkInputs.m会先检查你传入的数据维度和选项避免把错误的k或constCoeff带进优化。2.2 spec2param 与 param2spec 在“模型结构体”和“参数向量”之间做了什么我刚开始用这个包时最不习惯的一点是你不需要直接给 fmincon 写一个超长的参数向量而是先把模型结构写成一个结构体再由工具包内部的spec2param.m把它压成优化器可以处理的向量。反过来当优化器迭代时param2spec.m再把这个向量还原回具有业务含义的结构体。这样做的意义在于转移概率矩阵的每一行都必须满足“求和为 1 且每个元素在 (0,1) 内”的约束如果直接优化原始概率矩阵很容易在边界处卡住。spec2param会把转移概率矩阵的自由元素转成满足约束的变换参数constCoeff则通过 NaN 和具体数值来标记哪些参数自由估计、哪些被固定。我一般不会手动调用这两个函数但调试时我会在命令行输入help spec2param help param2spec查看输入输出签名然后手动把一个估计结果转成参数向量、再加小扰动转回结构体用来生成下一轮的初值。这比直接在 fmincon 上操作参数索引要安全得多。2.3 MS_Regress_Fit 的调用签名与输出结构运行一个最简单的两状态模型只需要准备好因变量、自变量含截距和状态数。假设你已经把压缩包解压到本地用如下方式启动% 把工具箱代码目录加入 MATLAB 搜索路径 addpath(m_Files); addpath(data_Files); % 读取工具箱自带的示例数据 load(data_Files/Example_FEX.txt); data Example_FEX; % 第一列作为因变量只有截距的自变量矩阵 dep data(:,1); indep ones(size(dep,1),1); % 两个状态对应“低均值低波动”和“高均值高波动” k 2; % 选项结构体 advOpt.distrib normal; advOpt.stdEst 1; % 主估计函数 [Spec_Out, Optim_Out] MS_Regress_Fit(dep, indep, k, advOpt);代码里的advOpt.distrib控制误差分布normal是高斯分布后续可以换成t或GEDadvOpt.stdEst取 1 表示估计后计算标准误取 0 可以加速收敛因为省去了海森矩阵数值计算。dep必须是T x 1列向量indep必须是T x n矩阵且第一列如果需要截距就直接放全 1 向量。MS_Regress_Fit返回两个结构体常见字段如下表具体字段名以who Spec_Out输出为准输出变量字段含义Spec_OutCoeff_1状态 1 下的回归系数Spec_OutCoeff_2状态 2 下的回归系数Spec_OutSigma_1/Sigma_2各状态的波动率Spec_OutPk x k转移概率矩阵Spec_OutSmoothedProb平滑概率T x kOptim_OutLL对数似然值Optim_Outparam最终参数向量Optim_Outhessian数值海森矩阵用于标准误估计SmoothedProb是最常被拿来画图的字段它表示每一个时点处于每个状态的概率在判断“哪一段时间属于状态 1、哪一段属于状态 2”时非常直观。Optim_Out.hessian不是每一次都会被计算只有当advOpt.stdEst1时才有否则里面是空矩阵。如果运行后出现“Undefined function or variable”第一步检查addpath是否覆盖了m_Files目录第二步检查当前目录是否还停留在压缩包解压的根目录而不是m_Files里。3. 用示例脚本跑通一次两状态估计GNP_Hamilton.txt、输出与分布假设3.1 从 Hamilton 数据开始读取与预处理压缩包里的data_Files/GNP_Hamilton.txt是 Hamilton 1989 年那篇论文用过的实际 GNP 数据。很多复现脚本喜欢直接load这个文件但需要注意它的格式并不一定是单列可能包含日期或样本编号。我通常先看一眼raw load(data_Files/GNP_Hamilton.txt); size(raw) plot(raw)如果raw是多列先用raw(:,1)或raw(:,end)画出序列确认哪一列是你要建模的 GNP 增长率。Hamilton 原始数据中的序列通常已经做过差分或对数差分所以不需要额外做平稳性处理。读取完成后构建因变量和截距dep raw(:,1); % 选一列作为观测序列 indep ones(size(dep,1),1); % 模型只包含状态相关的均值只估计均值而没有自回归项是最容易跑通也最容易解释的设定。如果你发现模型的残差仍然有很明显的自相关那就说明需要在indep里加入滞后项比如indep [ones(T,1), raw(1:end-1,1)]并让样本对齐。3.2 运行 Example_MS_Regress_Fit.m输出与结果解读压缩包里的Example_MS_Regress_Fit.m是官方给的入门脚本。使用它最稳妥的方式是在 MATLAB 中打开文件逐段运行而不是直接在命令行输入文件名。我把它最核心的逻辑整理成如下可改写的模板% 清理当前工作区 clc; clear; addpath(m_Files); addpath(data_Files); % 读取数据 load(data_Files/Example_FEX.txt); data Example_FEX; % 因变量和自变量 dep data(:,1); indep [ones(size(dep,1),1), data(:,2:3)]; % 模型设定 k 2; advOpt.distrib normal; advOpt.stdEst 1; % 估计 [Spec_Out, Optim_Out] MS_Regress_Fit(dep, indep, k, advOpt); % 画平滑概率和区制划分 doPlots(Spec_Out, Optim_Out, dep, indep);doPlots.m是包内自带的绘图函数会生成状态概率和拟合值的图。估计完成后直接在命令行打印比较关键的数% 转移概率矩阵 disp(Spec_Out.P) % 平滑概率最后几行 disp(Spec_Out.SmoothedProb(end-5:end,:)) % 对数似然和AIC disp(Optim_Out.LL)转移概率矩阵的对角元素如果接近 0.9 或更高说明模型识别出了非常强的状态持续性。平滑概率的每一行之和约等于 1数值越接近 0 或 1 越好如果长期在 0.4 到 0.6 之间徘徊说明两个状态区分度不足通常需要加入更多解释变量或换成分位表现更好的分布。3.3 把误差分布换成 Student-t 与 GED自由度参数与适用场景宏观和金融数据往往有比正态分布更厚的尾部MS_Regress_FEX直接支持两种替代分布。修改方式非常简单% 使用 Student-t 分布 advOpt.distrib t; [Spec_Out_t, Optim_Out_t] MS_Regress_Fit(dep, indep, k, advOpt); % 使用 GED 广义误差分布 advOpt.distrib GED; [Spec_Out_GED, Optim_Out_GED] MS_Regress_Fit(dep, indep, k, advOpt);换成t之后工具箱会额外估计一个“自由度”参数自由度越低表示尾部越厚换成GED则多一个形状参数取值为 2 时退化为正态分布小于 2 时尾部比正态厚。选择哪种分布不能只看似然值谁大因为自由度参数多的模型本来就会获得更高的似然。我一般的做法是同时运行三种分布然后用 AIC/BIC 比较% 手算 AIC LL_normal Optim_Out.LL; nParam_normal length(Optim_Out.param); AIC_normal -2*LL_normal 2*nParam_normal; LL_t Optim_Out_t.LL; nParam_t length(Optim_Out_t.param); AIC_t -2*LL_t 2*nParam_t;三种分布适用场景的差异整理如下分布额外参数典型场景注意点normal无宏观变量、对数产出估计最快但对异常值敏感t自由度nu金融收益率、波动率需要额外估计自由度似然面更平坦GED形状参数具有中等厚尾的序列自由度参数与形状参数含义容易混淆如果换分布后转移概率矩阵发生剧烈变化说明原来的“两状态”设定可能不够稳定此时优先检查数据是否包含异常突变点而不是继续增加分布复杂度。4. 多变量与 MSVAR 扩展constCoeff 约束、模拟和 MEX 加速4.1 从单方程到多变量MS_VAR_Fit 与 MultiVar 示例单方程模型能捕捉一个序列的区制切换但很多金融问题需要同时建模多个相关序列比如利率、通胀和产出之间的联动关系。压缩包里的Example_MS_Regress_Fit_MultiVar.m演示的就是多因变量的情况。多变量版本的核心与单变量一致只是dep变成T x n_series的矩阵而不是向量% 多变量模型示例 load(data_Files/Example_FEX.txt); data Example_FEX; % 取前三列作为三个因变量 dep data(:,1:3); indep ones(size(dep,1),1); k 2; advOpt.distrib normal; [Spec_MV, Optim_MV] MS_Regress_Fit(dep, indep, k, advOpt);多变量版本的输出字段会按序列展开比如Coeff_1变成矩阵每一列对应一个因变量。需要特别注意的是多变量模型的状态变量仍是同一个也就是说所有因变量共享同一个区制过程。如果你想要每个变量有自己的状态那已经不是MS_Regress_FEX的默认设计需要去改似然函数结构。4.2 用 constCoeff 锁定参数固定系数与部分约束在实际建模时经常会遇到“某个解释变量只在某一状态下显著”或者“我们希望固定某个参数以检验理论假设”的情况。MS_Regress_FEX用advOpt.constCoef来标记自由系数和固定系数规则是矩阵大小等于size(indep,2) x k元素为NaN表示该系数自由估计数值表示固定为指定值。下面这个例子把状态 2 的第一个解释变量系数固定为 0含义是该变量对状态 2 没有影响constCoef repmat(NaN, size(indep,2), k); constCoef(2, 2) 0; % 第2个解释变量在状态2中固定为0 advOpt.constCoef constCoef; [Spec_C, Optim_C] MS_Regress_Fit(dep, indep, k, advOpt);注意固定的系数仍然会占用参数向量的一部分吗答案是不会。build_constCoeff.m和checkSize_constCoeff.m会先检查constCoef的定义固定值不会进入可估参数集合所以Optim_C.param的长度会明显变短。这既是优点也是坑如果你固定了一个原本识别度很差的系数可能会让整个参数向量变成不可识别从而出现收敛到全局最优但似然面依然是平的。固定系数时要配合下一章的似然面扫描来判断该参数是否真的可识别。4.3 mex_MS_Filter.cpp 的编译与加速从 MATLAB 循环到 C状态概率滤波是这类模型最耗时的部分。工具箱里的mex_MS_Filter.cpp是专门用来加速滤波循环的官方脚本Example_MS_Regress_Fit_with_MEX.m演示了编译后的用法。首次使用前需要先编译cd(m_Files); mex -setup C mex mex_MS_Filter.cpp cd(..);如果已经装了 MinGW64 或 Visual Studio 编译器这一步一般不会报错。编译成功后运行带 MEX 的示例脚本工具包会自动优先加载编译好的 mex 函数。要注意的是mex_MS_Filter.cpp里用的接口是固定的如果你修改了模型结构比如增加了更多滞后项需要确保 C 代码里的数组尺寸与实际传入数据一致否则 MATLAB 可能直接崩溃。我自己实践下来四五参数的小模型加速不明显但三状态、五解释变量、样本量上千时加速效果通常是数量级的。调试时如果怀疑 MEX 版本有问题可以把 mex 文件临时改名让工具包退回纯 MATLAB 实现这样能更快定位是算法问题还是 C 接口问题。4.4 模拟数据验证估计MS_Regress_Sim 与 Simul_and_Fit 示例MS_Regress_Sim.m的作用是先生成一条符合已知参数的状态切换序列然后再把它当作观测数据去估计。这个功能对验证模型识别性非常关键。官方示例里的Example_MS_Regress_Simul_and_Fit_2_States.m和Example_MS_Regress_Simul_and_Fit_3_States.m都做了这个闭环测试。模拟数据的大致流程如下simOpt.beta [0.5 -0.5; 1.0 -2.0]; % 两状态各自两个解释变量的系数 simOpt.sigma [0.5; 0.8]; % 状态1和状态2的波动率 simOpt.P [0.9 0.1; 0.2 0.8]; % 转移概率矩阵 simOpt.nobs 500; [SimData, TrueStates] MS_Regress_Sim(simOpt);SimData是模拟出的观测值TrueStates是模拟器内部使用的真实状态序列。接下来直接用SimData作为dep、用对应的截距矩阵进行估计然后把估计出来的转移概率和模拟设定的P比较。如果差异很大只要模型本身没有 bug基本都是因为模拟样本量太小或初始值卡在了局部最优。模拟验证是我在使用任何新参数设定前必做的一步尤其当某个状态的时间占比低于 10% 时模型几乎不可能稳定识别出该状态的系数。5. 验证区制划分时的三个实用技巧初值扰动、平滑概率与似然面扫描5.1 从 Optim_Out.param 构造扰动初值绕过局部最优MS_Regress_Fit用的底层优化器是 fmincon它对初值非常敏感。最直接的验证方式就是多试几组初值。先把上一轮估计出的参数向量作为基准再用param2spec把它转成结构体加小扰动后放回新模型% 以第3章估计结果为基准 baseParam Optim_Out.param; for r 1:5 % 转换成结构体便于修改带业务含义的字段 initSpec param2spec(baseParam, Spec_Out); initSpec.Coeff_1 initSpec.Coeff_1 randn(size(initSpec.Coeff_1))*0.05; initSpec.Sigma_1 initSpec.Sigma_1 * (1 0.05*randn(1)); % 转回参数向量作为下一次估计的初值 x0 spec2param(initSpec); [Spec_r, Optim_r] MS_Regress_Fit(dep, indep, k, advOpt); end不同版本对advOpt里初值入口的命名略有差异有的叫initParam有的只能在调用结构里增加参数最稳妥的办法是看checkInputs.m的开头部分里面定义了所有可接受字段。比较多次运行的对数似然如果大部分结果都稳定在同一个LL附近说明这个模型是可靠的如果同一个数据换了初值就得到完全不同的系数和概率那就要缩短样本、减少状态数或增加解释变量。5.2 用平滑概率判断区制切换是否真的“锐利”一个好的区制切换模型平滑概率在大多数时候都应该非常接近 0 或 1。画出状态 2 的平滑概率plot(Spec_Out.SmoothedProb(:,2), LineWidth, 1.2); ylim([0 1]); grid on;如果图中看到大量时间段的概率在 0.4 到 0.6 之间徘徊说明状态 1 和状态 2 在观测方程层面几乎不可区分。此时不要急着加状态数更有效的方法是先检查是否有一个状态基本没有持续样本。我经常用一段简单逻辑统计“每个状态被平滑概率超过 0.7 锁定的时间占比”prob2 Spec_Out.SmoothedProb(:,2); locked_2 mean(prob2 0.7); locked_1 mean(prob2 0.3); fprintf(状态1占比 %.2f%%状态2占比 %.2f%%\n, locked_1*100, locked_2*100);如果两个比例加起来远低于 1说明模型有很大一部分时间处于模糊状态。另一个常见情况是某个状态只出现非常短的瞬间比如只有 3 到 5 个连续时点被分到状态 2那就需要考虑这个状态到底是在捕捉真实的区制转换还是在拟合异常值。一个实战经验先把异常值剔除掉重新估计如果状态 2 簇明显变宽说明原模型中的状态 2 就是为异常值准备的。5.3 用 constCoeff 做一维似然面扫描判断参数可识别性把某个系数固定为一系列不同的值然后让其他参数自由估计观察对数似然随该系数变化的曲线是判断参数是否可识别的有效手段。利用第 4 章介绍的constCoef就可以实现coefGrid -1:0.05:1; LL_grid zeros(size(coefGrid)); for i 1:length(coefGrid) cc repmat(NaN, size(indep,2), k); cc(2,1) coefGrid(i); % 固定状态1第2个解释变量系数 advOpt.constCoef cc; [~, Optim_temp] MS_Regress_Fit(dep, indep, k, advOpt); LL_grid(i) Optim_temp.LL; end plot(coefGrid, LL_grid, o-); grid on; xlabel(固定系数值); ylabel(对数似然);如果这条曲线在某个取值处有一个明显尖峰说明数据对系数是敏感的估计值可信如果曲线接近一条水平线说明不管该系数取 0 还是 1模型似然都变化不大那么点估计本身只是被优化器随便选了一个点不具备实际解释意义。扫描时要注意每次固定一个参数后重新估计其他参数都在新的约束下自由调整所以曲线本身的形状会比把整个模型在网格上重新估计更平滑。这个方法对状态转移概率的非对角元素尤其有效因为转移概率对似然的影响往往是非线性的。扫描之后如果发现转移概率矩阵中的某个自由参数似然面太平坦就应考虑把该元素固定为合理值避免它在优化过程中产生数值噪声。本文还有配套的精品资源点击获取