马尔可夫切换时变参数自回归模型的贝叶斯推断

Bayesian Inference for Markov-Switching Autoregressive Model with Time-Varying Parameter

摘要: 本文针对马尔可夫切换时变参数自回归模型构建一套完整的贝叶斯推断框架。传统极大似然估计在处理高维时变参数与潜状态转移耦合的问题时，常存在数值不稳定、易陷于局部最优及不确定性量化不足等局限。为应对这些挑战，本研究设计一种高效的马尔可夫链蒙特卡洛抽样算法，整合对数尺度下的隐状态前向滤波–后向抽样、时变系数的卡尔曼平滑，以及基于共轭先验的转移概率吉布斯抽样步骤，并通过引入平稳性先验约束确保动态过程的合理性。模拟实验表明，该算法收敛稳定、估计精度高：主要参数的后验均值相对偏差普遍低于5%，潜在状态路径恢复准确率达100%，样本内拟合效果良好。本研究为同时蕴含结构突变与参数渐进演化特征的时间序列分析提供一个计算可靠、推断完备的贝叶斯解决方案。

Abstract: This paper develops a comprehensive Bayesian inference framework for Markov-switching time-varying parameter autoregressive model. Traditional maximum likelihood estimation often encounters limitations such as numerical instability, a tendency to converge to local optima, and insufficient uncertainty quantification when dealing with high-dimensional time-varying parameters intertwined with latent regime shifts. To address these challenges, this study designs an efficient Markov chain Monte Carlo sampling algorithm that integrates log-scale forward filtering-backward sampling for latent states, Kalman smoothing for time-varying coefficients, and Gibbs sampling steps for transition probabilities based on conjugate priors. Stationarity-inducing prior constraints are also incorporated to ensure the plausibility of the dynamic process. Simulation experiments demonstrate that the proposed algorithm exhibits stable convergence and high estimation accuracy: the relative deviations of posterior means for key parameters are generally below 5%, the recovery accuracy of latent regime paths reaches 100%, and in-sample fitting performance is satisfactory. This research provides a computationally reliable and inferentially complete Bayesian solution for analyzing time series that exhibit both structural breaks and gradual parameter evolution.

1. 引言

时间序列数据中结构突变与参数渐变并存的现象，在宏观经济波动、金融市场状态切换以及环境系统演变等多个领域[1]-[3]广泛存在。马尔可夫切换时变参数自回归(MS-TVP-AR)模型[4]通过将状态切换的动态离散性与状态内参数的连续时变性相融合，为此提供了具有吸引力的理论框架，但其复杂的双重潜变量结构也为统计推断带来了实质性挑战。传统的极大似然估计方法依赖于非线性滤波技术[5]与数值优化算法。在应对MS-TVP-AR模型时，该方法存在三方面主要局限：高维参数空间与高度非凸的似然函数面使得优化过程极易陷入局部最优，估计结果对初值选择异常敏感；基于频率学框架的点估计难以对时变参数路径与状态识别结果提供完整的不确定性量化；该框架缺乏纳入先验理论信息的自然机制，限制了其在有限样本下的应用稳健性。

