简介：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_t=j|S_{t-1}=i,S_{t-2}=\dots)=P(S_t=j|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'或'GED'；advOpt.stdEst取 1 表示估计后计算标准误，取 0 可以加速收敛，因为省去了海森矩阵数值计算。dep必须是T x 1列向量，indep必须是T x n矩阵，且第一列如果需要截距就直接放全 1 向量。MS_Regress_Fit返回两个结构体，常见字段如下表，具体字段名以who Spec_Out输出为准：
| 输出变量 | 字段 | 含义 |
|---|---|---|
Spec_Out | Coeff_1 | 状态 1 下的回归系数 |
Spec_Out | Coeff_2 | 状态 2 下的回归系数 |
Spec_Out | Sigma_1/Sigma_2 | 各状态的波动率 |
Spec_Out | P | k x k转移概率矩阵 |
Spec_Out | SmoothedProb | 平滑概率，T x k |
Optim_Out | LL | 对数似然值 |
Optim_Out | param | 最终参数向量 |
Optim_Out | hessian | 数值海森矩阵，用于标准误估计 |
SmoothedProb是最常被拿来画图的字段，它表示每一个时点处于每个状态的概率，在判断“哪一段时间属于状态 1、哪一段属于状态 2”时非常直观。Optim_Out.hessian不是每一次都会被计算，只有当advOpt.stdEst=1时才有，否则里面是空矩阵。如果运行后出现“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，模型似然都变化不大，那么点估计本身只是被优化器随便选了一个点，不具备实际解释意义。扫描时要注意：每次固定一个参数后重新估计，其他参数都在新的约束下自由调整，所以曲线本身的形状会比把整个模型在网格上重新估计更平滑。这个方法对状态转移概率的非对角元素尤其有效，因为转移概率对似然的影响往往是非线性的。扫描之后如果发现转移概率矩阵中的某个自由参数似然面太平坦，就应考虑把该元素固定为合理值，避免它在优化过程中产生数值噪声。
本文还有配套的精品资源，点击获取
