Keywords— Markov-switching model, Time-varying transition probabilities, GAS model, Monte Carlo Simulation, Zero-coupon yield

1 Introduction↩︎

To capture nonlinear regime-dependent and cyclical dynamics in economics and finance, [2] 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 [3],[4], [5], [6], [7],[8],[9],[10]. 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. [11] extended the [2] model and allowed for transition probabilities depending on other economic indicators, making the regime change responsive to economic behavior. [12] and [13] 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. [12] 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 [11] and [12], the driver of the time-‐varying transition matrix are exogenous variables. [14] introduced a framework to update the parameters of the transition probability matrix based on the predictive likelihood score. [1] adopted the framework of [14]’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 [14] and [1] 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 [15] and the regime-switching models proposed by [16] 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 [17]; specify the volatility of crude oil futures prices[18]; document regime‑dependent relationships between output, money, and prices [19]; or perform the industrial electricity load forecasting [20].

Among different TVTP-MS models, the [1] 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 [1] 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 [21].

The paper is divided into follow sections: section 2 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 3 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 4 is the empirical analysis and the economic interpretation of the result. Section 5 is conclusion and suggestions for future research.

2 Regime-switching model with time-varying transition probabilities↩︎

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 \[\begin{align} \mathbb{P}\left(y_t|\psi,I_{t-1}\right) &=& \sum^{K}_{i = 1}\sum^{K}_{k = 1}\mathbb{P}(y_t|\theta_i,\psi)\pi_{ki}\mathbb{P}(z_{t-1} = k|\psi,I_{k-1})\nonumber\\ \end{align}\] 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 [22], [23] and [1].

Three types of TVTP models↩︎

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\}\): \[\begin{align} \tag{1} & {\rm Model\; (I):} & f_{t}=\alpha+\theta y_{t-1} \\ & {\rm Model\; (II):} & f_{t}=\beta+ \gamma X_{t-1}\tag{2} \\ & {\rm Model\; (III):} & f_{t}=w+As_{t-1}+Bf_{t-1}; \tag{3} \nonumber\\ && s_{t-1} = S_{t-1}\bigtriangledown_{t-1}, \bigtriangledown_{t-1}= \frac{\partial}{\partial f_{t-1}}\log p(y_{t-1} | \mathcal{I}_{t-1};f_{t-1},\theta) \end{align}\]

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 [11] and are quite straightforward. Model (III)¹ is proposed by[14] and investigated by [1] 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.

[14] 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}}\). [1] 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. [1] 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.

3 Monte Carlo simulation study↩︎

To validate the estimation procedure and assess parameter recovery performance under the three TVTP model specifications of Section ¿sec:sec:three32models?, we conduct a comprehensive Monte Carlo simulation study. The simulation covers the covariate-driven specifications of [11] (Models I and II), where transition probabilities depend on lagged observables, and the score-driven framework of [14] as implemented by [1] (Model III). While the Monte Carlo study in [1] 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).

3.1 Software implementation↩︎

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 https://github.com/smodee/multiregime-TVTP.

The package exposes a common interface for the four model families of Section ¿sec:sec:three32models?, 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 3.3. 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 4.

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 [1]. 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.

3.2 Simulation design↩︎

We consider nine data-generating processes (DGPs), summarised in Table 1, 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 2. 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 [1], 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.

3.3 Performance metrics↩︎

We evaluate estimation performance along four dimensions, following and extending the metrics reported in [1].

3.3.0.1 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.

3.3.0.2 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.

3.3.0.3 Forecast precision.

Following [1], 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}\).

3.3.0.4 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.

3.4 Simulation results↩︎

We first present parameter recovery and coverage results for correctly specified models, then compare forecast precision across estimation models to assess misspecification robustness.

Parameter recovery↩︎

Table 3 reports average bias and RMSE for each parameter group under correct specification, with the true parameter values given in Table 2. 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.

Coverage rates↩︎

Table 4 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.

Forecast precision and misspecification robustness↩︎

Table 5 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 [1] 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 6, where differences across models are substantial.

Filtered probability accuracy↩︎

Table 6 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.

Summary of Monte Carlo simulation↩︎

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 [1]. 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.

4 Empirical analysis↩︎

4.1 Data↩︎

We apply the three-regime TVTP models to U.S.Treasury zero-coupon yield data reconstructed by [21]. 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 [22] and produces a more nearly stationary series suitable for regime-switching analysis.

4.2 Model specification↩︎

For each maturity, we estimate a three-regime (\(K=3\)) Markov switching model where the conditional distribution of yield changes in regime \(i\) is \[(y_t \mid z_t = i) \sim \mathcal{N}(\mu_i, \sigma_i^2), \quad i = 1, 2, 3,\] 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:

• Constant: The baseline model with time-invariant transition probabilities, corresponding to the classical Hamilton framework ([2]) extended to three regimes.

• Model (I) – TVP: Time-varying transition probabilities driven by the lagged observation \(y_{t-1}\), as in equation (1 ). This specification allows the most recent yield change to inform the probability of transitioning between regimes.

• Model (II) – Exogenous: Time-varying transition probabilities driven by an exogenous covariate \(X_{t-1}\), as in equation (2 ). Here we set \(X_{t-1} = Y_{t-1}\), the lagged yield level, so that regime dynamics are linked directly to the level of interest rates rather than the yield change used in the TVP specification.

• Model (III) – GAS: Generalized Autoregressive Score transition probabilities, as in equation (3 ), where the transition parameters are updated using the scaled score of the conditional observation density.

The transition probabilities are parameterized via the logistic link function as described in Section ¿sec:sec:three32models?, 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.

4.3 Results↩︎

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 3.4), 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 7 reports the model fit statistics for each maturity and model specification. Table 8 presents the estimated regime-specific parameters.

Figure 1 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.

5 Conclusions and future work↩︎

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 et al., 2017). Also, the impact of the specific functional form for calculating the transition probabilities remains an interesting question for future work.

References↩︎