贝叶斯推断与马尔可夫链蒙特卡洛方法[6]的结合，为处理此类复杂潜变量模型提供了系统的解决方案。通过将状态序列、时变参数及超参数全部纳入随机变量框架并施加合理的先验分布，不仅能自然融合先验信息，更能通过后验抽样直接获得所有参数的完整分布特征，从而实现对估计不确定性的严谨度量。贝叶斯方法在时变自回归(TVP-AR)模型的估计与推断中已得到长足发展[7]-[9]。进一步地，马尔可夫切换时变参数自回归(MS-TVP-AR)模型，可被构建为一种包含双重潜变量的状态空间形式。尽管贝叶斯框架在状态空间模型[10]-[13]及其子类的估计中已较为成熟，但目前尚缺乏文献将其系统、完整地应用于兼具“状态切换”与“参数时变”双重属性的MS-TVP-AR模型。因此，本文填补这一方法应用上的空白。本研究旨在为MS(2)-TVP-AR(1)模型构建一套完整、稳健的贝叶斯MCMC推断框架。本文的方法论贡献主要体现在三个方面的技术整合：一是通过对数尺度前向滤波–后向抽样(FFBS)算法，解决序列分析中的数值下溢问题，显著提升状态变量抽样的效率与数值稳定性；二是将平稳性先验约束系统性地嵌入时变参数的贝叶斯更新过程中，通过截断正态先验与拒绝抽样机制，保障了模型在经济意义上的长期合理性；三是设计并实现完整的模块化Gibbs采样框架，保证抽样精度。通过模拟研究，我们验证了该框架在参数估计精度、状态识别准确率及计算稳定性方面的优良性能。

2. 模型

考虑一个可观测的单变量时间序列${\left\{y_t\right\}}_{t=1}^T$ ，本文主要考虑阶数$p=1$ ，状态切换数$K=2$ 的马尔可夫切换时变参数自回归模型[4]。观测方程建立了$y_t$ 与其滞后项、状态依赖的时变参数及误差项之间的联系：

$y_t=c_{s_t}+{\rho }_{t,s_t}{y}_{t-1}+{\epsilon }_t,{\epsilon }_t\sim \mathcal{N}\left(0,{\sigma }_{s_t}^2\right),t=1,\cdots ,T,$

其中$s_t\in \left\{1,2\right\}$ 表示时刻$t$ 所处的状态，$c_{s_t}$ 为状态依赖的截距项，${\rho }_{t,s_t}$ 为状态依赖的时变自回归系数，${\sigma }_{s_t}^2$ 为状态依赖的观测噪声方差。时变系数${\rho }_{t,s_t}$ 的演化由以下一阶自回归方程描述：

${\rho }_{t,s_t}=d_{s_t}+{\varphi }_{s_t}{\rho }_{t-1,s_{t-1}}+{\eta }_t,{\eta }_t~\mathcal{N}\left(0,{\tau }_{s_t}^2\right),t=1,\cdots ,T,$

其中$d_{s_t}$ 为漂移项，${\varphi }_{s_t}$ 为满足平稳性约束$\left|{\varphi }_{s_t}\right|<1$ 的持久性参数，${\tau }_{s_t}^2$ 为状态依赖的过程扰动方差。不可观测的状态序列$\left\{s_t\right\}$ 服从一阶齐次马尔可夫链，其转移概率矩阵为：

$P=\left[\begin{array}{cc}{p}_{11}& 1-p_{11}\\ 1-p_{22}& p_{22}\end{array}\right],p_{ij}=\text{Pr}\left(s_t=k|s_{t-1}=j\right),$

初始状态$s_0$ 设定为该链的平稳分布$\pi ={\left({\pi }_1,{\pi }_2\right)}^{\top }$ ，其中${\pi }_1=\left(1-p_{22}\right)/\left(2-p_{11}-p_{22}\right)$ ， ${\pi }_2=\left(1-p_{11}\right)/\left(2-p_{11}-p_{22}\right)$ 。该模型同时捕捉两种动态切换：离散的马尔可夫链$\left\{s_t\right\}$ 反映数据的结构性突变，而连续的时变系数$\left\{{\rho }_{t,s_t}\right\}$ 刻画同一状态内回归关系的渐进演化，从而为分析兼具突变与渐变特征的复杂时间序列提供了一种贴合实际的建模框架。

注1：本文中，我们使用$k,j\in \left\{1,2\right\}$ 作为状态索引，用以标识所有依赖于状态的模型参数集合，例如$c_k,{\sigma }_k^2$ 。

3. 贝叶斯推断

