Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
52,754 characters · 18 sections · 53 citation commands
Multi-regime Markov-switching models with time-varying transition probabilities: An application to U.S. Treasury yields
To capture nonlinear regime-dependent and cyclical dynamics in economics and finance, hamilton1989new introduced Hamilton’s Markov Switching model which assumes change between unobserved regimes or states---such as recession vs expansion; high volatility vs. low volatility, bull market vs. bear market---follows a Markov chain of first order with constant transition probabilities. This model has since then gained its popularity and has been widely applied in macroeconomics and finance, including detecting recessions; identifying different levels of volatility and classifying the stock market regime Cai1994A,chauvet1995econometric, gray1996modeling, garcia1996analysis, Bai2011Conditional,Guidolin2011Markov,Doornik2013A,Berentsen2022Modelling. However, constant transition probabilities could be too restrictive in empirical macro- and finance applications, as the probability of moving among states should depend on the economic situation, duration in different states, or other economic variables. diebold1994regime extended the hamilton1989new model and allowed for transition probabilities depending on other economic indicators, making the regime change responsive to economic behavior. filardo1994business and Filardo1998Choosing embeds time‑varying transition probabilities (TVTP) in the nonlinear Hamilton model, using observed “information variables” (economic or financial covariates) to drive the probability of moving between expansion and contraction in business cycles. filardo1994business also shows that the Hidden Markov Model with TVTP can track business cycles more closely than constant probability models, and macro variables can help to predict the turning point. In the TVTP - Markov Switching (MS) models proposed by diebold1994regime and filardo1994business, the driver of the time-‐varying transition matrix are exogenous variables. creal2013generalized introduced a framework to update the parameters of the transition probability matrix based on the predictive likelihood score. bazzi2017time adopted the framework of creal2013generalized's score-driven model and formulated conditions for the estimated time-varying probabilities. Their empirical implementation is carried out on US industrial production growth data. The models in creal2013generalized and bazzi2017time belong to the Generalized Autoregressive Score (GAS) family of models, while the driver of the time‑varying transition matrix is the score of the conditional log‑likelihood.
Based on how the transition probabilities and Markov switching process are modeled, other types of TVTP -MS models are proposed, such as time-varying and second-order Markov Models proposed by Neale2016Regime and the regime-switching models proposed by Li2020Asymptotic whereas the driver of the transition probability is a latent autoregressive factor via a threshold rule. As TVTP-MS preserves Hamilton’s regime‑switching structure while allowing regime persistence and switching risks to respond to economic information or the previous data-stream, they can therefore further improve turning points detection, forecasting and model interpretability. Thus, TVTP-MS models are widely utilized, e.g., to identify the recession/expansion in business cycles Diebold2020Regime; specify the volatility of crude oil futures pricesFong2002A; document regime‑dependent relationships between output, money, and prices ravn1995stylized; or perform the industrial electricity load forecasting Berk2018Probabilistic.
Among different TVTP-MS models, the bazzi2017time GAS models are intuitive and economically interpretable. The model parameters evolve over time in response to new information, and that information is measured by the score; thus the transition probabilities can adjust directly to the observations, and the updating is proportional to the information content in the recent data. This model has also demonstrated strong performance in Monte Carlo simulations and outperforms constant probability specifications in likelihood, information criteria, and forecasting. Although the implementation in bazzi2017time is carried out to reveal the three regimes of industrial production growth data with three means and two variances, their simulation is carried out by Monte Carlo study for the two-regime models with two different means and common variance. As in empirical macroeconomics, the common situation of the regimes is divided into three regimes, whereas the underlying states being recession, stability, and growth with separate means and variances, it is necessary to investigate the performance of estimation result of the time-varying transition probabilities models by simulations for the situation of more than two regimes with different means and their corresponding variances, and our paper will fill this gap. In our simulations, the transition probability will follow three types of dynamic changing as: (a)- transition probability changes with the lag value of the observations; (b) - transition probability changes with exogenous explanatory process; (c) - transition probabilities evolves based on the score of the predictive likelihood.Those three types of dynamics will correspond to model (I), (II) and (III) in the paper. We develop an open source R package multiregimeTVTP for the implementation of simulation in the general case of $K, K\geq 2$ regimes. Different performance metrics are used to evaluate estimation performance and mis-specification analysis is carried out to investigate the forecasting robustness to the choice of TVTP specification. The empirical implementation is carried out on U.S.\ Treasury zero-coupon yield data sorted by liu2021reconstructing.
The paper is divided into follow sections: section (ref) is the introduction of general structure of the regime-switching model with time-varying transition probabilities (TVTP) and the description of three types of TVTP models whereas the transition probabilities follow three different dynamics. Section (ref) presents the estimation result based on comprehensive simulation studies under the assumption of $K, K\geq2$ regimes whereas each regime has its own mean and variance. Section (ref) is the empirical analysis and the economic interpretation of the result. Section (ref) is conclusion and suggestions for future research.
Let $Y = (y_1,y_2, \dots, y_T )$ denote a time series of $T$ univariate observations of a stochastic process ${y_t}$, whereas its distribution depends on the realizations of the hidden discrete Markov process ${z_t}$ with finite state space ${1,2,...K}$ and the transition probability between different spaces is $\pi_{ij} = P(z_{t}=j|z_{t-1}=i); i, j = 1,2,...,K;\sum_{j=1}^{K} \pi_{ij}=1$. We further assume that the random variables ${y_t}$ are conditionally independent given ${z_t}$ with density $ (y_{t}|z_{t}=i)\sim \mathcal{N}\left(\mu_{i},\sigma_{i}\right)$. Let $\psi$ be the non-regime-specific parameters and $\theta_i = \left(\mu_{i},\sigma_{i}\right)$ be the regime-specific parameters. Given the observed information available at time $t$ denoted as $I_{t-1}$ = $\left\{y_{t-1}, y_{t-2},... \right\}$, the conditional density of $y_t$ is
The Hamilton filter is used evaluate the likelihood recursively and MLE is used to estimate the parameters ($\theta_i$,$\psi$). More detail of estimation can refer to hamilton1994time, Fruhwirth-Schnatter1999State-Space and bazzi2017time.
We assume that transition probabilities vary with time as $\pi_{ij,t}= \frac{exp(-{f_{ij,t}})}{1+\sum^{K}_{j = 1}exp(-{f_{ij,t}})}$ whereas the logistic function is used to map dynamic process $f_{ij,t}$ to the transition probability $f_{ij,t}$ which lies between 0 and 1. Let $f_{t}$ being $K(K-1)\cdot 1$ vector that collects the time-varying parameters $f_{ij,t}, i=1,...,K; j=1,...,K-1$, this paper will utilize following three models to determine the dynamic process $\left\{f_{{t}}, t\geq 0\right\}$:
In the above three models, $\alpha, \beta, \theta, \gamma, w$ are $K(K-1)\cdot 1$ coefficient vectors, $A$ and $B$ are diagonal coefficient matrices. In model (I), $f_{{t}}$ is a linear function of lagged observation $y_{t-1}$. In model (II), $f_{{t}}$ is a linear function of observations from exogenous explanatory variables $\boldsymbol{X}$. Both models were proposed by diebold1994regime and are quite straightforward. Model (III)\footnote{For numerical implementation we use the equivalent mean-reverting representation $f_t = \omega + A\, s_{t-1} + B\,(f_{t-1}-\omega)$ with $\omega = w/(1-B)$; see Section (ref).} is proposed bycreal2013generalized and investigated by bazzi2017time where $f_{{t}}$ is linear function of the scales score process $\left\{s_{{t}}, t\geq 0\right\}$. More specifically, $s_{{t}}$ is the scaled score of the conditional observation density $p (y_{t-1} | \mathcal{I}_{t-1};f_{t-1},\theta)$ with respect to $f_{{t}}$. By setting the scaling matrix $S_{{t}}$ being the square root matrix of the inverse Fisher information matrix, the scaled score function $s_{{t}}$ has a unit variance.
creal2013generalized pointed out as the score depends on the density function of dataset, and it defines the steepest ascent direction for improving the model's local fit in terms of the likelihood at time $t$, thus it is intuitive to use score to update $f_{{t}}$. bazzi2017time gave out also detailed interpretation of how the transition probability can be updated by the conditional observation density. In other words, model (III) can incorporate the information embedded in the conditional observation densities to the dynamics of transition probability. bazzi2017time carried out a comprehensive simulation study and showed that this model can indeed adequately track the dynamic patterns in the transition probabilities, even if the underlying dynamics themselves are possibly misspecified.
All the above three models of $f_{{t}}$ can link eventually the transition probability with informative observations of either own dependent variable or certain exogenous explanatory variables, and determine the time varying patterns of the transition probability. This paper will therefore utilize all the three models of $f_{{t}}$ to capture the dynamics of the transition probability in both the monte-carlo simulation and empirical pricing studies.
To validate the estimation procedure and assess parameter recovery performance under the three TVTP model specifications of Section (ref), we conduct a comprehensive Monte Carlo simulation study. The simulation covers the covariate-driven specifications of diebold1994regime (Models I and II), where transition probabilities depend on lagged observables, and the score-driven framework of creal2013generalized as implemented by bazzi2017time (Model III). While the Monte Carlo study in bazzi2017time was limited to the two-regime case ($K=2$) with a common variance and different means for each regime, our study extends the analysis to three regimes ($K=3$) with regime-specific variances. In addition to correctly specified estimation, we systematically evaluate the estimation robustness to model misspecification by cross-estimating the three TVTP data-generating processes under all four model types (constant, Model I, Model II, and Model III).
We found no publicly available implementation covering Models I--III jointly under the general $K$-regime, regime-specific-variance setting studied here, and therefore developed an open-source R package which provides a unified implementation of all four specifications together with the Monte Carlo infrastructure used to produce the results of this section. The package, which is called multiregimeTVTP, is available at \url{https://github.com/smodee/multiregime-TVTP}.
The package exposes a common interface for the four model families of Section (ref), with each model having matched data<Model>CD() simulators, Rfiltering_<Model>() filters returning the log-likelihood and filtered probabilities, and estimate_<model>_model() estimators that wrap optim with multi-start initialization. All specifications support arbitrary $K \geq 2$ and both the diagonal and off-diagonal parameterizations of the transition matrix. Wald standard errors on the original scale are recovered from the numerical Hessian (numDeriv) via the delta method.
For the Monte Carlo study, higher-level analysis functions coordinate the full DGP $\times$ sample-size $\times$ estimation-model grid, the $R$ replications, and the $n$ random starts per fit, and aggregate replication-level estimates into the bias, RMSE, coverage, forecast, and filtered-probability summaries of Section (ref). Multi-start optimization is parallelized through the future/future.apply backends, so that replications and starts run concurrently across cores. The four filtering routines, that are called at every likelihood evaluation, are additionally backed by a compiled C implementation (src/filtering.c) that yields a 10--56$\times$ speedup per likelihood call relative to the pure-R filter, with the largest gains on Models I and III. The backend is selected automatically when available and can be toggled for benchmarking, which made the full Monte Carlo study tractable (days of compute time instead of weeks/months) on standard desktop hardware while leaving the pure-R filters available as a reference implementation. The same estimation, filtering, and analysis routines are reused without modification in the empirical application of Section (ref).
For Model III we implement the algebraically equivalent mean-reverting form $f_t = \omega + A\, s_{t-1} + B\,(f_{t-1}-\omega)$, where $\omega = w/(1-B)$ is the unconditional mean of $f_t$, following bazzi2017time. The filter is initialized at $f_1 = \omega$, so $A \to 0$ collapses $f_t$ to the constant $\omega$ and the model reduces exactly to constant transition probabilities $\pi_{ij} = \text{logistic}(\omega)$. We exploit this by falling back to the constant filter whenever $\max_{ij}|A_{ij}|$ drops below a small numerical threshold, which also avoids near-flat likelihood regions during optimization.
We consider nine data-generating processes (DGPs), summarised in Table (ref), each evaluated at two sample sizes $T\in\{500,\, 1{,}000\}$, giving 18 scenario--sample-size combinations in total. The DGPs are arranged in two parallel blocks: a $K=2$ block (DGPs 1--4) with diagonal transition parameterization and common variance, and a $K=3$ block (DGPs 6--9) with off-diagonal parameterization and regime-specific variances. Within each block, the first DGP (1 and 6) is a constant-transition baseline, estimated only under the constant specification; the remaining three (2--4 and 7--9) are generated by TVTP Models I, II, and III respectively and estimated under all four model types, yielding the cross-TVTP misspecification analysis. Fitting Models I--III to the constant baselines would only test whether the dynamic coefficients shrink to zero, which is a separate question from the cross-TVTP comparison targeted here. DGP 5 is an additional $K=2$ scenario that keeps Model I as the data-generating process but adopts the off-diagonal, regime-specific-variance parameterization of the $K=3$ block. It is estimated only under the constant and Model I specifications, to isolate the effect of the parameterization change rather than repeat the misspecification comparison already covered by DGPs 2--4.
The true parameter values used to generate the data are reported in Table (ref). Within the $K=2$ block, DGPs 1--4 share regime means $\mu = (-1, 1)$ and a common variance $\sigma^2 = 0.5$; DGP 5 retains the same means but uses regime-specific variances $\sigma^2 = (0.3, 0.7)$. Within the $K=3$ block (DGPs 6--9), all DGPs share regime means $\mu = (-2, 0, 2)$ and regime-specific variances $\sigma^2 = (0.3, 0.5, 0.8)$. The TVTP coefficient magnitudes are moderate, representing empirically plausible effect sizes that allow the transition probabilities to vary meaningfully over time without dominating the regime dynamics.
For DGPs with TVTP dynamics, each dataset is estimated not only under the correctly specified model but also under the remaining model types. This misspecification design, analogous to that of bazzi2017time, allows us to assess how Models I--III perform when the true dynamics follow a different specification. For Model II estimation on data not generated by an exogenous process, we supply an independent standard normal series as the exogenous variable, representing the case where the covariate carries no information about the true regime dynamics.
Each DGP--sample-size combination is replicated $R=50$ times. For every replication, the model is estimated by maximum likelihood using $n=10$ random starting points, keeping the best result (by log-likelihood among converged runs). The filtering procedure uses a burn-in of $B=100$ and a cut-off of $C=10$ observations. To address the label-switching problem inherent in mixture and regime-switching models, we align estimated regimes to the true ordering by finding the permutation of regime labels that minimizes the sum of absolute deviations between estimated and true regime means.
We evaluate estimation performance along four dimensions, following and extending the metrics reported in bazzi2017time.
\paragraph{Parameter recovery.} For each model parameter $\vartheta \in \{\mu_i, \sigma_i^2, \pi_{ij}, A_{ij}, \ldots\}$, with true value $\vartheta_0$, we compute the bias $\text{Bias}(\hat\vartheta) = R^{-1}\sum_{r=1}^{R}(\hat\vartheta_r - \vartheta_0)$ and root mean squared error $\text{RMSE}(\hat\vartheta) = \bigl[R^{-1}\sum_{r=1}^{R}(\hat\vartheta_r - \vartheta_0)^2\bigr]^{1/2}$ across replications.
\paragraph{Coverage rates.} Standard errors are obtained from the numerical Hessian of the log-likelihood (evaluated at the MLE in the transformed parameter space) via the numDeriv package, with the delta method applied for parameters subject to log or logit transformations. We report the empirical coverage rate of nominal 95% Wald confidence intervals.
\paragraph{Forecast precision.} Following bazzi2017time, we compute the one-step-ahead conditional forecasts $\hat{y}_{t|t-1} = \sum_{i=1}^{K} \hat\mu_i \, \hat{P}(z_t = i \mid I_{t-1})$ and evaluate them using the mean absolute forecast error (MAFE), mean squared forecast error (MSFE), and their standardized counterparts MASFE and MSSFE, where standardization is by the conditional standard deviation $\hat\sigma_{t|t-1}$.
\paragraph{Filtered probability accuracy.} For correctly specified estimations where the true filtered probabilities can be recovered by running the filter at the true parameter values, we compute the mean squared error and mean absolute error of $\hat\pi_{ij,t}$ relative to $\pi_{ij,t}^{*}$, averaged over all transition elements and time periods.
We first present parameter recovery and coverage results for correctly specified models, then compare forecast precision across estimation models to assess misspecification robustness.
Table (ref) reports average bias and RMSE for each parameter group under correct specification, with the true parameter values given in Table (ref). For DGPs 1 and 6 (constant transition probabilities), the regime means $\mu_i$, variances $\sigma_i^2$, and transition probabilities $\pi_{ij}$ are all recovered with small bias and RMSE that decrease as the sample size doubles from $T=500$ to $T=1{,}000$, confirming consistency. For Model I (DGPs 2, 5, 7), the distribution parameters $\mu_i$ and $\sigma_i^2$ are similarly well recovered, but the TVTP coefficient vector $A$ exhibits substantially larger RMSE, particularly for $K=3$ regimes. At $T=500$, 4 of 44 convergent replications in DGP 7 exhibit extreme $A$ estimates ($|\hat{A}_{ij}|>10$, against true values of order $0.05$), indicating optimizer breakdown on a flat likelihood region rather than genuine poor estimation. These replications are excluded from the reported $A$ statistics; after exclusion the trimmed RMSE is $1.550$, and increasing the sample to $T=1{,}000$ reduces it further to $0.549$. The off-diagonal parameterization of DGP 5 also yields elevated $A$ RMSE at $T=500$ ($1.228$), which improves to $0.319$ at $T=1{,}000$.
For Model II (DGPs 3 and 8), parameter recovery shows a similar pattern: the $\gamma$ coefficients (reported in the $A$ column) have larger RMSE than the distribution parameters, with clear improvement at $T=1{,}000$. In the two-regime setting (DGP 3), the $A$ RMSE drops from $0.311$ to $0.189$; in the three-regime setting (DGP 8), it drops from $0.815$ to $0.298$.
For Model III (DGPs 4 and 9), the distribution parameters $\mu_i$ and $\sigma_i^2$ and the transition probabilities $\pi_{ij}$ are well recovered, with bias and RMSE comparable to the other TVTP models. However, the score coefficient matrix $A$ is not reliably identified. The RMSE equals the magnitude of the true $A$ values and shows no improvement from $T=500$ to $T=1{,}000$ in either DGP, pointing to a statistical identifiability problem inherent to the GAS specification rather than a computational artefact. Profile likelihood analysis confirms this, showing that the 1D NLL is lower at $A=0$ than at the true value and that the MLE is at $A\approx 0$ regardless of the data-generating parameter. Examining the joint $(\sigma^2, A)$ space through a 2D profile further reveals a pronounced ridge in the likelihood surface, arising because the GAS score scaling couples $\sigma^2$ (through the Fisher information) and $A$ (as a direct multiplier) such that many parameter combinations yield nearly identical filtered transition probabilities and log-likelihoods.
Table (ref) reports the average empirical coverage rates of nominal 95% Wald confidence intervals. For the constant-probability DGPs (1 and 6), coverage rates for $\mu$ and $\sigma^2$ are close to the nominal level across both sample sizes, though $\sigma^2$ coverage in the three-regime case (DGP 6) drops to $0.867$ at $T=500$ before recovering to $0.960$ at $T=1{,}000$.
Under TVTP specifications, the coverage for transition probability parameters $\pi$ is systematically below the nominal 95% level, particularly for $K=3$ models. In DGPs 7 and 8 at $T=500$, the average $\pi$ coverage is only $0.675$ and $0.713$, respectively. While coverage improves at $T=1{,}000$ ($0.753$ and $0.817$), it remains substantially below nominal. Coverage for $\mu$ remains robust across all settings, generally exceeding $0.90$. For Model III (DGP 4), $\sigma^2$ coverage reaches $1.000$ at $T=500$, consistent with a numerically ill-conditioned Hessian along the $(\sigma^2, A)$ ridge identified in the parameter recovery discussion above: the near-flat likelihood inflates the computed standard errors, producing confidence intervals wide enough to contain the true value in virtually every replication.
Table (ref) compares forecast precision across estimation models for each DGP, assessing both correct-specification performance and robustness to misspecification. Across all DGPs and sample sizes, the forecast metrics are remarkably stable: MAFE and MSFE differ by less than 1% whether the correctly specified or a misspecified model is fitted. This finding extends the result of bazzi2017time from two regimes to three, and its explanation is the same. The one-step-ahead forecast $\hat{y}_{t|t-1} = \sum_i \hat{\mu}_i \hat{P}(z_t = i \mid I_{t-1})$ is dominated by the regime means $\hat{\mu}_i$, which are robustly recovered under all specifications. Since transition dynamics have only a marginal effect on the predicted probabilities over a single step, misspecifying them has negligible impact on point forecast accuracy.
This result should be interpreted with the scope of the evaluation in mind. The study considers only one-step-ahead point forecasts, which is precisely the setting where transition dynamics matter least. At longer horizons the transition matrix is applied repeatedly and misspecification would likely compound into visible differences. Density and interval forecasts would also discriminate more sharply between specifications, since a correctly specified TVTP model captures time-varying uncertainty that a misspecified one cannot. The practical value of correct specification is more directly visible in the filtered probability accuracy reported in Table (ref), where differences across models are substantial.
Table (ref) reports the accuracy of filtered transition probabilities for correctly specified models, measured by mean squared error (MSE) and mean absolute error (MAE) relative to the true filtered probabilities. For the $K=2$ constant model (DGP 1), filtered probability recovery is excellent (MSE $< 0.001$), and increases only modestly for the TVTP specification of DGP 2.
Moving to $K=3$, the filtered probability accuracy deteriorates substantially: MSE values range from $0.115$ (DGP 6, $T=500$) to $0.180$ (DGP 7, $T=500$), reflecting the inherent difficulty of distinguishing among three regimes. The off-diagonal transition parameterization in the $K=3$ block implies $K(K-1) = 6$ time-varying parameters per time step, considerably increasing the estimation burden. Per-regime analysis (not in the table) reveals that the middle regime is consistently the hardest to identify, with approximately twice the MSE of the extreme regimes. In contrast to the point forecast results, where all specifications perform equivalently, filtered probability accuracy is where the choice of TVTP specification has the most tangible effect: correctly specifying the transition dynamics yields meaningfully more accurate regime probability estimates, particularly in the three-regime setting.
The Monte Carlo results demonstrate that the estimation procedure reliably recovers the distribution parameters ($\mu_i$, $\sigma_i^2$) and transition probabilities across both $K=2$ and $K=3$ settings, with performance improving at larger sample sizes as expected. The TVTP driving coefficients ($A$ for Models I and II) are the most challenging parameters to estimate, requiring $T=1{,}000$ observations for the three-regime case to achieve acceptable RMSE levels. For Model III (GAS), the score coefficient $A$ is statistically non-identifiable due to a ridge in the joint $(\sigma^2, A)$ likelihood surface; the remaining parameters are nevertheless well recovered. Coverage rates for transition probabilities under TVTP models are systematically below the nominal 95%, particularly for $K=3$, indicating that alternative standard error procedures (e.g., bootstrap or sandwich estimators) may be warranted in practice.
The misspecification analysis shows that one-step-ahead point forecast accuracy is robust to the choice of TVTP specification, consistent with bazzi2017time. This robustness is a mechanical consequence of regime means dominating short-horizon forecasts, and should not be read as implying that model choice is inconsequential. The filtered probability results show meaningful accuracy differences across specifications, and longer-horizon or density-based forecast evaluations would be expected to further differentiate the models. The primary benefit of correct TVTP specification lies in accurately characterising the time-varying regime dynamics, not in improving point forecast accuracy.
We apply the three-regime TVTP models to U.S.\ Treasury zero-coupon yield data reconstructed by liu2021reconstructing. The dataset contains monthly annualized continuously-compounded zero-coupon yields for maturities ranging from 1 to 360 months. Our sample spans June 1961 to December 2024, providing $T = 763$ monthly observations. We select four representative maturities along the yield curve: 1 month (short end), 12 months (short-to-medium), 36 months (medium), and 72 months (long end).
Following standard practice for yield curve modelling, we work with first differences of the yield series, i.e.\ $y_t = Y_t - Y_{t-1}$, where $Y_t$ denotes the yield level at time $t$. This transformation produces $T = 762$ observations of monthly yield changes for each maturity. First-differencing removes the strong persistence in yield levels hamilton1994time and produces a more nearly stationary series suitable for regime-switching analysis.
For each maturity, we estimate a three-regime ($K=3$) Markov switching model where the conditional distribution of yield changes in regime $i$ is
with regime-specific means $\mu_i$ and variances $\sigma_i^2$. The three regimes are intended to capture distinct yield curve dynamics: a high-volatility regime associated with large yield movements, a moderate regime representing normal market conditions, and a low-volatility regime reflecting periods of relative stability.
We consider four model specifications for the transition probability dynamics:
The transition probabilities are parameterized via the logistic link function as described in Section (ref), with $K(K-1) = 6$ unconstrained parameters governing the off-diagonal elements of the transition matrix. Each model is estimated by maximum likelihood using $n = 100$ random starting points with a burn-in period of $B = 100$ observations. The best result across starting points (by log-likelihood among converged runs) is retained.
We attempted to estimate the GAS specification (Model III) alongside the other models. However, across all four maturities, the GAS model exhibited severe convergence difficulties: out of 100 random starting points per maturity, only a single start converged in each case, and in every instance the estimated score coefficients $A_{ij}$ collapsed to exactly zero, reducing the model to the constant-probability specification. The remaining 99 starts with non-zero initial $A_{ij}$ values uniformly failed to converge, suggesting that the GAS likelihood surface is essentially flat or ill-conditioned in the score-coefficient dimensions for these yield curve series. These symptoms are consistent with the identifiability issue diagnosed in the Monte Carlo study (Section (ref)), where a ridge in the joint $(\sigma^2, A)$ likelihood drove $\hat A$ to zero regardless of the true parameter. We therefore exclude the GAS specification from the results tables below and focus our comparison on the Constant, TVP, and Exogenous models.
Table (ref) reports the model fit statistics for each maturity and model specification. Table (ref) presents the estimated regime-specific parameters.
Figure (ref) displays the filtered regime classifications from the Exogenous model across all four maturities. At each point in time, the background shading indicates the most probable regime, with regimes ordered by variance: blue corresponds to the low-volatility regime, salmon to moderate volatility, and red to the high-volatility regime.
The paper studies time-varying transition probability (TVTP) Markov-switching (MS) models with various dynamics in the transition probability matrix. The regime switching models are univariate Markov-switching models with $K$ regimes, each with its own mean and variance and regime changes follow a first‑order Markov chain. A unified R package, multiregimeTVTP, is developed to simulate data, and estimate the model parameters for arbitrary $K\geq 2$. Monte Carlo simulation shows that regime means/variances and average transition probabilities are well estimated, but TVTP coefficients---especially in GAS models and three‑regime settings---are hard to identify. Short‑horizon point forecasts are robust to TVTP misspecification, yet regime probability estimates and empirical fit, particularly for yield curves, clearly benefit from correctly modeling the time‑varying, covariate‑driven transition probabilities. Apart from more empirical applications of the TVTP-MS model with more than two regimes, further analysis for different processes governing the transition probabilities and the price process can be interesting extension. For example, a mean-reverting process might be more suitable than a GBM for firms connected to commodity markets, or transition probabilities could be cyclical (Bazzi {\it et al.}, 2017). Also, the impact of the specific functional form for calculating the transition probabilities remains an interesting question for future work.