Nonlinear Modeling via Discrete Regime Change
2026-07-24
We consider:
Both settings allow for unobserved regimes and dynamic specifications (lagged dependent variables, autocorrelated errors).
Let y_t follow a process governed by an unobserved discrete state variable s_t \in \{1, \ldots, M\}.
Each regime m has its own regression model. Given regressors X_t (switching coefficients) and Z_t (common coefficients):
\mu_{t,m} = X_t'\beta_m + Z_t'\gamma
Errors are normally distributed with (potentially) regime-specific variance:
y_t \mid s_t = m \;\sim\; \mathcal{N}\!\left(\mu_{t,m},\; \sigma_m^2\right)
The likelihood contribution for observation t weights each regime density by its one-step-ahead probability:
L_t = \sum_{m=1}^{M} \frac{1}{\sigma_m} f\!\left(\frac{y_t - \mu_{t,m}}{\sigma_m}\right) P(s_t = m \mid \mathcal{F}_{t-1};\, \delta)
The full log-likelihood is:
\ell(\beta, \gamma, \sigma, \delta) = \sum_{t=1}^{T} \log L_t
Maximized via iterative numerical methods — BFGS is standard across implementations.
In simple switching, regime probabilities are independent across time (no state dependence).
Probabilities are parameterized via a multinomial logit:
P(s_t = m \mid G_{t-1};\, \delta) = \frac{\exp(G_{t-1}'\delta_m)}{\sum_{j=1}^{M} \exp(G_{t-1}'\delta_j)}
with the normalization \delta_M = 0.
Note: Valid only for models with uncorrelated errors. Serial correlation requires a different approach.
After observing y_t, we can update our regime probability estimates using Bayes’ theorem:
P(s_t = m \mid \mathcal{F}_t) = \frac{\frac{1}{\sigma_m}f\!\left(\frac{y_t - \mu_{t,m}}{\sigma_m}\right) P(s_t = m \mid \mathcal{F}_{t-1})}{\sum_{j=1}^{M} \frac{1}{\sigma_j} f\!\left(\frac{y_t - \mu_{t,j}}{\sigma_j}\right) P(s_t = j \mid \mathcal{F}_{t-1})}
For simple switching without AR terms, filtered = smoothed.
Under Markov switching, the regime probability depends on the previous state:
P(s_t = j \mid s_{t-1} = i) = p_{ij,t}
Collected into an M \times M transition matrix:
\mathbf{P}_t = \begin{pmatrix} p_{11,t} & \cdots & p_{1M,t} \\ \vdots & \ddots & \vdots \\ p_{M1,t} & \cdots & p_{MM,t} \end{pmatrix}
The (i,j) element is the probability of transitioning from regime i at t-1 to regime j at t.
Each row of \mathbf{P}_t is parameterized by a separate multinomial logit:
p_{ij,t} = P(s_t = j \mid s_{t-1} = i;\, G_{t-1}, \delta_i) = \frac{\exp(G_{t-1}'\delta_{ij})}{\sum_{s=1}^{M}\exp(G_{t-1}'\delta_{is})}
with \delta_{iM} = 0 for identification.
Special cases:
| Case | Description |
|---|---|
| G_{t-1} \equiv 1 | Constant transition probs (Hamilton 1989) |
| G_{t-1} = leading indicator | Time-varying probs (Filardo 1994; Diebold et al. 1994) |
Starting from filtered probabilities P(s_{t-1} = m \mid \mathcal{F}_{t-1}), each period involves four steps:
1. Predict one-step-ahead regime probabilities: P(s_t = m \mid \mathcal{F}_{t-1}) = \sum_{j=1}^{M} p_{jm,t}\, P(s_{t-1} = j \mid \mathcal{F}_{t-1})
2. Form joint densities of data and state: f(y_t, s_t = m \mid \mathcal{F}_{t-1}) = \frac{1}{\sigma_m} f\!\left(\frac{y_t - \mu_{t,m}}{\sigma_m}\right) P(s_t = m \mid \mathcal{F}_{t-1})
3. Evaluate the likelihood contribution (marginalise over states): L_t = \sum_{j=1}^{M} f(y_t, s_t = j \mid \mathcal{F}_{t-1})
4. Filter — update regime probabilities via Bayes’ rule: P(s_t = m \mid \mathcal{F}_t) = \frac{f(y_t, s_t = m \mid \mathcal{F}_{t-1})}{L_t}
Smoothed probabilities use the full-sample information set \mathcal{F}_T rather than only \mathcal{F}_t.
Kim’s (1994) efficient backward recursion:
P(s_t = i, s_{t+1} = j \mid \mathcal{F}_T) = P(s_t = i \mid \mathcal{F}_t) \cdot \frac{p_{ij,t+1}\, P(s_{t+1} = j \mid \mathcal{F}_T)}{P(s_{t+1} = j \mid \mathcal{F}_t)}
Marginalize over s_{t+1} to obtain smoothed probabilities:
P(s_t = i \mid \mathcal{F}_T) = \sum_{j=1}^{M} P(s_t = i, s_{t+1} = j \mid \mathcal{F}_T)
The Markov filter requires initialization P(s_0 = m \mid \mathcal{F}_0). Common approaches:
| Method | Description | Default in |
|---|---|---|
| Ergodic (steady-state) | Implied by transition matrix | Stata, most packages |
| Uniform | 1/M for each regime | Manual override |
| User-specified | Analyst-supplied prior | All packages |
| Estimated | Treated as free parameters | depmixS4 |
Ergodic initialization is arbitrary when transition probabilities are time-varying.
Two strategies for introducing dynamics into switching models:
| MSDR | MSAR / Hamilton Model | |
|---|---|---|
| Full name | Markov Switching Dynamic Regression | Markov Switching Autoregressive |
| Also known as | MSI (regime-invariant lags) | MSM, Hamilton model |
| Dynamics via | Lagged y_t as regressors | Serially correlated errors |
| Mean depends on | Current state s_t only | Current and lagged states |
| State vector | M-dimensional | M^{p+1}-dimensional |
| Standard filter | Applies directly | Requires Hamilton filter |
Include lagged dependent variables directly as regressors:
y_t = X_t'\beta_{s_t} + Z_t'\gamma + \sum_{r=1}^{p} \phi_r^{(s_t)}\, y_{t-r} + \sigma_{s_t} \varepsilon_t
where X_t'\beta_{s_t} + Z_t'\gamma is the regime mean from regressors and \phi_r^{(s_t)} are the (optionally regime-varying) lag coefficients.
Deviations of y_t from the regime mean follow an AR(p) process:
\left(1 - \sum_{r=1}^{p} \rho_r L^r\right)(y_t - \mu_{s_t}) = \sigma_{s_t} \varepsilon_t
Rearranged (expanding the lag operator):
y_t = \mu_{s_t} + \sum_{r=1}^{p} \rho_r\left(y_{t-r} - \mu_{s_{t-r}}\right) + \sigma_{s_t} \varepsilon_t
Key points:
mswitchmswitch has been part of official Stata since Stata 14. Two sub-commands map directly onto the two dynamic specifications:
| Command | Model | Dynamics |
|---|---|---|
mswitch dr |
MSDR | Lagged y_t as regressors |
mswitch ar |
MSAR | AR errors (Hamilton filter) |
norswitching(varlist) specifies regime-invariant regressors Z_t.
| Option | Purpose |
|---|---|
states(#) |
Number of regimes M |
ar(numlist) |
AR lags for MSAR (e.g., ar(1/4)) |
norswitching(varlist) |
Regime-invariant regressors Z_t |
varswitch |
Regime-specific error variances \sigma_m^2 |
tvtp(varlist) |
Time-varying transition probabilities |
from(matrix) |
User-supplied starting values |
technique(bfgs) |
Optimization algorithm (default BFGS) |
Use
tvtp(c l.dlogidx)for Filardo-style TVP: Stata’sl.prefix lags within the option, consistent with the data being at period t.
* Transition matrix and expected durations
estat transition
* Regime probabilities (all three types available)
predict p1 p2, pr // one-step-ahead
predict f1 f2, filterpr // filtered
predict s1 s2, smoothpr // smoothed ← Kim smoother
* Plot smoothed P(recession)
tsline s2, yline(0.5, lpattern(dash)) ///
title("Smoothed P(Regime 2)") ///
xtitle("") ytitle("Probability")
* Transition probability restrictions: supply a matrix
* with fixed values; missing (.) means estimated
matrix RESTR = (., 0 \ ., .) // p_12 = 0 forced
mswitch ar y, states(2) ar(1) constraint(RESTR)* Data: quarterly log-difference of US GNP, 1951q1–1984q4
use gnp_hamilton.dta, clear
tsset date
* Estimate MSAR(4), 2 states, common AR and variance
mswitch ar g, states(2) ar(1/4)
* Inspect
estat transition // transition matrix + expected durations
estimates store msar4
* Extract and plot smoothed recession probability
predict smprob1 smprob2, smoothpr
tsline smprob2, ///
title("P(S(t)=2 | Full Sample)") ///
ytitle("Probability") xtitle("")Expected durations: ~10.4 quarters in expansion, ~4.1 in recession.
Three packages cover most switching regression needs:
| Package | Strengths | Typical use |
|---|---|---|
MSwM |
Simple API; wraps lm/glm |
Quick MSDR / MSAR |
depmixS4 |
EM algorithm; TVP; non-Gaussian | Flexible HMM / MSAR |
MSGARCH |
Regime-switching GARCH | Financial volatility |
Neither package replicates every feature from the chapter out of the box. For production work, check each package’s vignette and verify the log-likelihood against a known benchmark.
MSwM — MSAR WorkflowmsmFit() wraps a fitted lm object. The p argument adds AR lags (MSAR); sw controls which parameters are regime-specific.
library(MSwM)
# Baseline model (intercept only — regime switching applies to mean)
base <- lm(g ~ 1, data = gnp_df)
# 2-state MSAR(4): intercept switches, AR coefs and sigma common
# sw: logical vector over (intercept, sigma) — length = npar + 1
ms <- msmFit(base,
k = 2, # number of regimes
p = 4, # AR order (MSAR)
sw = c(TRUE, FALSE), # intercept switches; sigma does not
control = list(parallel = FALSE))
summary(ms)ms@transMat — transition matrix; ms@Coef — regime-specific coefficients.
MSwM — Regime Probabilities & Plots# Filtered probabilities: T × k matrix
filt <- ms@Fit@filtProb # P(s_t = k | F_{t-1})
# Smoothed probabilities
smoo <- ms@Fit@smoProb # P(s_t = k | F_T)
# Built-in plot
plotProb(ms, which = 2) # smoothed P(regime 2)
# Or with ggplot2
library(ggplot2)
df_prob <- data.frame(date = gnp_df$date,
p_rec = smoo[, 2])
ggplot(df_prob, aes(date, p_rec)) +
geom_line(colour = "steelblue") +
geom_hline(yintercept = 0.5, linetype = "dashed") +
labs(title = "Smoothed P(Recession)", x = NULL, y = "Prob.") +
theme_minimal()Add NBER shading with geom_rect() using dates from the {usrecessions} package.
depmixS4 — Flexible HMMdepmixS4 is better suited to TVP models, non-Gaussian responses, and models with covariates on the transition matrix.
library(depmixS4)
# Specify: Gaussian response, 2 states, TVP via transition covariate
mod <- depmix(
response = g ~ 1, # switching mean
data = gnp_df,
nstates = 2,
family = gaussian(),
transition = ~ dlogidx # TVP: Filardo specification
)
# Fit via EM (try multiple starts to avoid local optima)
set.seed(42)
fit <- fit(mod, verbose = FALSE)
summary(fit)
# Posterior state probabilities
post <- posterior(fit) # columns: state, S1, S2, ...getpars(fit) extracts all parameters; logLik(fit) gives the maximised log-likelihood.
Switching models often have multiple local optima. Strategies in R:
# Run msmFit from multiple random starts (MSwM)
best <- NULL
for (i in 1:25) {
set.seed(i)
try({
m <- msmFit(base, k = 2, p = 4, sw = c(TRUE, FALSE),
control = list(parallel = FALSE))
if (is.null(best) || logLik(m) > logLik(best)) best <- m
})
}
# depmixS4: EM with random initialisation
fits <- replicate(25, {
set.seed(sample.int(1e5, 1))
try(fit(mod, verbose = FALSE), silent = TRUE)
}, simplify = FALSE)
best_depmix <- fits[[which.max(sapply(fits, function(x)
tryCatch(logLik(x), error = function(e) -Inf)))]]| Feature | Stata mswitch |
R MSwM |
R depmixS4 |
|---|---|---|---|
| MSDR | ✓ | ✓ | ✓ |
| MSAR | ✓ | ✓ | Limited |
| TVP transitions | ✓ (tvtp()) |
✗ | ✓ (transition=~) |
| Regime-specific \sigma | ✓ (varswitch) |
✓ | ✓ |
| > 2 regimes | ✓ | ✓ | ✓ |
| Transition restrictions | ✓ | ✗ | Partial |
| Smoothed probabilities | ✓ | ✓ | ✓ (posterior) |
| Non-Gaussian response | ✗ | Limited | ✓ |
| Multiple random starts | Manual | Manual | Manual |
| Publication-ready plots | Via coefplot |
Via ggplot2 |
Via ggplot2 |
Model: Two-state Markov switching mean for U.S. GNP growth, with regime-invariant AR(4) errors.
Specification:
Regime-specific means — Hamilton’s “fast” and “slow” growth:
| Regime 1 (Expansion) | Regime 2 (Recession) | |
|---|---|---|
| C | 1.164 (0.000) | −0.359 (0.189) |
Transition matrix:
\hat{\mathbf{P}} = \begin{pmatrix} 0.904 & 0.096 \\ 0.245 & 0.755 \end{pmatrix}
Expected durations: ~10.4 quarters in expansion, ~4.1 quarters in recession.
Predicted probabilities of the low-output state align closely with NBER recession dates.
Model: MSAR(4) for industrial production (DLOGIP) with time-varying transition probabilities driven by a leading indicator (DLOGIDX).
Specification: 1948M01–1991M04; switching mean, common AR(1)–AR(4) and \sigma; TVP logit on G_{t-1} = \Delta \log \text{IDX}_{t-1}.
Use the lagged indicator so that the period-t transition probability depends on information available at t-1.
Regime-specific means:
| Regime 1 (Expansion) | Regime 2 (Contraction) | |
|---|---|---|
| C | 0.517 (0.000) | −0.866 (0.000) |
Time-varying transition parameters:
| Parameter | Estimate | Interpretation |
|---|---|---|
| P11-C | 4.36 | Baseline persistence in expansion |
| P11-DLOGIDX(−1) | 1.77 | Higher growth → higher prob. of staying in expansion |
| P21-C | −1.65 | Baseline exit from contraction |
| P21-DLOGIDX(−1) | 0.99 | Higher growth → higher prob. of moving to expansion |
Model: Three-state Markov switching model of stock return volatility — mean zero, regime-specific variances \sigma_1^2 < \sigma_3^2 < \sigma_2^2.
Specification: Monthly CRSP excess returns, 1926M01–1986M12; intercept fixed at zero; 3 regimes; constant transition probabilities.
Estimation challenge: Many local roots — run at least 50 random starts and select the solution with the highest log-likelihood.
Implied standard deviations by regime:
| Regime | \log(\hat{\sigma}) | \hat{\sigma} | Interpretation |
|---|---|---|---|
| 1 | −3.352 | 0.035 | Low volatility |
| 2 | −1.736 | 0.176 | High volatility |
| 3 | −2.760 | 0.063 | Medium volatility |
Expected durations: ~38 months (low), ~20 months (high), ~32 months (medium).
Key finding: Transition probability from low to high volatility is essentially zero — warranting a restriction p_{12} = 0.
Impose via a 3 \times 3 restriction matrix with 0 in cell (1,2) and NA elsewhere (except fixed cells).
mswitch) and R (MSwM, depmixS4) provide mature implementations; multiple random starts are essential in both to escape local optima
ESSSSDA26-2J: One or More Time Series