本节构建MS(2)-TVP-AR(1)模型的贝叶斯推断框架和MCMC抽样。模型包含三类未知量：固定参数集$\Theta ={\left\{c_k,{\sigma }_k^2,d_k,{\varphi }_k,{\tau }_k^2,p_{11},p_{22}\right\}}_{k=1}^2$ 、潜在状态序列$s=\left(s_0,s_1,\cdots ,s_T\right)$ 以及时变参数路径$\rho =\left\{{\rho }_{t,k}\right\}$ 。状态变量$s_t\in \left\{1,2\right\}$ 服从齐次马尔可夫链，其平稳分布为${\pi }_1=\left(1-p_{22}\right)/\left(2-p_{11}-p_{22}\right)$ ，${\pi }_2=1-{\pi }_1$ 。对参数设定如下独立共轭先验：

$\begin{array}{ll}{c}_k~\mathcal{N}\left({\mu }_k^c,{\omega }_k^2\right),\hfill & {\sigma }_k^2~\text{IG}\left({\alpha }_k^{\sigma },{\beta }_k^{\sigma }\right),\hfill \\ d_k~\mathcal{N}\left({\mu }_k^d,{\nu }_k^2\right),\hfill & {\varphi }_k~\mathcal{N}\left({\mu }_k^{\varphi },{\kappa }_k^2\right)\mathbb{I}\left(\left|{\varphi }_k\right|<1\right),\hfill \\ {\tau }_k^2~\text{IG}\left({\alpha }_k^{\tau },{\beta }_k^{\tau }\right),\hfill & p_{kk}~\text{Beta}\left({\alpha }_{kk},{\beta }_{kk}\right),\hfill \\ {\rho }_{0,k}~\mathcal{N}\left({\mu }_k^{{\rho }_0},{\zeta }_k^2\right),\hfill & k=1,2.\hfill \end{array}$

其中正态分布与逆伽玛分布的选取基于共轭性考虑，平稳性约束$\mathbb{I}\left(\left|{\varphi }_k\right|<1\right)$ 保证了时变参数演化过程的经济合理性。给定观测数据 $y$ ，联合后验分布正比于：

$p\left(\Theta ,s,\rho |y\right)\propto p\left(y|s,\rho ,\Theta \right)p\left(\rho |\Theta \right)p\left(s|\Theta \right)p\left(\Theta \right),$

其中：

$\begin{array}{l}p\left(y|\cdot \right)={\displaystyle \prod _{t=1}^T\mathcal{N}}\left(y_t|c_{s_t}+{\rho }_{t,s_t}{y}_{t-1},{\sigma }_{s_t}^2\right),\\ p\left(\rho |\Theta \right)={\displaystyle \prod _{k=1}^2\left[\mathcal{N}\left({\rho }_{0,k}|{\mu }_k^{{\rho }_0},{\zeta }_k^2\right){\displaystyle \prod _{t=1}^T\mathcal{N}}\left({\rho }_{t,k}|d_k+{\varphi }_k{\rho }_{t-1,k},{\tau }_k^2\right)\right]},\\ p\left(s|\Theta \right)={\pi }_{s_0}{\displaystyle \prod _{t=1}^T{p}_{s_{t-1},s_t}}.\end{array}$

为从联合后验分布$p\left(\Theta ,s,\rho \mid y\right)$ 中高效抽取样本，本文设计了一套基于吉布斯采样的MCMC算法。该算法采用模块化迭代框架，依次对状态序列$s$ 、时变参数路径$\rho$ 及固定参数$\Theta$ 进行条件抽样。这一综合设计有效应对模型因同时包含离散潜状态与连续时变参数而带来的计算复杂性。

注2：在算法描述中，符号$k$ 和$j$ 将继续用于遍历状态标签($\left\{1,2\right\}$ )。在参数更新步骤中，$k$ 指代目标参数所在的状态；在状态空间模型的概率计算中，$k$ 则表示对“当前状态取值为$k$ ”这一假设的遍历。

