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.
86,906 characters · 24 sections · 42 citation commands
Regime-Switching Models for Disaggregated Data
\\[1em]
}
\thispagestyle{empty}\setcounter{page}{0}
Regime-switching models have become a common framework for capturing nonlinear dynamics in macroeconomic and financial data since hamilton1989.\footnote{See Ang2012 and hamilton2016 for surveys.} Most applications fit such models to just one series. For example, hamilton1989 and hamilton2006 studied U.S. aggregate output; masson2001 examined exchange rates; hamilton1988, hamilton2005, and Ang2002b modeled interest rates; davig2004 analyzed the debt-output ratio; and Ang2002a considered equity returns. A small number of studies looked at small-scale multivariate settings: see diebold1996, kim1998, chauvet1998, and hamilton2006, who modeled sales, income, employment, and industrial production using a latent factor with regime switching; Chauvet2002, who applied a similar model to U.S. unemployment rates across seven age groups; and krolzig2000, who estimated a Markov-switching vector error correction model, for six industrial production indices for U.K. business cycles.
Even for datasets with many series, most analyses treat series separately (e.g., OwyangPigerWall2005) rather than formulate a joint model. Only a few papers have studied high-dimensional models: LiuChen2016 propose a factor model with regime-dependent means, loadings, and covariance structure, estimated via eigenanalysis and the Viterbi algorithm. UrgaWang2024 develop a quasi-maximum likelihood estimator for factor models with regime-switching loadings, combining an EM algorithm with principal component analysis. BarigozziMassacci2025 introduce an approximate factor model with regime-switching loadings driven by a latent Markov process, estimated by principal components followed by an EM step. In a different approach, HamiltonOwyang2012 develop a Bayesian latent-class model for state-level recessions, modeling recession probabilities through logistic functions tied to national and regional clusters; their results show cross-state differences in recession timing but little evidence that states avoid national downturns. These contributions broaden the scope of regime-switching analysis and suggest potential benefits can come from modeling many series simultaneously.
This study considers regime switching in high-dimensional settings, with a different motivation from prior work. Our starting point is the observation that many macroeconomic variables are built from disaggregated components. For example, GDP combines personal consumption expenditures, business investment, government spending, and imports and exports, each containing subcomponents.\footnote{Source: U.S. Bureau of Economic Analysis, via \href{https://fred.stlouisfed.org/release/tables?rid=53&eid=13498\#snid=13511} {Federal Reserve Economic Data (FRED)}, Federal Reserve Bank of St.\ Louis. Percent Change From Preceding Period in Real Gross Domestic Product, Expanded Detail. } The same applies to industrial production, capacity utilization, employment, and hours worked, among many others. In practice, one can work with either aggregate or disaggregate data when modeling regime switching. This raises the question of whether disaggregate data can better identify regime switching than the corresponding aggregates. We study this question analytically and through simulations; both suggest that disaggregate data can lead to sharper identification, with gains increasing in the heterogeneity of signal-to-noise ratios across series. Given the prevalence of aggregation in economic data, this suggests that models and estimation methods designed for disaggregated series can play a useful role and should not be overlooked in practice.
To illustrate this point empirically, we apply a Hamilton (1989) style model to U.S. quarterly real GDP growth from 1972:II to 2019:III. The sample is restricted to the pre-COVID period, so that the pandemic does not dominate the comparison. Let $Y_t$ denote GDP growth. We estimate
where $S_t$ follows a two-state Markov chain, with $S_t=1$ and $S_t=2$ being recession and expansion, respectively; $\mu_2$ is the expansion mean; $\delta_1$ is the difference between the recession and expansion means; and $e_t \sim \text{i.i.d. } N(0,\sigma_1^2)$ if $S_t=1$ and $e_t \sim \text{i.i.d. } N(0,\sigma_2^2)$ if $S_t=2$. Figure (ref)(a) shows the smoothed recession probabilities, with NBER-defined recessions shown as shaded bars. The recession probabilities are broadly consistent with the NBER chronology; however, the turning points are often not sharply identified, and the model essentially misses the recession in the early 2000s. These patterns also appear in the recession probabilities published by the St.\ Louis Fed, reproduced in Appendix Figure (ref) for reference.\footnote{Chauvet, Marcelle, and Jeremy M.\ Piger, Smoothed U.S. Recession Probabilities \href{https://fred.stlouisfed.org/series/RECPROUSM156N}{[RECPROUSM156N]}, retrieved from Federal Reserve Economic Data (FRED), Federal Reserve Bank of St.\ Louis, December 31, 2025.} Next, we fit the same model to the components of GDP to see whether there are gains from using disaggregate data. We consider a system of 30 equations, each for a component of GDP (excluding government spending because it is not procyclical):
where $S_t=1$ and $S_t=2$ are as before; $e_{i,t} \sim \text{i.i.d. } N(0,\sigma_{i,1}^2)$ if $S_t=1$ and $e_{i,t} \sim \text{i.i.d. } N(0,\sigma_{i,2}^2)$ if $S_t=2$. The errors $e_{i,t}$ are specified as cross-sectionally independent. Figure (ref)(b) reports the smoothed recession probabilities. All NBER recessions, including the one in the early 2000s, are identified, and the turning points are sharper than in the aggregate case. Since the two models share the same functional form and are estimated using the same sampling algorithm and prior, the difference in results is likely due to disaggregation. In Section (ref) we show that still sharper results can be obtained by incorporating additional disaggregated time series in the estimation.
The above comparison shows that disaggregated data can improve regime detection even without changing the model in an essential way. The reason is that the components of GDP differ in their signal-to-noise ratios with respect to regime switching. Aggregation ignores this heterogeneity and fails to downweight noisier series, while disaggregated data allow the likelihood to adapt accordingly. In Section (ref), using a simple two-variable example, we show analytically that aggregation preserves all information about regime switching if and only if the signal-to-noise ratios are equal across the two series. When they differ, as in the present application, the information loss can be substantial. We also quantify the gains through simulations.
Our proposed models are generalizations of (ref) that allow for autoregressive dynamics and grouped heterogeneity. In a typical application, the model contains hundreds of parameters, and standard maximum likelihood often fails to converge. To address this, we develop a Gibbs sampling algorithm for estimation, inference, and forecasting, which incorporates a new multi-step sampler for latent state variables. We show that the algorithm can handle large macroeconomic datasets (several hundred time series) in reasonable computation time.
We apply the models to a dataset of 94 disaggregated series, covering components of real GDP, industrial production, capacity utilization, employment, and hours worked. The models detect all NBER-defined recessions over 1972-2024. The smoothed recession probabilities are close to being binary, taking values near zero or one with few intermediate values. This contrasts with typical estimates in the literature, such as those in Figure (ref), which tend to linger at intermediate values for extended periods. The estimated recession starting dates either coincide exactly with the NBER dates, or in a few cases, lag them by one to three quarters. Meanwhile, the estimated recession ending dates closely track the starts of NBER expansions. This asymmetry may reflect the nature of business cycles: downturns build gradually across sectors, while recoveries tend to generate a more synchronized signal across disaggregated series. These findings suggest that methods designed around a wide set of disaggregated variables can be useful for business cycle research, particularly for accurate regime identification.
The issue of cross-sectional aggregation has been studied mainly for linear models such as ARMA models and for forecasting; see, e.g., lutkepohl1987, Granger1988, Luetkepohl2009. To our knowledge, this is the first formal analysis of aggregation in regime-switching models. Here, the object of interest is a common regime component shared across series rather than a linear forecast. Our results suggest that disaggregation gives clear gains in this setting. We conjecture that the same insight extends to other settings where a shared latent component is of interest, such as nonlinear factor models or models with common stochastic volatility. In those cases, working with less aggregated series may produce similar gains.
The paper is organized as follows. Section 2 introduces the models and examines gains from disaggregation. Section 3 presents a Gibbs sampler for estimation, filtering, and forecasting. Section 4 studies the methods through simulations. Section 5 presents an empirical application, and Section 6 concludes. Additional figures and details are available in the online appendix.
We consider $N$ time series partitioned into $J$ groups. Denote the series in group $j$ by \[ y_{j,t} = (y_{j,1,t}, y_{j,2,t}, \ldots, y_{j,N_j,t})^{\prime}, \] where $N_j$ is the number of series in group $j$, and $N = \sum_{j=1}^{J} N_j$. If each group contains only a single series, then $J = N$. The full data vector at time $t$ (for $t = 1,\ldots,T$) is
Let \( Y_{j,t} = \sum_{i=1}^{N_j} y_{j,i,t} \) be the sum of the series in group $j$. The vector of aggregated series at $t$ is
Our focus is on the disaggregate series in (ref); the aggregate series in (ref) are included for comparison. Throughout, uppercase letters denote aggregate variables.
We propose two models. Each consists of a common regime-switching component and an autoregressive component that captures additional dynamics. Some coefficients are equal within each group to capture grouped heterogeneity. Both the mean and variance can switch.
The first model is a generalization of hamilton1989 to the disaggregated setting:
where $(j,i,t)$ indexes the $i$-th unit in the $j$-th group at time $t$; $k$ is the lag order; and $(S_{t-k},\ldots,S_t)$ is a vector of latent variables determining the state of the system at $t$. In general, $S_t$ may take values in a finite state space $S_t\in\mathcal{S}=\{1,2,\ldots,S\}$, with known $S$; in the leading case, $S_t$ is binary, representing recession and expansion. We assume $S_t$ follows a Markov process with transition probabilities $p_{ss'}$ for any $s,s'\in\mathcal{S}$. The errors are assumed to be normally distributed: \[ e_{j,i,t}\sim N(0,\sigma_{j,i,S_t}^2). \] The unknown parameters consists of the group-level autoregressive coefficients $\phi_{j,m}$, the regime-dependent means $\mu_{j,i,S_{t}}$, and the regime-dependent variances $\sigma_{j,i,S_t}^2$.
A known feature of regime-switching models is label switching: without an ordering restriction, the regimes are not identified. We impose the following ordering restriction on the means: \[ \mu_{j,i,1} \le \mu_{j,i,2} \le \cdots \le \mu_{j,i,S}, \qquad \text{for any } j \text{ and } i . \] Because $S_t$ is a scalar process, restricting the mean of a single series is in fact sufficient for identification; e.g., one could require the first series in the first group to satisfy $\mu_{1,1,1} \le \mu_{1,1,2} \le \cdots \le \mu_{1,1,S}$. In our applications, however, we find that imposing the restriction for all series improves computational performance. For variables such as the unemployment rate, whose values are higher during recessions, we multiply by $-1$ so that they satisfy these inequalities.
In the two-regime case, it is convenient to write \[ \mu_{j,i,S_t} = \delta_{j,i,1}\mathbf{1}_{\{S_t=1\}}+\mu_{j,i,2}, \qquad \delta_{j,i,1}\le 0 \ \text{for all } i,j . \] Under this parameterization, $\mu_{j,i,2}$ is unrestricted, and each ordering restriction involves a single coefficient rather than two. The reduction is particularly useful in high-dimensional settings, as it removes $N$ parameters from the set of coefficients subject to restrictions during estimation.
The aggregate model implied by Model A is
where $e_{j,t}\sim N(0,\sigma_{j,S_t}^2$), \( \mu_{j,S_t} =\sum_{i=1}^{N_j}\mu_{j,i,S_t}$, and $\sigma_{j,S_t}^2 =\sum_{i=1}^{N_j}\sigma_{j,i,S_t}^2. \)
The second model generalizes the univariate models in hansen1992, Cho2007, and qu2021:
The error terms $e_{j,i,t}$, the latent state $S_{t}$, and the parameters ($\phi_{j,m}$, $\mu_{j,i,S_{t}}$, $\sigma_{j,i,S_t}^2$) are the same as in Model A. The corresponding aggregate model is
where $\mu_{j,S_t}$, $\phi_{j,m}$, and $e_{j,t}$ are defined as before.
Models A and B share similar second-order dynamics but are not observationally equivalent. For example, when $k=1$, both models have an ARMA(1,1) representation via the AR(1) representation of the latent Markov chain. They differ in how the series reacts to regime changes. Model A implies an immediate jump toward the new regime, while Model B implies a more gradual transition as the lagged level $y_{t-1}$ carries information from the previous regime. As a result, they imply different higher-order moments, as in the univariate setting.
Models A and B are designed for situations where average growth rates determine regimes. In other settings, regimes may depend on other parameters, such as a set of regression coefficients. Model B can be generalized to such applications. For example, one may consider \[ y_{j,i,t} = x_{j,i,t}'\beta_{j,i} + z_{j,i,t}'\delta_{j,i,S_t} + e_{j,i,t}, \] where $x_{j,i,t}$ and $z_{j,i,t}$ are random vectors. The coefficients $\beta_{j,i}$ remain constant and may exhibit a grouped structure, while $\delta_{j,i,S_t}$ are regime dependent. The error term $e_{j,i,t}$ and the latent states are defined as in Model B. The proposed algorithm can be adapted to estimate this model.
The following assumptions are familiar in the Markov-switching literature and ensure that the models are identified and that estimation is feasible. Their main role is to help assess whether the models are suitable for the intended applications.
This assumption implies a common cycle affects all data units with the same timing. It is nonrestrictive when all series are coincident indicators of the business cycle, such as changes in non-farm employment, industrial production, capacity utilization, real personal income, and real manufacturing and trade sales. It is violated when coincident indicators are combined with leading indicators, such as stock market indices or the term structure of interest rates. It implies that the disaggregate variables should be chosen so that they respond to $S_t$ with approximately the same timing. Economic theory and empirical evidence can often guide this choice.
This assumption implies that every state can be reached from every other state with positive probability. This property ensures the existence of a unique stationary distribution and is necessary for identification and inference based on the likelihood function.
This assumption addresses the label-switching property of regime-switching models. Requiring at least one strict inequality ensures that the regimes are distinguishable in the data without restricting heterogeneity across series and groups.
The normality assumption reflects the parametric nature of the regime-switching model and is standard in the literature. Setting the off-diagonal elements of the covariance matrix to zero implies that the joint likelihood factorizes as the product of the marginal likelihoods of the individual series. This is analogous to a composite likelihood approach: cross-sectional dependence in the errors is intentionally not modeled, which may sacrifice estimation efficiency but reduces computational costs and the risk of misspecification.\footnote{See varin2011overview for a review of composite likelihoods.}
This assumption places mild restrictions on the parameters. It rules out degenerate cases by requiring finite means and positive variances. The stability condition on the lag polynomials implies stationary autoregressive dynamics within each regime. These conditions are needed for the numerical stability of the proposed estimation algorithm.
We consider a simple setting, with two disaggregated time series and known parameters, to examine the effects of aggregation on the identification of recessions. This can be viewed as an asymptotic exercise in which all parameters are assumed to have been consistently estimated. In finite samples, the disaggregated model involves more parameters and greater estimation uncertainty. Simulations later in the paper examine this issue.
Regime switching in mean and in variance each carries useful information. We consider the two cases separately. Lagged dependent variables (as in Model B) do not affect regime updating, since the autoregressive coefficients do not switch; they are omitted.
\paragraph{Case 1: Switching in the mean only.}
We consider two time series given by
The aggregated series is \( Y_{t}=\mu_{1}\mathbf{1}_{\{S_{t}=1\}}+\mu_{2}\mathbf{1}_{\{S_{t}=2\}}+e _{t}, \ \ e _{t}\sim N(0,\sigma ^{2}), \) where $Y_{t}=y_{1,t}+y_{2,t}$, $\mu_{1}=\mu _{1,1}+\mu _{2,1}$, $\mu_{2}=\mu _{1,2}+\mu _{2,2}$, and $\sigma ^{2}=\sigma _{1}^{2}+\sigma _{2}^{2}$. Assume $\mu _{1,1}\leq\mu _{1,2}$ and $\mu _{2,1}\leq\mu _{2,2}$.
For the disaggregated model, let $\hat{S} _{t|t-1}^{1}\equiv E(\mathbf{1}_{\{S_{t}=1 \}}|y_{t-1},y_{t-2},\ldots)$ be the predicted probability of recession, with conditional mean squared error $\hat{P} _{t|t-1}^{1}\equiv E\{(\mathbf{1}_{\{S_{t}=1\}}-\hat{S} _{t|t-1}^{1})^{2}|y_{t-1},y_{t-2},\ldots\}=\hat{S}_{t|t-1}^{1}(1-\hat{S} _{t|t-1}^{1})$. Similarly, for the aggregated model, define $\Breve{S}_{t|t-1}^{1}\equiv E(\mathbf{1}_{\{S_{t}=1 \}}|Y_{t-1},Y_{t-2},\ldots)$ and $\Breve{P} _{t|t-1}^{1}\equiv E\{(\mathbf{1}_{\{S_{t}=1\}}-\Breve{S} _{t|t-1}^{1})^{2}|Y_{t-1},Y_{t-2},\ldots\}$. We use $\hat{S}^1_{t\mid t}$ and $\Breve{S}^1_{t\mid t}$ for filtered probabilities, and $\hat{S}^1_{t\mid T}$ and $\Breve{S}^1_{t\mid T}$ for smoothed probabilities. The following result characterizes when the aggregated and disaggregated models are equivalent, and when they are not.
The proposition implies that unless the signal-to-noise ratios are equal as in (ref), the two filters differ with probability one. This is because filtering is driven by the log-likelihood ratio comparing the two regimes. In the disaggregated model, the ratio is \[ \ell_t^{(d)} = \text{const} + \frac{\delta_1}{\sigma_1^2}\,(y_{1,t}-\mu_{1,1}) + \frac{\delta_2}{\sigma_2^2}\,(y_{2,t}-\mu_{2,1}), \qquad \delta_i=\mu_{i,1}-\mu_{i,2}, \] where $y_{1,t}$ and $y_{2,t}$ provide regime information through two separate weights $\delta_1/\sigma_1^2$ and $\delta_2/\sigma_2^2$. In the aggregate model, the ratio is \[ \ell_t^{(a)} = \text{const} + \frac{\delta_1+\delta_2}{\sigma_1^2+\sigma_2^2} \bigl[(y_{1,t}+y_{2,t})-(\mu_{1,1}+\mu_{2,1})\bigr], \] and regime evidence enters only through the sum, i.e., with equal weights on $y_{1,t}$ and $y_{2,t}$. If condition (ref) holds, the two disaggregated weights are equal, and the disaggregated likelihood ratio reduces to a function of $Y_t = y_{1,t}+y_{2,t}$ alone. Otherwise, there exist movements in $(y_{1,t},y_{2,t})$ that change the disaggregated likelihood ratio but leave the aggregate unchanged, and aggregation discards information useful for filtering. The Appendix contains the formal proof.
The next result shows that the disaggregated filter is weakly more precise than the aggregate filter, and strictly more precise when (ref) does not hold. The precision is measured in terms of the mean squared error: for the disaggregated model $E\{(\mathbf{1}_{\{S_{t}=1\}}-\hat{S} _{t|t}^{1})^{2}\}=E\{\hat{S}_{t|t}^{1}(1-\hat{S} _{t|t}^{1})\}$, and for the aggregrated model $E\{(\mathbf{1}_{\{S_{t}=1\}}-\Breve{S} _{t|t}^{1})^{2}\}=E\{\Breve{S}_{t|t}^{1}(1-\Breve{S} _{t|t}^{1})\}$.
To supplement the theoretical result, we consider a simulation experiment on how relative filtering precision depends on the signal-to-noise ratio in ((ref)). We report the mean squared error (MSE) and mean absolute error (MAE) of the filtered probabilities for $S_t=1$. For the disaggregated model, these are given by $ E\{(\hat{S}_{t\mid t}^1-\mathbf{1}_{\{S_t=1\}})^2\} $ and $E\{\lvert \hat{S}_{t\mid t}^1-\mathbf{1}_{\{S_t=1\}}\rvert\}, $ and for the aggregated model by $ E\!\{(\Breve{S}_{t\mid t}^1-\mathbf{1}_{\{S_t=1\}})^2\} $ and $E\{\lvert \Breve{S}_{t\mid t}^1-\mathbf{1}_{\{S_t=1\}}\rvert\}$. All quantities are averaged over 1{,}000 simulation replications. To cover a range of cases, we set $\sigma_{1}^{2}=1$, $\mu_{1,2}-\mu_{1,1}=\mu_{2,2}-\mu_{2,1}=2$, and let $\sigma_{2}^{2}$ vary from $0.1$ to $10$. Two extreme cases, $\sigma_{2}^{2}=0.01$ and $\sigma_{2}^{2}=100$, are also considered. The transition probabilities are set to $p_{11}=0.75$ and $p_{22}=0.95$, consistent with typical estimates for U.S. business cycles.
The results are reported in Table (ref). When $\sigma_2^2=1$, the two models perform identically, as condition (ref) predicts. Their differences increase monotonically as $\sigma_2^2$ decreases from 1. When $\sigma_2^2=0.5$, the MSE under disaggregation is about 18% lower, with a similar reduction in MAE. When $\sigma_2^2=0.2$, the improvement increases to 78% under both loss measures; when $\sigma_2^2=0.1$, it reaches 97%. When $\sigma_2^2>1$, the disaggregated model also does better, where the MSE reduction increases from 11% at $\sigma_2^2=2$ to 46% at $\sigma_2^2=10$, and reaches 60% when $\sigma_2^2=100$. Even with $\sigma_2^2=10$, the standard deviation of the second series is only $\sqrt{10}$ times that of the first, a difference that remains empirically relevant. In the extreme case $\sigma_2^2=100$, the second series contains essentially no regime information. The disaggregated filter adapts to this by downweighting the noisy series, and its MSE is only twice the $\sigma_2^2=1$ benchmark. The aggregated filter becomes completely uninformative. The results show that aggregation can lead to substantial loss when signal-to-noise ratios differ across series.
\paragraph{Case 2: Switching in the variance only.}
We consider two time series given by
The aggregate model is $ Y_{t}=\mu+e _{t}$, with $e _{t}\sim N(0,\sigma _{1}^{2}\mathbf{1} _{\{S_{t}=1\}}+\sigma _{2}^{2}\mathbf{1}_{\{S_{t}=2\}}), $ where $Y_{t}=y_{1,t}+y_{2,t}$, $\mu=\mu _{1}+\mu _{2}$, $\sigma _{1}^{2}=\sigma _{1,1}^{2}+\sigma _{2,1}^{2}$, and $\sigma _{2}^{2}=\sigma _{1,2}^{2}+\sigma _{2,2}^{2}$. Assume $\sigma_{1,1}\geq \sigma_{1,2}$ and $\sigma_{2,1}\geq \sigma_{2,2}$.
The proposition implies that the filters are different whenever variances switch. In the disaggregated model, the log-likelihood ratio comparing the two regimes is \[ \ell_t^{(d)} = \text{const} + \tilde y_{1,t}^{2} \left( \frac{1}{2\sigma_{1,1}^{2}}-\frac{1}{2\sigma_{1,2}^{2}} \right) + \tilde y_{2,t}^{2} \left( \frac{1}{2\sigma_{2,1}^{2}}-\frac{1}{2\sigma_{2,2}^{2}} \right), \] where $\tilde y_{i,t}=y_{i,t}-\mu_i$, and each series provides regime information through its own squared deviation $\tilde y_{i,t}^2$. In the aggregated model, the log-likelihood ratio is \[ \ell_t^{(a)} = \text{const} + (\tilde y_{1,t}+\tilde y_{2,t})^{2} \left( \frac{1}{2(\sigma_{1,1}^{2}+\sigma_{2,1}^{2})} - \frac{1}{2(\sigma_{1,2}^{2}+\sigma_{2,2}^{2})} \right), \] which depends on the data only through the squared deviation of the aggregate: $(Y_t-\mu)^2=(\tilde y_{1,t}+\tilde y_{2,t})^2$. For the two filters to coincide, the two terms in $\ell_t^{(d)}$ must combine into a multiple of $(\tilde y_{1,t}+\tilde y_{2,t})^2$. However, this cannot hold except in the trivial case of no-switching: a weighted sum $a\,\tilde y_{1,t}^2 + b\,\tilde y_{2,t}^2$ matches $c(\tilde y_{1,t}+\tilde y_{2,t})^2$ for all $(\tilde y_{1,t}, \tilde y_{2,t})$ only when $c = 0$ (since the latter contains an additional cross term $2c\,\tilde y_{1,t}\tilde y_{2,t}$ which must be set to zero), and this forces $a = b=0$. Even when the aggregate variance is constant across regimes, i.e., $\ell_t^{(a)}$ carries no regime information, the disaggregated likelihood continues to distinguish states through the relative scaling of individual squared deviations. Equivalence happens only in the degenerate case in which variances do not switch at all. The next result confirms that the disaggregated filter is strictly more accurate whenever at least one variance is switching.
To supplement the theoretical result, we run a simulation experiment similar to the switching-mean case. We set $\sigma_{1,1}^{2}=\sigma_{1,2}^{2}=\sigma_{2,1}^{2}=1$ and vary $\sigma_{2,2}^{2}$ from 0.1 to 10. Two extreme cases $\sigma_{2,2}^{2}=0.01$ and $\sigma_{2,2}^{2}=100$ are also included. Transition probabilities are the same as in the switching-mean case. Table (ref) reports the MSE and MAE of filtered regime probabilities.
The benchmark case $\sigma_{2,2}^2=1$ corresponds to no variance switching; the two filters perform identically in terms of MSE and MAE, as predicted by Proposition 2. When $\sigma_{2,2}^2<1$ (Panel A), the overall variance of the second series decreases. The difference in MSE and MAE between the two filters grows monotonically as $\sigma_{2,2}^2$ moves away from one, reaching about 80% improvement at $\sigma_{2,2}^2=0.01$. A similar pattern is seen under $\sigma_{2,2}^2>1$ (Panel B), although the gains are more modest. In this case, the overall variance of the second series increases, which spreads out probability mass and flattens the likelihood, making regime switching harder to detect. In summary, the results confirm that aggregation loses no information only in the trivial case of no switching; otherwise the disaggregated filter is strictly more precise. In addition, switches into a lower-variance regime are easier to detect, and aggregation costs more in that case.
Let $\theta$ denote the vector of model parameters, including the means $\mu_{j,i,s}$, variances $\sigma^{2}_{j,i,s}$, autoregressive coefficients $\phi_{j,m}$, and transition probabilities $p_{ss'}$, where $j$ indexes groups, $i$ indexes series within group $j$, $m$ denotes the lag order, and $s$ indexes regimes. The dimension of $\theta$ is \( M = 2S\sum_{j=1}^{J} N_j + kJ + S(S-1), \) which exceeds 200 when there are 50 series. In such settings, likelihood-based optimization often fails to converge. Instead, we use a Bayesian approach that samples $\theta$ and the latent state sequence $\{S_t\}_{t=1}^T$ jointly. The exposition focuses on the disaggregate specification in (ref); differences that arise under (ref) will be noted where relevant.
Our sampler adapts the sampling algorithms in albert1993, kim1998, and kim1999 to a high-dimensional setting. The goal is to sample from the following joint posterior distribution conditional on $y_{1:T}=(y_{1}^{\prime},y_{2}^{\prime},\ldots,y_{T}^{\prime})^{\prime}$:
where $S_{1:T}$ is the latent state sequence over the sample, $l$ is the forecast horizon, and $y_{T+1},\ldots,y_{T+l}$ and $S_{T+1},\ldots,S_{T+l}$ are future observables and latent states, respectively. As often done in the literature, we factor this joint distribution as a conditional and a marginal distribution:
Once the within-sample latent states and parameter values have been sampled, future states and observables can be generated conditional on these draws. For the distribution $f\!\left(S_{1:T},\theta \mid y_{1:T}\right)$, we consider its full set of conditionals and sample from them iteratively:
Here $\mu_{j,i,s}$, $\sigma_{j,i,s}^2$, and $\phi_{j,m}$ denote the full parameter collections $\{\mu_{j,i,s}\}_{j,i,s}$, $\{\sigma_{j,i,s}^2\}_{j,i,s}$, and $\{\phi_{j,m}\}_{j,m}$.
Next we provide details for each step, with an emphasis on the aspects related to the high-dimensional setting. We assume there are two states, $S_t \in \{1,2\}$. Note that for a given \(t\), the log likelihood conditional on \(S_{t-k:t}\) is
where $ \varepsilon_{j,i,t} = y_{j,i,t} - \mu_{j,i,S_t} - \sum_{m=1}^{k} \phi_{j,m} \bigl(y_{j,i,t-m}-\mu_{j,i,S_{t-m}}\bigr)$ for $t=k+1,\ldots,T . $ Summing over \(t=k+1,\ldots,T\), we get the log likelihood for the full sample:
In each iteration, we first sample $\mu_{j,i,s}$, followed by $\sigma_{j,i,s}^{2}$, and then $\phi_{j,m}$.
The regime-dependent means are sampled from $ f\!( \mu_{j,i,s} \mid y_{1:T}, S_{1:T}, p_{ss'},\sigma_{j,i,s}^{2}, \phi_{j,m})$ at once. Because the sampling is conditional on \(\{\phi_{j,m}\}_{m=1}^k\) and \(S_{1:T}\), Model A can be rewritten as \[ y_{j,i,t} - \sum_{m=1}^{k}\phi_{j,m}y_{j,i,t-m} = \delta_{j,i,1} \bigl( \mathbf{1}_{\{S_{t}=1\}} - \sum_{m=1}^{k}\phi_{j,m}\mathbf{1}_{\{S_{t-m}=1\}} \bigr) + \mu_{j,i,2} \bigl( 1-\sum_{m=1}^{k}\phi_{j,m} \bigr) + e_{j,i,t}. \] For each \(t\), this represents a system of \(N\) seemingly unrelated regression equations, with outcome \(y_{j,i,t}-\sum_{m=1}^{k}\phi_{j,m}y_{j,i,t-m}\), regressors \( \mathbf{1}_{\{S_{t}=1\}} - \sum_{m=1}^{k}\phi_{j,m}\mathbf{1}_{\{S_{t-m}=1\}} \quad\text{and}\quad 1-\sum_{m=1}^{k}\phi_{j,m}, \) coefficients \(\delta_{j,i,1} = \mu_{j,i,1} - \mu_{j,i,2}\le 0\) and \(\mu_{j,i,2}\), and normal error terms with a diagonal covariance matrix. Let $\delta_{1}$ be a vector containing all $\delta_{j,i,1}$, and $\mu_{2}$ all \(\mu_{j,i,2}\). Under the truncated normal priors that we use, the posterior distribution of \(\delta_{1}\) is a truncated multivariate normal, and that of \(\mu_{2}\) is a standard multivariate normal (see Zellner1971). Our sampling relies on vectorized operations and exploits independence across equations, so that the required matrix inversion has dimension at most \(2N\), regardless of the number of observations in the system. Sampling for the truncated normal uses the minimax tilting method of Botev2017, specifically designed for high-dimensional data. Details on these aspects can be found in the Appendix. Monte Carlo experiments show that drawing \(\delta_{1}\) and \(\mu_{2}\) jointly is feasible for up to a few hundred time series. When the cross-sectional dimension is even larger, the elements of these vectors can be sampled in blocks, provided that the priors are independent across parameters.
The variances are sampled from $f\!( \sigma_{j,i,s}^{2} \mid y_{1:T}, S_{1:T}, p_{ss'},\mu_{j,i,s}, \phi_{j,m}).$ The procedure uses independence across time series and samples the variances one at a time. The explicit sampling steps are identical to those in kim1999. The variances are expressed as
where \(h_{j,i,2}=\sigma_{j,i,1}^{2}/\sigma_{j,i,2}^{2}-1\). With independent inverse-Gamma priors on \(\sigma_{j,i,2}^{2}\) and \(1+h_{j,i,2}\), the conditional posterior distributions are also inverse-Gamma and can be sampled directly. Each step is a one-dimensional problem and does not require any matrix inversion. Implementation details are included in the Appendix.
The autoregressive coefficients are sampled from \( f\!( \phi_{j,m} \mid y_{1:T}, S_{1:T}, p_{ss'}, \mu_{j,i,s}, \sigma_{j,i,s}^{2}). \) These are drawn separately from $\mu_{j,i,s}$ to avoid nonlinear interactions between them. Following the same approach as for $\mu_{j,i,s}$, we write Model A as \[ y_{j,i,t}-\mu_{j,i,S_t} = \sum_{m=1}^{k}\phi_{j,m}\bigl(y_{j,i,t-m}-\mu_{j,i,S_{t-m}}\bigr) + e_{j,i,t}, \qquad e_{j,i,t}\mid S_t \sim N(0,\sigma_{j,i,S_t}^2). \] Stacking across groups gives a system with dependent variables $y_{j,i,t}-\mu_{j,i,S_t}$, regressors $y_{j,i,t-m}-\mu_{j,i,S_{t-m}}$, common coefficients $\phi_{j,m}$ within group, and diagonal residual covariance. Under a normal prior on $\phi_{j,m}$, the conditional posterior distribution is normal and can be sampled directly. As in the case of sampling $\mu$, we adopt a vectorized representation so that the required matrix inversion has dimension at most $Jk$, regardless of the cross-sectional dimension. Implementation details are provided in the Appendix.
The latent states are sampled from $ f\!\left( S_{1:T} \mid y_{1:T}, \theta\right)$. We build on the multi-move Gibbs sampler of kim1998, which generates the entire sequence of latent states at once. It does so by factoring the joint posterior distribution into conditionals and sampling backward in time: \[ f\!\left(S_{1:T}\mid y_{1:T},\theta\right) = f\!\left(S_T\mid y_{1:T},\theta\right) \prod_{t=1}^{T-1} f\!\left(S_t \mid S_{t+1:T},\, y_{1:T},\,\theta\right). \] Let \(k\ge 0\) denote the number of state lags entering the observation density. In the Appendix, we show that the following relationship between smoothing and filtering holds for Model A:
for any \(k>0\) and $p\!\left(S_t \mid S_{t+1:T},\, y_{1:T},\, \theta\right) \;\propto\; p\!\left(S_{t:t+1}\mid y_{1:t+1},\, \theta\right)$ for $k=0.$ This delivers an exact sampler for the smoothed distribution (left-hand side) as a by-product of the state filtering step in the filter of hamilton1989 (right-hand side). The sampler of kim1998 is exact only if \(k=0\); when \(k>0\), it is approximate because the dependence on some states is omitted.
In implementation, the right hand side of ((ref)) is computed forward recursively as follows. Given \(p(S_{t-k-1:t-1}\mid y_{1:t-1},\theta)\), the Markov property of \(\{S_t\}\) implies the following prediction:
Then, given $y_t$, the state probability can be updated as:
This update is carried out using log weights because the observation density involves products over many series, which can become numerically small. Specifically, for each configuration \((S_{t-k},\ldots,S_t)\), we compute the unnormalized log weight as the sum of the log observation density in (ref) and the log predictive probability from (ref). To improve numerical stability, the maximum log weight across all configurations is subtracted before exponentiation, after which the resulting weights are normalized to sum to one. Given the filtered joint distributions \(\{p(S_{t-k:t}\mid y_{1:t},\theta)\}\), sampling proceeds backward starting at \(T\). At time \(t\), one draws from \( p\!\left(S_t,\, S_{t+1:t+k}\mid y_{1:t+k},\, \theta\right), \) which is obtained directly from the filtered probabilities at time \(t+k\). When \(t+k>T\), the conditioning set \(S_{t+1:t+k}\) is truncated accordingly. These operations do not involve any matrix operations, and the sampling involves drawing Bernoulli random variables. The procedure remains feasible for large cross-sections. To our knowledge, this exact state sampler is new.
The elements of the transition matrix are sampled from $ f\!\left( p_{ss'} \mid y_{1:T}, S_{1:T}, \mu_{j,i,s}, \sigma_{j,i,s}^{2}, \phi_{j,m} \right)= f\!\left( p_{ss'} \mid S_{1:T} \right)$, as in kim1998. Since \(S_t\) is a scalar process, the sampling is independent of the model's dimension. We use independent Beta priors for \(p_{11}\) and \(p_{22}\), which implies the posterior distributions are also Beta. The Appendix provides the exact expressions.
The values of $y_{T+1:T+l}$ for any $l\geq1$ are sampled from $ f(y_{T+1:T+l}, S_{T+1:T+l}\mid y_{1:T}, S_{1:T}, \theta). $ Since $ f(y_{T+1:T+l}, S_{T+1:T+l}\mid y_{1:T}, S_{1:T}, \theta) = \prod_{j=1}^{l} f(y_{T+j}, S_{T+j}\mid y_{1:T+j-1}, S_{1:T+j-1}, \theta), $ the joint distribution can be sampled sequentially from the $l$ conditional distributions on the right-hand side. For each $j=1,\ldots,l$, we draw $S_{T+j}$ from the transition distribution $ f(S_{T+j}\mid S_{T+j-1}, \theta), $ and then draw $y_{T+j}$ from $ f(y_{T+j}\mid y_{1:T+j-1}, S_{1:T+j}, \theta). $ Repeating this procedure for $j=1,\ldots,l$ gives a draw from $ f(y_{T+1:T+l}, S_{T+1:T+l}\mid y_{1:T}, S_{1:T}, \theta). $ Applying this step to each posterior draw of $(\theta, S_{1:T})$ produces predictive simulations of future observations.
This part is essentially the same as in low-dimensional regime-switching models. It is tractable in the high-dimensional case because the state $S_{T+j}$ is common across all series. Only one draw is needed for the state, and the observations are then drawn independently across series.
The following prior distributions are used for the simulation experiments and the empirical application. The priors are assumed to be independent unless stated otherwise.
\paragraph{Regime-dependent means.} For each series, the regime-dependent means are assigned priors $ \delta_{j,i,1}=\mu_{j,i,1}-\mu_{j,i,2} \sim N(-0.5,\,50^2)$ and $\mu_{j,i,2} \sim N(0,\,50^2). $ The prior variances are set to be diffuse, reflecting the growth rates can differ substantially across disaggregated series.
\paragraph{Innovation variances.}
When the variances are restricted to be constant across regimes, we set $\sigma_{j,i}^2 \sim IG(0,0)$, which corresponds to the improper prior $p(\sigma_{j,i}^2)\propto (\sigma_{j,i}^2)^{-1}$. When the variances are allowed to switch across regimes, we still assign $\sigma_{j,i,2}^2 \sim IG(0,0)$, and as in kim1999, we write $\sigma_{j,i,1}^2 = (1+h_{j,i,2})\sigma_{j,i,2}^2$ and adopt a conjugate inverse-gamma prior for $1+h_{j,i,2}$. We set $1+h_{j,i,2} \sim IG\!\left(\tfrac{T_{1}}{2},\,\tfrac{T_{1}+2}{2}\right)$, where $T_{1}$ denotes the number of observations assigned to regime 1 under the current draw of $S_{1:T}$. This represents a state-dependent prior rather than a fixed prior. Under this specification, the mode of $1+h_{j,i,2}$ is equal to 1 (the prior for the recession variance peaks at the expansion variance), and the amount of shrinkage increases with $T_{1}$. The shrinkage is helpful because the recession regime typically contains few observations. The prior also allows for right-skewed dispersion in the variance ratio.
\paragraph{AR coefficients.} For $k>0$, the AR coefficients (by group) are assigned Gaussian prior \( N(0,\,0.5^2), \) with the stationarity restriction enforced by rejection sampling, i.e., draws are accepted only if the roots of the AR polynomial lie outside the unit circle.
\paragraph{Transition probabilities.} For the Markov chain, the diagonal transition probabilities $p_{11}$ (recession) and $p_{22}$ (expansion) are assigned priors $ p_{11}\sim \mathrm{Beta}(2,2)$ and $p_{22}\sim \mathrm{Beta}(30,2), $ with means $0.50$ and $0.94$, and standard deviations $0.22$ and $0.04$. These values are informed by previous estimates of regime-switching models for U.S. business cycles. The tight prior on $p_{22}$ reflects the consistent empirical evidence and reinforces the separation between regimes.
This section examines finite-sample effects of disaggregation on state recovery. We begin with a two-variable setting in which the signal-to-noise ratio varies and compare aggregated and disaggregated filters. This analysis complements the asymptotic characterization in Section 2.2 by incorporating parameter estimation uncertainty. We then consider a setting where the signal-to-noise ratio is fixed while the cross-sectional dimension increases, and examine how filtering and smoothing change as additional series are included. The sample size set to 190, matching the empirical sample size used to estimate model ((ref)). Parameters are estimated via MCMC using the priors in Section 3.5. All results are averages over 100 independent replications.
We use the DGP in ((ref)) with the same parameter values. The estimated model has the same specification as the DGP, with no lags, and assumes constant variances across regimes. We report filtered probabilities, as those in Table (ref); results based on smoothed probabilities are similar.
Table (ref) shows that, relative to Table (ref), where parameters are known, the MSE and MAE are generally larger for both filters, reflecting the uncertainty from parameter estimation. At the benchmark $\sigma_2^2 = 1$, the differences in MSE and MAE between the aggregated and disaggregated filters are negligible, in line with the theory. Gains from disaggregation emerge as $\sigma_2^2$ moves away from 1 in either direction. For $\sigma_2^2 < 1$, the improvement reaches about $50\%$ at $\sigma_2^2 = 0.3$ and close to $100\%$ for $\sigma_2^2 \leq 0.10$. For $\sigma_2^2 > 1$, gains grow more slowly but are still substantial -- about $30\%$ by $\sigma_2^2 = 4$ and above $45\%$ at $\sigma_2^2 = 10$. The monotonic pattern seen in Table (ref) is preserved under parameter estimation uncertainty.
Therefore, although parameter uncertainty increases the MSE and MAE for both filters, it does not alter the main finding: aggregation attenuate regime signals, and the gains from disaggregation can be large when heterogeneity is substantial.
We calibrate the DGP to empirical estimates from model ((ref)). These estimates provide 30 sets of regime-specific means and variances, along with a common transition matrix. For simulation, we first generate a common latent state vector of length 190 using this transition matrix. We then draw one of the 30 parameter sets at random, scale the regime difference $\delta_{i,1} = \mu_{i,1} - \mu_{i,2}$ by one-half, and use it together with the latent state vector to generate a time series of length 190 from model ((ref)). This procedure is repeated to obtain $N$ series, with $N \in \{5,10,20,40,60,80,100,200,400\}$. The estimated model matches the specification of ((ref)), with no lags and variances allowed to switch across regimes.
Here the regime difference $\delta_{i,1}$ is multiplied by 0.5 to keep regime classification nontrivial even in large cross sections. Otherwise, the advantage of disaggregation remains, but the regimes are estimated with little uncertainty in both the aggregated and disaggregated models once $N$ exceeds 60, making this experiment uninformative about relative performance at higher dimensions.
Table (ref) reports the MSE and MAE of inference on $S_t$, with panel A for filtering and panel B for smoothing. For filtering, the two models perform comparably at $N = 5$. With $N = 10$ series, disaggregation reduces the filtered MSE from 0.1300 to 0.0872, an improvement of about 32%. At $N = 20$, the filtered MSE is 0.1171 for the aggregated model and 0.0416 for the disaggregated model, a reduction of about 65%. At $N = 100$, the filtered MSE is 0.0346 for the aggregated model and 0.0023 for the disaggregated model, a reduction of about 93%. Smoothing shows similar percentage reductions. As $N$ further increases from 100, MSEs for both approaches become small. This reflects the design: since the signal-to-noise ratio of each series is held fixed, the total information increases with $N$. These results suggest that disaggregation improves state inference across a wide range of cross-sectional sample sizes in this empirically calibrated setting.
The NBER identifies US expansions and recessions based on a broad set of economic indicators, and its dating decisions are not tied to any single econometric model. Markov regime-switching model is arguably the most important econometric model for identifying business cycles. It provides a probabilistic statement about the state of the economy, including the likelihood of a recession. Among existing studies, the number of time series analyzed is often limited compared with the NBER's information set. Also, in most cases, only aggregate variables are included in the analysis. This section departs from the existing literature by applying the proposed models to disaggregated time series. As shown below, this approach yields sharper inference on the latent states than is typically available in the existing literature.
This subsection describes the data and previews the model. We use five sets of disaggregated series: components of real GDP, industrial production by industry, capacity utilization by industry, nonfarm payroll employment by sector, and average weekly hours by sector. There are 94 series in total. The detailed sectoral and industry breakdowns are reported in the Appendix. We use seasonally adjusted quarterly data. Log differences are used to obtain percentage changes.
Because some disaggregate series begin in 1972 (e.g., capacity utilization for wood products and for computer and electronic products), the sample spans 1972–2019. We first exclude the COVID period so that it does not dominate the estimation. We then extend the sample to include it, using winsorization to limit its influence on estimation (more details are given below).
For illustration, Figure (ref) plots the cross-sectional distribution of the disaggregated series as a fan chart at each point in time, with shaded areas corresponding to NBER recessions. The series exhibit clear comovement, as shown by the quantiles relative to the median, indicating the presence of common cycles. The variability of the series is noticeably higher during NBER recessions than during expansions, and is particularly high during severe recessions. This confirms the importance of allowing the variances of the series to switch across regimes. From an identification viewpoint, it implies that changes in variances contain useful information for identifying business cycles, a point that is not always recognized in the literature.
We first estimate models without autoregressive components (Models A and B with lag order set to zero), first on the pre-COVID sample (1972:II–2019:III) and then on the full sample through 2024:IV. Then we consider models with lagged dependent variables, with the lag order set to one, and data divided into two groups: GDP components and the rest. The variance is always allowed to switch. As a benchmark, Figure (ref) reports smoothed recession probabilities from FRED (series RECPROUSM156N) based on the model of Chauvet (1998), which is representative of estimates typically reported in the literature. These probabilities are at the quarterly frequency, constructed from monthly data using the end-of-period method.
Figure (ref)(a) shows the smoothed recession probabilities together with NBER recession dates (peak to trough quarter) for the 1972-2019 sample. The model detects all NBER-defined recessions. With very few exceptions, the probabilities are either zero or one, showing high posterior certainty about the prevailing regime. This contrasts with the probabilities reported by the FRED series, which can linger at intermediate values for multiple periods. This difference confirms the informativeness of the 94 disaggregate series in identifying the latent state.
At recession onsets, the model dates are typically close to the NBER dates, with four of the six onsets identified within one quarter: the 1980, 1981--82, and 1990--91 recessions are each identified one quarter after the NBER onset; the 2001 recession onset is dated to 2001 Q1 by both. The two exceptions are the 1973--75 recession, where the probability switches to one in 1974 Q3, three quarters after the NBER onset of 1973 Q4, and the 2008--09 recession, where the probability reaches one in 2008 Q2, two quarters after the NBER onset of 2007 Q4. These onset lags may reflect the nature of the underlying data: the early quarters of a downturn are often gradual, and the model achieves high certainty only when a sufficient number of the 94 series have moved in a consistent direction.
The model provides sharp identification of recoveries. In five of the six recessions, the probability of the recession regime drops to (near) zero exactly when the NBER expansion begins. The only exception is the 1975 recovery, where the probability remains at one in 1975 Q2 and drops to zero in 1975 Q3, a one-quarter lag relative to NBER dating. This asymmetry (sometimes moderate delay at onsets, sharp identification of recoveries) is consistent with the view that recessions build gradually across sectors and series, while recoveries tend to be more synchronized and broad-based, which generates a clearer signal in the cross-section of disaggregate data.
We next include the COVID period in the sample. This period is associated with movements in many of the 94 disaggregate series that were extreme by historical standards. This raises the issue of how to handle extreme values. To address this, we winsorize observations in the post-2019 sample that fall outside the 1972–2019 range to the corresponding historical maximum or minimum. That is, if an observation in the extended period exceeds the historical maximum, it is replaced by that maximum; if it falls below the historical minimum, it is replaced by that minimum; otherwise, it is unchanged. The motivation is that for recession detection, extreme values carry little additional information beyond what is already conveyed by a large observation: once a series has moved far enough to be clearly consistent with a recession regime, trimming it has no practical effect on regime classification. At the same time, the adjustment can prevent a small number of pandemic observations from dominating the likelihood, so that parameter estimates would not be pulled away from values that fit the majority of the history.
The results are shown in Figure (ref)(b). The results for the 1972-2019 period are essentially the same. For the COVID period, the NBER dates the recession from 2019 Q4\footnote{This is based on the NBER's quarterly dating. The committee noted that the monthly peak (February 2020) occurred in a different quarter (2020 Q1) than the quarterly peak; thus the dating differs across frequencies. See \url{https://www.nber.org/news/business-cycle-dating-committee-announcement-june-8-2020}.} through 2020 Q2. The model dates the recession from 2020 Q1 to 2020 Q2, a one-quarter onset lag in line with most earlier episodes. The post-2020 period through 2024 Q4 is assigned probability zero throughout, including the period of elevated inflation and monetary tightening in 2022-23 which some contemporaneous observers flagged as recession risk.
In summary, the model identifies all NBER-defined recessions in the sample, including the mild recessions of 1990--91 and 2001, the 2008--09 Great Recession, and the COVID recession. Given the simple model specification, this performance likely reflects the informativeness of the disaggregated data. We next turn to richer models that allow for autoregressive components.
We first consider model ((ref)) and then ((ref)), each on the subsample and then the full sample.
The smoothed recession probabilities for the 1972--2019 sample are shown in Figure (ref)(a) and can be compared with the baseline Figure (ref)(a). Again the model detects all NBER-defined recessions with high posterior certainty, and the onset and exit features are similar to the baseline. This reinforces the view that the cross-sectional information from the 94 disaggregate series is the main driver of regime identification.
There are only three differences from the baseline case. First, the current model produces a recession probability of 0.982 in 1978 Q1, with no counterpart in the baseline. This likely reflects that, with the autoregressive dynamics absorbing some variation in the data, the model becomes more sensitive to short-term movements in the data. In 1978 Q1, industrial production momentum was somewhat subdued relative to 1975-77, and it was interpreted by the model as a brief recession signal. Second, the 2001 onset is somewhat slower: whereas the baseline model assigns probability one immediately in 2001 Q1, the current model produces only 0.431 that quarter, reaching one in 2001 Q2. Third, the 2008–09 recession exit is sharper: the probability drops to zero in 2009 Q3, with no mass around 0.1 as in the baseline case. Other than these, the onset and exit dates coincide with the baseline.
The results for the full sample are reported in Figure (ref)(b) and can be compared with the baseline results in Figure (ref)(b). The smoothed probabilities are close for most of the sample, with only a few differences. The 2008-09 recession exit is now sharper: the baseline assigns a residual probability of 0.10 in 2009 Q3, while the current model assigns zero, matching the start of the NBER expansion. The 2001 onset is one quarter slower than the baseline, with probability 0.02 in 2001 Q1, reaching one in 2001 Q2. Similarly, the COVID-period onset is one quarter slower: the baseline reaches one in 2020 Q1 (a one-quarter lag relative to the NBER onset of 2019 Q4), while the current model assigns probability 0.01 in 2020 Q1 and reaches one in 2020 Q2. Overall, adding the AR dynamics leaves the recession dating conclusions largely unchanged.
We next consider Model B with lag order one. The results for the 1972-2019 sample are reported in Figure (ref)(a) and can be compared with the baseline and the Model A results discussed above.
The estimates remain close to the baseline, with all NBER-defined recessions identified and broadly similar timing. Compared with the baseline, Model B is more responsive to cyclical signals and the onset lags relative to NBER are shorter in several cases. The 1980 recession is identified in its first NBER quarter (1980 Q1) with probability 0.85, compared with near-zero in the other specifications. For 2008-09, although the recession probability still reaches one in 2008 Q2 (as do the other specifications), the model assigns probability 0.38 in 2008 Q1, closer to the NBER dating.
The trade-off is a tendency to assign elevated probabilities outside NBER recession dates: 0.56 in 1990 Q2 and 0.99 in 1981 Q1 during the short inter-recession expansion. The pattern suggests that Model B trades some posterior certainty for shorter onset lags over this period.
The results for the full sample are reported in Figure (ref)(b). The estimates are close to the no-AR baseline (Figure (ref)(b)), with onset lags and exit dates nearly identical across recessions. The only exception is the 1980 exit, where the probability drops to zero in 1980 Q3, one quarter earlier than the baseline (which drops to zero in 1980 Q4, matching the start of the NBER expansion). Notably, the earlier-onset signals and elevated mid-expansion probabilities seen in the short sample (1974 Q1, 1981 Q1, 1990 Q2, and 2008 Q1) do not appear here. The more responsive behavior of Model B in the short sample appears to be sample-dependent: with the longer sample, the estimates move toward the no-AR baseline in Figure (ref)(b).
The above results show that across all model specifications and sample periods, the smoothed probabilities are close to zero or one with few intermediate values. This contrasts with typical recession probability estimates, which tend to linger at intermediate values for extended periods.
Models with autoregressive dynamics show higher sensitivity to cyclical signals, which in the short sample produces more variation: for Model A the elevated probability of 0.98 in 1978 Q1, and for Model B elevated probabilities of 0.83 in 1974 Q1, 0.56 in 1990 Q2, and 0.99 in 1981 Q1 during the short inter-recession expansion. Over the full sample, the estimates across specifications are quite consistent, with differences limited to only isolated quarters. The consistency across specifications confirms that the disaggregate data, rather than the model dynamics, drive regime identification.
The model dates either coincide with NBER onsets or lag them by one to three quarters depending on the episode, while exits coincide closely with the NBER expansion start in most cases. This asymmetry may be a reflection of the nature of recessions: downturns build gradually across sectors, while recoveries tend to generate a more synchronized signal across disaggregate series. The degree of agreement achieved is interesting: the model detects all NBER-defined recessions across more than five decades, despite the fact that NBER dating itself involves a significant judgmental component, which no purely statistical model is designed to replicate.
In terms of computational costs (on an Intel Xeon CPU at 2.80 GHz using 8 cores, with 8 chains producing a total of 0.2 million draws), the computational times are: 0.6 hours for the no-lag specification, 3.1 hours for Model A with one lag, and 0.6 hours for Model B with one lag. Model A takes the most time due to the additional state dependence. Together with the simulation cases reported above, these values confirm that the procedure is computationally feasible for empirically relevant model and sample sizes.
We showed that cross-sectional aggregation attenuates regime-switching signals, and that working with disaggregate data improves state inference under empirically relevant conditions. We then presented models and estimation methods for uncovering regime switching directly from disaggregated data. Applied to a U.S. macroeconomic dataset of 94 series, the models produce nearly binary recession probabilities across all specifications and sample periods, detecting all NBER-defined recessions from 1972 through 2024. The estimated recession onset dates either coincide with or lag the NBER by one to three quarters, while exits closely track the start of the NBER expansion. Applying the models to other datasets (such as European data) or adapting the sampler to related high-dimensional regime-switching models could be useful directions for future research.
\baselineskip=22pt
\setcounter{section}{0} \setcounter{figure}{0} \setcounter{table}{0} \setcounter{equation}{0}