状态序列的更新通过对数尺度的前向滤波–后向抽样(FFBS)算法[14]完成，以实现从条件后验分布$p\left(s|y,\rho ,\Theta \right)$ 中的高效稳健抽样。该算法的核心在于将对数滤波概率${\ell }_t\left(k\right)=\ln \text{Pr}\left(s_t=k|y_{1:t},\Theta \right)$ 的递推计算置于对数空间，从而避免长序列中因概率连乘导致的数值下溢。此对数尺度FFBS算法在数学上等价于标准算法，但通过数值稳定的递推公式，为长序列下的状态识别提供了可靠保障，是本文推断框架的基石。前向滤波过程始于平稳分布${\ell }_0\left(k\right)=\ln {\pi }_k$ ，并按以下步骤递推($t=1,\cdots ,T$ )：

(1) 预测步：计算${\ell }_{t|t-1}\left(k\right)=\ln \left[{\displaystyle {\sum }_{j=1}^2\exp \left({\ell }_{t-1}\left(j\right)+\ln p_{jk}\right)}\right]$ ，其中对数求和通过Log-Sum-Exp (LSE)技巧保持稳定。

(2) 更新步：纳入观测似然，${\ell }_t\left(k\right)={\ell }_{t|t-1}\left(k\right)+\ln \mathcal{N}\left(y_t;c_k+{\rho }_{t,k}{y}_{t-1},{\sigma }_k^2\right)$ ，其中$\ln \mathcal{N}(\cdot )$ 为正态分布的对数密度。

(3) 归一化：对${\ell }_t\left(k\right)$ 进行中心化与归一化处理，确保数值表示稳定。

后向抽样阶段首先依据$\text{exp}\left({\ell }_T\left(k\right)\right)$ 抽取$s_T$ ，随后逆向递推($t=T-1,\cdots ,1$ )：在给定$s_{t+1}=j$ 的条件下，依概率$\text{Pr}\left(s_t=k|\cdot \right)\propto \text{exp}\left({\ell }_t\left(k\right)\right)\cdot p_{kj}$ 抽取$s_t$ 。

对于时变参数路径的更新，在给定状态序列$s$ 后，转化为对两个独立线性高斯状态空间模型的抽样，利用状态空间模型的结构实现整条路径的高效采样。对于每个状态$k$ ，我们采用标准的卡尔曼滤波与蒋森平滑(RTS平滑器)算法[15] [16]，从条件后验分布$p\left({\rho }^{\left(k\right)}|y,s,\Theta \right)$ 中联合抽取整条路径 ${\rho }^{\left(k\right)}={\left({\rho }_{0,k},\cdots ,{\rho }_{T,k}\right)}^{\top }$ 。该算法通过前向滤波获得状态的条件分布，再通过后向平滑递归进行抽样，其平滑步的核心公式为：

${\rho }_{t,k}|{\rho }_{t+1,k},\cdots ~\mathcal{N}\left(a_t+J_t\left({\rho }_{t+1,k}-a_{t+1|t}\right),P_t-J_t^2{P}_{t+1|t}\right),$

其中$J_t={\varphi }_k{P}_t/P_{t+1|t}$ ，$a_t$ 和$P_t$ 为前向滤波的输出。

在贝叶斯推断框架下，固定参数$\Theta$ 的条件后验分布得益于共轭先验设计，具有解析形式，从而能够实现高效的精确抽样，无需依赖近似算法。观测方程参数$c_k$ 和${\sigma }_k^2$ 的条件后验分布可直接推导。对于每个状态$k$ ，在给定状态序列$s$ 和时变参数${\rho }_k$ 后，定义属于该状态的观测集合$T_k=\left\{t:s_t=k\right\}$ 及调整后的观测值${\tilde{y}}_t=y_t-{\rho }_{t,k}{y}_{t-1}$ ($t\in T_k$ )。观测误差方差${\sigma }_k^2$ 的条件后验为逆伽马分布：

${\sigma }_k^2|\cdots ~\text{Inv-Gamma}\left({\alpha }_k^{\sigma }+\frac{n_k}{2},{\beta }_k^{\sigma }+\frac{1}{2}{\displaystyle \sum }_{t\in T_k}{\left({\tilde{y}}_t-c_k\right)}^2\right),$

其中$n_k=\left|T_k\right|$ 。给定${\sigma }_k^2$ ，截距项$c_k$ 的条件后验为正态分布：

$c_k|{\sigma }_k^2,\cdots ~\mathcal{N}\left({\mu }_k^{c\text{*}},V_k^{c\text{*}}\right),$

其后验均值和方差分别为${\mu }_k^{c\text{*}}=V_k^{c\text{*}}\left({\mu }_k^c/{\omega }_k^2+{\displaystyle {\sum }_{t\in T_k}{\tilde{y}}_t/{\sigma }_k^2}\right)$及$V_k^{c\text{*}}={\left(1/{\omega }_k^2+n_k/{\sigma }_k^2\right)}^{-1}$ 。状态方程参数$d_k$ 、${\varphi }_k$ 和${\tau }_k^2$ 的条件后验可通过将系数演化方程重写为线性回归形式得到。定义向量$Y_k={\left({\rho }_{1k},\cdots ,{\rho }_{T,k}\right)}^{\top }$ 和设计矩阵$X_k=\left[1,{\left({\rho }_{0,k},\cdots ,{\rho }_{T-1,k}\right)}^{\top }\right]$ 。过程噪声方差${\tau }_k^2$ 的条件后验为：

${\tau }_k^2|\cdots ~\text{Inv-Gamma}\left({\alpha }_k^{\tau }+\frac{T}{2},{\beta }_k^{\tau }+\frac{1}{2}{\left(Y_k-X_k{\beta }_k\right)}^{\top }\left(Y_k-X_k{\beta }_k\right)\right),$

其中${\beta }_k={\left(d_k,{\varphi }_k\right)}^{\top }$ 。给定${\tau }_k^2$ ，回归系数${\beta }_k$ 的条件后验为正态分布：

${\beta }_k|{\tau }_k^2,\cdots ~\mathcal{N}\left({\mu }_k^{\beta \text{*}},{\Sigma }_k^{\beta \text{*}}\right),$

其中${\Sigma }_k^{\beta \text{*}}={\left({\left({\Sigma }_k^{\beta }\right)}^{-1}+X_k^{\top }{X}_k/{\tau }_k^2\right)}^{-1}$ ，${\mu }_k^{\beta \text{*}}={\Sigma }_k^{\beta \text{*}}\left({\left({\Sigma }_k^{\beta }\right)}^{-1}{\mu }_k^{\beta }+X_k^{\top }{Y}_k/{\tau }_k^2\right)$ ，且${\mu }_k^{\beta }={\left({\mu }_k^d,{\mu }_k^{\varphi }\right)}^{\top }$ ，${\Sigma }_k^{\beta }=\text{diag}\left({\nu }_k^2,{\kappa }_k^2\right)$ 为先验矩。为满足平稳性要求，对${\varphi }_k$ 的抽样需进行拒绝抽样，仅接受满足$\left|{\varphi }_k\right|<1$ 的样本。转移概率参数$p_{11}$ 和$p_{22}$ 的条件后验基于状态序列的转移计数。令$n_{kj}=\#\left\{t:s_{t-1}=k,s_t=j\right\}$ 为从状态$k$ 到$j$ 的转移次数。结合独立的Beta先验，其后验分布为共轭的Beta分布，可直接精确抽样：

$p_{11}|s~\text{Beta}\left({\alpha }_{11}+n_{11},{\beta }_{11}+n_{12}\right),p_{22}|s~\text{Beta}\left({\alpha }_{22}+n_{22},{\beta }_{22}+n_{21}\right).$

平稳分布的更新是算法迭代中的必要步骤。每次获得新的转移概率$p_{11}^{\left(m\right)}$ 和$p_{22}^{\left(m\right)}$ 后，需重新计算对应马尔可夫链的平稳分布${\pi }^{\left(m\right)}={\left({\pi }_1^{\left(m\right)},{\pi }_2^{\left(m\right)}\right)}^{\top }$ ，其解析表达式为：

${\pi }_1^{\left(m\right)}=\frac{1-p_{22}^{\left(m\right)}}{2-p_{11}^{\left(m\right)}-p_{22}^{\left(m\right)}},{\pi }_2^{\left(m\right)}=\frac{1-p_{11}^{\left(m\right)}}{2-p_{11}^{\left(m\right)}-p_{22}^{\left(m\right)}}.$

更新后的平稳分布将用于下一轮MCMC迭代中状态序列$s$ 的前向滤波初始化步骤。

以上所有步骤构成了完整的Gibbs采样循环，各参数块的条件后验均具有标准形式，确保了算法的高效性与稳定性。MS(2)-TVP-AR(1)模型的完整MCMC抽样算法流程总结如下：

(1) 初始化：设定参数初值${\Theta }^{\left(0\right)}$ 、状态序列$s^{\left(0\right)}$ 及时变系数${\rho }^{\left(0\right)}$ 。设定先验超参数、总迭代次数$M$ 与预烧期$B$ 。

(2) 迭代抽样(对于$m=1,\cdots ,M$ )：

(a) 更新状态序列：使用对数尺度FFBS算法，基于$y,{\rho }^{\left(m-1\right)},{\Theta }^{\left(m-1\right)}$ 抽取$s^{\left(m\right)}$ 。

(b) 更新时变参数：对于$k=1,2$ ，使用卡尔曼平滑算法，基于$y,s^{\left(m\right)},{\Theta }^{\left(m-1\right)}$ 抽取${\rho }^{\left(k\right)\left(m\right)}$ 。

(c) 更新测量方程参数：对于$k=1,2$ ，从共轭后验中抽取${\sigma }_k^{2\left(m\right)}$ 与$c_k^{\left(m\right)}$ 。

(d) 更新转移方程参数：对于$k=1,2$ ，从共轭后验中抽取${\tau }_k^{2\left(m\right)}$ 、$d_k^{\left(m\right)}$ 与${\varphi }_k^{\left(m\right)}$ 。

(e) 更新转移概率和平稳分布：抽样$p_{11}^{\left(m\right)},p_{22}^{\left(m\right)}$ ；抽样${\pi }_1^{\left(m\right)},{\text{π}}_2^{\left(\text{m}\right)}$ 。

(3) 输出：丢弃前$B$ 次预烧样本，保留剩余的样本集合${\left\{{\Theta }^{\left(m\right)},s^{\left(m\right)},{\rho }^{\left(m\right)}\right\}}_{m=B+1}^M$ 用于后续的后验统计推断。

该算法通过融合对数尺度FFBS、卡尔曼平滑与Gibbs采样模块，为MS(2)-TVP-AR(1)模型提供了一个计算稳定、效率良好且易于实现的贝叶斯推断完整解决方案。

尽管本文聚焦于基本形式，但所构建的贝叶斯推断框架具备向更一般的MS(K)-TVP-AR(p)模型扩展的理论基础与算法模块。当状态数$K$ 与自回归阶数$p$ 的增加，时变参数向量的维度需要扩展，状态方程需以向量化的形式表示，观测方程需扩展为包含$p$ 阶滞后项。状态序列抽样仍可采用对数尺度前向滤波—后向抽样算法，但转移概率矩阵$P$ 的维度变为$K\times K$ 。时变参数的抽样仍可使用卡尔曼平滑算法，但需处理更高维的状态向量与设计矩阵。大多数固定参数的共轭更新步骤形式不变，需调整参数维度。扩展后模型参数数量显著增加，对MCMC抽样的混合效率与计算资源提出更高要求。同时，模型可能存在更严重的似然函数平坦区与状态标签互换等识别问题，需通过施加更细致的先验约束来解决。

4. 模拟实验

为全面评估所提MS(2)-TVP-AR(1)模型的贝叶斯MCMC算法的性能，本节设计了模拟实验。实验采用经济学领域常见的季度GDP增长率数据作为参考，生成符合MS(2)-TVP-AR(1)模型的模拟时间序列，并对所提方法进行多维度评估。代码实现基于R语言，核心算法模块包括对数尺度FFBS、卡尔曼滤波平滑及参数的条件共轭抽样。

参考GDP增长率的典型特征，模型设定两个状态：状态1表示扩张期，状态2表示收缩期。真实参数设定如表1所示，模拟数据时间序列和状态切换如图1所示。

Table 1. True parameters of simulated data

表1. 模拟数据真实参数

参数模块	参数	状态1	状态2	经济学解释
观测方程	$c_k$	1.200	−0.500	状态下的平均增长率
	${\sigma }_k^2$	0.010	0.015	观测波动性
状态方程	$d_k$	0.100	−0.050	时变系数的长期漂移
	${\varphi }_k$	0.850	0.750	时变系数演化持续性，满足$\left|{\varphi }_k\right|<1$
	${\tau }_k^2$	0.002	0.005	时变系数演化不确定性
转移机制	$p_{kk}$	0.950	0.900	状态持续性

Figure 1. Time series and state transition diagram

图1. 时间序列和状态切换图

基于上述参数，模拟25年季度数据，生成长度$T=100$ 的时间序列。生成的数据中，扩张期实际占比61.0%，出现了5个完整的扩张阶段，平均长度为12.2季度；收缩期占比39.0%，出现了5个阶段，平均长度为7.8季度。生成的GDP增长率序列范围在[−0.71%, 4.34%]之间，均值为1.38%，标准差为1.58%，与设定的参数特征相符。

采用两条独立的MCMC链进行后验推断，每条链运行8000次迭代，舍弃前3000次作为烧退期，之后每5次迭代保留一个样本，最终每条链获得1000个有效样本。算法为所有未知参数设置了如表2的共轭或弱信息先验分布。

Table 2. Prior distribution of parameters

表2. 参数的先验分布

参数	先验分布	超参数设定
$c_k$	$\mathcal{N}\left({\mu }_k^c,{\omega }_k^2\right)$	${\mu }^c=\left(1.25,-0.55\right),\omega =\left(1.0,1.0\right)$
${\sigma }_k^2$	$\text{Inv-Gamma}\left({\alpha }_k^{\sigma },{\beta }_k^{\sigma }\right)$	${\alpha }^{\sigma }=\left(2.5,2.5\right),{\beta }^{\sigma }=\left(0.015,0.011\right)$
$d_k$	$\mathcal{N}\left({\mu }_k^d,{\nu }_k^2\right)$	${\mu }^d=\left(0.05,-0.05\right),\nu =\left(0.3,0.3\right)$
${\varphi }_k$	$\mathcal{N}\left({\mu }_i^{\varphi },{\kappa }_k^2\right)\mathbb{I}\left(\Vert {\varphi }_k\Vert <1\right)$	${\mu }^{\varphi }=\left(0.8,0.7\right),\kappa =\left(0.2,0.2\right)$
${\tau }_k^2$	$\text{Inv-Gamma}\left({\alpha }_k^{\tau },{\beta }_k^{\tau }\right)$	${\alpha }^{\tau }=\left(11,15\right),{\beta }^{\tau }=\left(0.02,0.005\right)$
$p_{11}$	$\text{Beta}\left({\alpha }_{11},{\beta }_{11}\right)$	${\alpha }_{11}=19,{\beta }_{11}=1$
$p_{22}$	$\text{Beta}\left({\alpha }_{22},{\beta }_{22}\right)$	${\alpha }_{22}=9,{\beta }_{22}=1$

采用Gelman-Rubin统计量$\widehat{R}$ 对两条MCMC链的收敛性进行综合评估。所有参数的$\widehat{R}$ 值均在1.000至1.006之间，远低于1.1的临界标准，且其97.5%分位数上限也均低于1.1，表明两条MCMC链已充分混合，收敛至同一后验分布。基于合并两条链的后验样本，计算各参数的后验均值、标准差及95%置信区间，并与真实值进行对比，结果如表3所示。各个参数的后验分布直方图如图2所示。

Table 3. Parameter estimation

表3. 参数估计

参数	真实值	后验均值	后验标准差	95%置信区间	区间包含真实值
$c_1$	1.200	1.233	0.037	[1.159, 1.303]	是
$c_2$	−0.500	−0.541	0.027	[−0.592, −0.489]	是
${\sigma }_1^2$	0.010	0.010	0.003	[0.006, 0.018]	是
${\sigma }_2^2$	0.015	0.011	0.003	[0.006, 0.019]	是
$d_1$	0.100	0.129	0.019	[0.092, 0.167]	是
$d_2$	−0.050	−0.057	0.019	[−0.092, −0.018]	是
${\varphi }_1$	0.850	0.810	0.035	[0.741, 0.876]	是
${\varphi }_2$	0.750	0.772	0.055	[0.646, 0.870]	是
${\tau }_1^2$	0.002	0.002	0.001	[0.001, 0.003]	是
${\tau }_2^2$	0.005	0.004	0.002	[0.001, 0.009]	是
$p_{11}$	0.950	0.937	0.027	[0.875, 0.978]	是
$p_{22}$	0.900	0.878	0.047	[0.769, 0.953]	是

表3的参数估计结果显示，在12个主要参数中，只有参数$c_2$ 的偏差略高于5%，绝大多数的后验均值相对偏差低于5%，表明模型整体上能够准确恢复参数。例如，核心参数$p_{11}$ 估计为0.937，真实值为0.950，偏差仅1.4%；${\varphi }_1$ 估计为0.810，真实值为0.850，偏差为4.7%。图2展示了各参数的MCMC采样后验分布，红色实线表示真实值，黑色虚线表示后验均值，灰色虚线分别表示95%置信区间的上限和下限。所有真实值均落在其对应的95%后验置信区间内，覆盖率为100%，表明贝叶斯方法能够合理地量化参数估计的不确定性。时变参数的自回归系数${\varphi }_1$ 和${\varphi }_2$ 的后验估计均满足$\left|{\varphi }_k\right|<1$ 的平稳性约束。该模型对模拟数据的样本内拟合如图3所示，由图可以看出拟合效果良好，RMSE为0.079%，MAE为0.062%，MAPE为8.394%。通过计算每个时间点上状态的后验概率，并取概率众数作为状态估计，模型对隐状态序列的恢复准确率达到100%，能准确识别经济周期切换时点。后验状态概率在真实状态发生转换的时刻呈现清晰、急剧的跳跃，表明模型能精准捕捉经济周期的转折点。

Figure 2. Histogram of posterior distribution

图2. 后验分布直方图

Figure 3. Growth rate fitting plot

图3. 增长率拟合图

5. 总结

本文提出了MS(2)-TVP-AR(1)模型的贝叶斯推断方法。算法通过引入对数尺度FFBS解决状态序列抽样的数值下溢问题，采用平稳性先验约束保障时变参数的合理演化，并利用共轭先验实现了转移概率和其他参数的高效更新。基于经济周期特征的模拟实验表明，该方法收敛稳健，所有参数的Gelman-Rubin统计量$\widehat{R}<1.01$ 。参数估计表现良好，核心参数的估计偏差普遍较小，95%后验置信区间均覆盖真实值，且对隐含状态序列的识别准确率达到100%。本文实现该模型首个完整的贝叶斯计算方案，其系统架构与验证结果为该模型在实证研究中的推广应用提供可靠的技术基础。