2.1. A regime-switching VAR

Consider a simple VAR model where the variables of interest, $y_{t}$ , exhibit regime-switching behavior:

(1) \begin{align} \begin{array}{c} y_{t}=\beta _{S_{t}}^{\left(0\right)}+\beta _{S_{t}}^{\left(1\right)}y_{t-1}+\cdots +\beta _{S_{t}}^{\left(p\right)}y_{t-p}+u_{t},\quad u_{t}\sim N\left(0, \Sigma_{{S_{t}}}\right),\end{array} \end{align}

where $y_{t}$ is an $n\times 1$ vector of data at time $t ;\, \beta _{S_{t}}^{(0)}$ is an $n\times 1$ vector of constants and $\beta _{S_{t}}^{(l)}$ is an $n\times n$ matrix associated with the lagged terms; and the error term, $u_{t}$ is an $n$ -dimensional zero mean white noise process with covariance matrix, $\Sigma _{{S_{t}}}$ . Stacking the dependent variables into a $T\times n$ matrix, $y$ , whose $t$ ’th row is $y_{t}^{\prime}$ and the regressors into a $T\times (np+1)$ matrix, $X$ , where the $t$ ’th row is $X_{t}=(1,y_{t-1}^{\prime},\cdots, y_{t-p}^{\prime})$ , respectively, we have

(2) \begin{align} \begin{array}{c} y=X\beta _{{S_{t}}}+u, \end{array} \end{align}

where $u$ is a $T\times n$ matrix of innovations in which the $t$ ’th row is $u_{t}^{\prime}$ and $\beta _{{S_{t}}}=[\beta _{S_{t}}^{(0)},{} \beta _{S_{t}}^{(1)},\cdots, \beta _{S_{t}}^{(p)} ]^{\prime}$ . The evolution of latent ( $R$ -state) Markov state variable $S_{t}$ , is governed by the transition probability matrix, $P$ . Specifically, the $(r,r^{\prime})$ element of $P$ , denoted as $P_{(r,r^{\prime})}$ , represents the transition probability from state $r$ at time $t-1$ to state $r^{\prime}$ at time $t$ , that is, $P_{(r,r^{\prime})}=p(S_{t}=r^{\prime}|S_{t-1}=r)$ .

A novel feature of this model is that the key properties of each regime and the associated (regime-specific) parameters (i.e., $\beta _{{S_{t}}}$ and $\Sigma _{{S_{t}}}$ ) are characterized by their own distinctive spectral properties. Specifically, by postulating a few regimes in the frequency domain such as low-frequency, BC-frequency, and/or high frequency, we assume that the dynamics of data at each time can be represented as one or a combination of these regimes.

Our main objective is to make an inference about the unobserved states as well as the model parameters. To do this, we rely on Bayesian inference, which provides a natural framework for incorporating the frequency-dependent behavior of the data. The following subsections describe our approach to prior specification and outline the estimation algorithm.

2.2. Prior

This subsection discusses the (regime-dependent) priors imposed on the model parameters. We adopt a normal-inverse Wishart prior for the VAR model parameters ( $\beta _{{S_{t}}}$ and $\Sigma _{{S_{t}}})$ , along with the Minnesota prior, which provides shrinkage toward the specified prior means in each frequency regime. In the following, for the sake of notational simplicity, we drop the regime subscript, $S_{t}$ , unless it causes confusion.

Eliciting priors for the VAR coefficients largely draws on the framework of the unobserved components model, widely used in the trend-cycle decomposition literature (Harvey, Reference Harvey1985; Watson, Reference Watson1986; Clark, Reference Clark1987). This approach is useful in that it allows us to impose priors that appropriately capture the key characteristics of different frequency regimes with a small set of hyperparameters.

In a typical unobserved component model, (log of) a macro time series is assumed to consist of three components, trend ( $y_{t}^{\tau }$ ), cyclical ( $y_{t}^{c}$ ), and irregular components ( $y_{t}^{\varepsilon }$ ). First, the trend is usually assumed to evolve as follows:

(3) \begin{align} y_{t}^{\tau }=y_{t-1}^{\tau }+\mu _{t}, \\[-24pt] \nonumber \end{align}

(4) \begin{align} \mu _{t}=\mu _{t-1}+\delta _{t},\quad \delta _{t}\sim N\left(0,\sigma _{\delta }^{2}\right), \end{align}

which implies that the trend follows an integrated random walk, or ARIMA(0, 2, 0):

(5) \begin{align} \begin{array}{c} \Delta ^{2}y_{t}^{\tau }=y_{t}^{\tau }-2y_{t-1}^{\tau }+y_{t-2}^{\tau }=\delta _{t}. \end{array} \end{align}

Building on this, we impose an I(2) process for the VAR coefficients in the low-frequency regime ( $S_{t}=1$ ).Footnote ⁴

Next, the prior for the dynamics of the BC-frequency regime ( $S_{t}=2$ ) is specified using trigonometric functions. Specifically, following Planas et al. (Reference Planas, Rossi and Fiorentini2008) and Jarociński and Lenza (Reference Jarociński and Lenza2018), we assume:

(6) \begin{align} \begin{array}{c} \left(1-2\rho \cos \left(2\pi /\tau \right)L+\rho ^{2}L^{2}\right)y_{t}^{c}=u_{t}, \end{array} \end{align}

where $L$ is the lag operator, $\rho$ denotes the amplitude, and $\tau$ represents the periodicity.Footnote ⁵

Finally, we set a simple white noise process for the fluctuations in the high-frequency regime ( $S_{t}=3$ ) (when additionally considered).Footnote ⁶

To summarize, the prior means for the VAR coefficients are set as follows:

(7) \begin{align} \begin{array}{c} \underline{\beta }_{{S_{t}}}=\begin{cases} \left[0_{\left(1\times n\right)} ;\,2I_{n} ;\, -I_{n} ;\, 0_{\left(n-2\right)p\times n} \right], & \text{for }S_{t}=1;\\[3pt] \left[0_{\left(1\times n\right)} ;\, 2\rho \cos \left(2\pi /\tau \right)I_{n} ;\, -{\rho ^{2}}I_{n} ;\, 0_{\left(n-2\right)p\times n} \right], & \text{for }S_{t}=2;\\[3pt] \left[0_{\left(\left(np+1\right)\times n\right)}\right], & \text{for }S_{t}=3. \end{cases} \end{array} \end{align}

We implement these priors by adding the following (regime-dependent) dummy observations, grouped into several blocks, to the actual data (Sims and Zha, Reference Sims and Zha1998).

The first block is specified as:

(8) \begin{align} \underset{\left(np\times n\right)}{\underbrace{ y_{1}^{D}}}=\left[\underline{\beta }^{\left(1\right)}\odot \underline{S}\cdot 1^{{\lambda _{2}}};\cdots ;\underline{\beta }^{\left(p\right)}\odot \underline{S}\cdot p^{{\lambda _{2}}}\right]/\lambda _{1}, \end{align}

and

(9) \begin{align} \underset{\left(np\times \left(np+1\right)\right)}{\underbrace{ X_{1}^{D}}}=\left[0_{\left(np\times 1\right)}, \text{diag}\left(1^{{\lambda _{2}}},\cdots, p^{{\lambda _{2}}}\right)\otimes \underline{S}/\lambda _{1}\right], \end{align}

where $\underline{\beta }^{(l)}$ is an $n\times n$ matrix of (regime-dependent) prior means of the VAR coefficients at lag $l$ , whose $(i,j)$ element represents the coefficient for variable $j$ in equation $i$ at lag $l ;\, \underline{S}$ is the diagonal scale matrix whose elements are set to the standard deviations of error terms from OLS estimates of the autoregression for each variable in the model; and $\odot$ is the element-by-element multiplication operator.

Since this block implies that, for variable $j$ in equation $i$ at lag $l$ ,

(10) \begin{align} \beta _{\left(i,j\right)}^{\left(l\right)}\sim N\left(\underline{\beta }_{\left(i,j\right)}^{\left(l\right)}, \left(\frac{\lambda _{1}}{l^{{\lambda _{2}}}}\times \frac{\Sigma _{\left(i,i\right)}}{\underline{S}_{\left(i\right)}}\right)^{2}\right), \end{align}

it follows that $\lambda _{1}$ governs the overall tightness of the prior, while $\lambda _{2}$ controls the shrinkage of more distant lags. Specifically, a smaller $\lambda _{1}$ indicates a tighter prior while a larger value of $\lambda _{2}$ implies greater shrinkage for longer lags.

The second block is defined as:

(11) \begin{align} \underset{\left(n\times n\right)}{\underbrace{ y_{2}^{D}}} = \underline{\tilde{\beta }}\overline{y}/\lambda _{3} \quad \mathrm{and} \underset{\left(n\times \left(np+1\right)\right)}{\underbrace{ X_{2}^{D}}}=\left[0_{\left(n\times 1\right)},\overline{y}/\lambda _{3},\cdots, \overline{y}/\lambda _{3}\right], \end{align}

where $\underline{\tilde{\beta }}$ is an $n\times n$ diagonal matrix whose elements are the sum of its own coefficients over lags, that is, $\underline{\tilde{\beta }}_{(i,i)}=\sum _{l=1}^{p}\beta _{(i,i)}^{(l)}$ ; and $\overline{y}$ is an $n\times n$ diagonal matrix whose $i$ ’th element is equal to the sample mean of the variable $y_{i}$ (as in Sims and Zha, Reference Sims and Zha1998). This prior implies:

(12) \begin{align} \sum _{l=1}^{p}\beta _{\left(i,i\right)}^{\left(l\right)}\sim N\left(\underline{\tilde{\beta }}_{\left(i,i\right)}, \left(\frac{\lambda _{3}}{\overline{y}_{\left(i\right)}}\right)^{2}\Sigma _{\left(i,i\right)}\right), \end{align}

indicating that $\lambda _{3}$ controls the tightness of the sum of its own coefficients. This prior, a generalization of the random walk prior, allows for flexible forms of persistence (through autocorrelation), with a smaller $\lambda _{3}$ corresponding to a more informative prior.

The third block is given by:

(13) \begin{align} \underset{\left(n\times n\right)}{\underbrace{ y_{3}^{D}}}=\text{diag}\left(\overline{y}\right)/\lambda _{4}\quad \mathrm{and} \underset{\left(n\times \left(np+1\right)\right)}{\underbrace{ X_{3}^{D}}}=\left[0_{\left(n\times 1\right)}, y_{3}^{D},\cdots, y_{3}^{D}\right], \end{align}

and this prior pushes the variables toward their unconditional mean, ensuring that the series in the VAR share a common stochastic trend. Thus, this prior pertains to co-persistence or co-integration. The tightness of this prior is controlled by the hyperparameter $\lambda _{4}$ , where a smaller $\lambda _{4}$ indicates a tighter prior.

The last two blocks are:

(14) \begin{align} \underset{\left(n\times n\right)}{\underbrace{ y_{4}^{D}}}=0_{\left(n\times n\right)} \quad \text{and }\underset{\left(n\times \left(np+1\right)\right)}{\underbrace{ X_{4}^{D}}}=\left[1_{\left(n\times 1\right)}/\lambda _{5}, 0_{\left(n\times np\right)}\right], \end{align}

and

(15) \begin{align} \underset{\left(n\times n\right)}{\underbrace{ y_{5}^{D}}}=\underline{S}\quad \mathrm{and} \underset{\left(n\times \left(np+1\right)\right)}{\underbrace{ X_{5}^{D}}}=0_{\left(n\times \left(np+1\right)\right)}, \end{align}

and these two blocks relate to the priors on the constant terms and the covariances of the error terms, respectively. All these blocks are collected in $y^{D} = [y_{1}^{D};\, y_{2}^{D};\,y_{3}^{D};\,y_{4}^{D};\,y_{5}^{D}]$ and $X^{D} =[X_{1}^{D};\,X_{2}^{D};\,X_{3}^{D};\,X_{4}^{D};\,X_{5}^{D}]$ , both with a length of $T_{D}$ .

The dummy observation priors are characterized by a set of hyperparameters, $\lambda =(\lambda _{1},\cdots, \lambda _{5}).$ In setting their values, we follow the hierarchical approach of Giannone et al. (Reference Giannone, Lenza and Primiceri2015), by treating these hyperparameters as additional regime-dependent parameters to be estimated.Footnote ⁷ In addition to its flexibility, this approach is useful since there seems no consensus on the appropriate values for frequency-specific priors (especially in models involving regime-switching, as in this study). In particular, given that the degree and pattern of shrinkage and persistence may differ across frequencies, it will be interesting to observe whether there are significant differences in the tightness of these priors across different frequency regimes.

Note that by assuming a normal-inverse Wishart distribution, our prior specification treats each equation symmetrically. To address this potentially restrictive property, several approaches involving asymmetric conjugate priors have been proposed (e.g., Carriero et al. Reference Carriero, Clark and Marcellino2019; Carriero et al. Reference Carriero, Chan, Clark and Marcellino2022; Chan, Reference Chan2022). While these approaches offer greater flexibility and computational efficiency, it is unclear whether this feature would be beneficial in the context of our model. On the one hand, allowing for a flexible structure through an asymmetric prior may allow the data to better reflect equation- or variable-specific information. On the other hand, this flexibility may complicate the estimation of the joint switching behavior of the variables in the frequency domain. We revisit this issue later by examining a model in which the symmetric property of the prior is relaxed.

For the transition probabilities, we impose a Dirichlet distribution, as a conjugate prior, for each column of the transition matrix, $P$ . Specifically, for column $r$ , its density depends on $\alpha _{r}=(\alpha _{r1},\cdots, \alpha _{rR})$ and the prior mean of the transition probability (from regime $r$ to $r^{\prime}$ ) is

(16) \begin{align} P_{(r,r^{\prime})}=\frac{\alpha _{rr^{\prime}}}{\sum _{r^{\prime}=1}^{R}\alpha _{rr^{\prime}}}, \end{align}

Collecting all the parameters in $\theta =(\beta, \Sigma, P,\lambda )$ , the prior structure (for each regime) can be written as:

(17) \begin{align} p\left(\theta \right)=p\left(\beta |\Sigma, \lambda \right)p\left(\Sigma |\lambda \right)p\left(\lambda \right)p\left(P\right). \end{align}

2.3. Regime inference, posterior, and likelihood

This section briefly describes the algorithm for regime inference, posterior distribution of model parameters, and the computation of likelihood. In terms of notation, a variable with a superscript $t$ denotes the entire history of the variable up to time t, for example, $y^{t}=(y_{1},\cdots, y_{t})^{\prime}$ , and variables with a subscript $S_{t}$ denote the actual observations (along with dummy observations) corresponding to the respective regime.

In essence, regime inference amounts to computing the probabilities of each regime at each time $t$ , based on the actual observations, $y^{t}$ ,

(18) \begin{align} \xi _{t|t}=p\left(S_{t}|y^{t},\theta \right), \end{align}

where the $t$ ’th element in $\xi _{t|t}$ represents the filtered probability of regime $r$ , $\xi _{t|t}(r)=p(S_{t}=r|y^{t},\theta )$ . Starting with appropriate initial regime probabilities, $\xi _{0|0}$ , this inference can be made recursively, through prediction and updating steps, as follows.

At the start of time $t$ (conditional on the data up to time $t-1$ ), based on the regime inference for the previous period, $\xi _{t-1|t-1}$ , the predicted probabilities for time $t$ , $\xi _{t|t-1}$ (of each regime) are given by:

(19) \begin{align} \xi _{t|t-1}=P\xi _{t-1|t-1}. \end{align}

Then, updating the information contained in $y_{t}$ yields the filtered probabilities:

(20) \begin{align} \xi _{t|t}=\frac{\xi _{t|t-1}\odot \eta _{t}}{\mathbf{1}^{\prime}\left(\xi _{t|t-1}\odot \eta _{t}\right)}, \end{align}

where the symbol $\odot$ denotes element-by-element multiplication, $\mathbf{1}$ is a vector with all elements equal to 1, and $\eta _{t}$ is a vector whose $r$ ’th element contains the conditional likelihood for regime $r$ :

(21) \begin{align} \eta _{t}\left(r\right) & = p\left(y_{t}|S_{t}=r,y^{t-1},\theta \right) \nonumber \\[3pt] & =\left(2\pi \right)^{-n/2}\left| \Sigma _{{S_{t}}}\right| ^{-1/2}\exp \left\{-0.5\left(y_{t}-\beta _{S_{t}}^{\prime}X_{t}\right)^{\prime} \Sigma _{\mathrm{S}_{\mathrm{t}}}^{-1}\left(y_{t}-\beta _{S_{t}}^{\prime}X_{t}\right)\right\}. \end{align}

As a byproduct of this algorithm, the log likelihood function (of all the observed data) can be calculated as follows:

(22) \begin{align} \mathrm{Log}L\left(y^{T}|\theta \right)=\sum _{t=1}^{T}\log f\left(y_{t}|y^{t-1},\theta \right), \end{align}

where $f(y_{t}|y^{t-1},\theta )=\mathbf{1}^{\prime}(\xi _{t|t-1}\odot \eta _{t})$ .

In addition to these filtered probabilities, we can compute the smoothed probabilities for each regime, $p(S_{t}=r|y^{T},\theta )$ , by incorporating the full sample information, $y^{T}$ . The recursive algorithm for this is summarized as follows:

(23) \begin{align} \xi _{t|T}=\xi _{t|t}\odot \left[P^{\prime}\xi _{t+1|T}\left(\div \right)\xi _{t+1|t}\right], \end{align}

where (÷) denotes element by element division.Footnote ⁸

We simulate the posterior distribution of the parameters, using the algorithm of random walk Metropolis-Hastings within Gibbs. For the hyperparameters, we first draw $\lambda ^{*}$ from the normal proposal distribution, $N(\lambda ^{\prime},c_{\lambda }H_{\lambda }^{-1})$ , for each regime, where $\lambda ^{\prime}$ is the previous draw of $\lambda$ , $H_{\lambda }$ is the Hessian of the negative of the log posterior of the hyperparameters at the peak, and $c_{\lambda }$ is a scaling constant. Set $\lambda =\lambda ^{*}$ with probability $\alpha _{\lambda }$ ; otherwise reject the draw and set $\lambda =\lambda ^{\prime}$ , where the acceptance probability is

(24) \begin{align} \alpha _{\lambda }=\min \left(1,\frac{p_{\lambda }\left(y_{{S_{t}}}|\lambda ^{*}\right)}{p_{\lambda }\left(y_{{S_{t}}}|\lambda ^{\prime}\right)}\right), \end{align}

and the marginal likelihood, $p_{\lambda }(y)$ , is as given by Del Negro and Schorfheide (Reference Del Negro, Schorfheide, Del Negro and Schorfheide2011). The scaling parameter $c_{\lambda }$ is set so that acceptance probability is around 0.3 to 0.4.

Given the result of regime inference, we can select and group the sub-sample of actual data that pertains to a particular regime, $y_{{S_{t}}}$ and $X_{{S_{t}}}$ . Then, the conditional posterior distribution for the VAR parameters (which is normal-inverse Wishart due to the natural conjugate prior) is given by:

(25) \begin{align}vec(\beta _{{S_{t}}})|\Sigma _{{S_{t}}},\lambda, S^{T},y^{T}\sim N\left(vec(\overline{\beta }_{{S_{t}}}), \Sigma _{{S_{t}}}\otimes \left(X_{S_{t}}^{*^{\prime}}X_{S_{t}}^{*}\right)^{-1}\right)\end{align}

and

(26) \begin{align}\Sigma _{{S_{t}}}|\lambda, S^{T},y^{T}\sim iW\left(S_{S_{t}}^{*},T_{S_{t}}^{*}\right),\end{align}

where $y_{S_{t}}^{*}$ and $X_{S_{t}}^{*}$ are the regime-specific actual data, augmented with the dummy observations, that is, $y_{S_{t}}^{*}=[y_{{S_{t}}};\,y_{S_{t}}^{D}]$ and $X_{S_{t}}^{*}=[X_{{S_{t}}};\,X_{S_{t}}^{D}] ;\, \overline{\beta }_{{S_{t}}}=(X_{S_{t}}^{*^{\prime}}X_{S_{t}}^{*})^{-1}X_{S_{t}}^{*^{\prime}}y_{S_{t}}^{*} ;\, S_{S_{t}}^{*}=(y_{S_{t}}^{*}-X_{S_{t}}^{*}\beta _{{S_{t}}})(y_{S_{t}}^{*}-X_{S_{t}}^{*}\beta _{{S_{t}}})^{\prime}$ ; and $T_{S_{t}}^{*}$ is the length of $y_{S_{t}}^{*}$ . In each draw of $\beta _{{S_{t}}}$ , we discard all unstable draws and keep only stable ones.

Finally, the transition probabilities (whose posterior is, conditional on $S^{T}$ , independent of the data and other parameters in the model) are sampled from the following Dirichlet distribution:

(27) \begin{align} p\left(P_{r}|S^{T}\right)\sim Dir\left(\alpha _{r1}+\varphi _{r1},\cdots, \alpha _{rR}+\varphi _{rR}\right), \end{align}

where $r=1,2,\cdots, R$ refers to the column number of the transition probability matrix. The parameter $\varphi _{rr^{\prime}}$ refers to the number of times that regime $r$ is followed by regime $r^{\prime}$ , which can be counted using the draw of $S^{T}$ .

The full MCMC estimation algorithm is summarized as follows.

• • Step 0. Initialize parameters and latent state variables, $\theta ^{(0)}$ and $\xi _{0|0}$ .

• • Step 1. Sample $S^{T}$ from $p(S^{T}|\beta, \Sigma, \lambda, y^{T})$ , as described in equations (20) and (23).

• • Step 2. Draw $P$ from $p(P|S^{T})$ , as described in equation (27).

• • Step 3. Draw $\lambda$ from the normal proposal distribution, $N(\lambda ^{\prime},c_{\lambda }H_{\lambda }^{-1})$ , as described above.

• • Step 4. Draw $\Sigma$ from $p(\Sigma |\lambda, y^{T},S^{T})$ , as described in equation (26).

• • Step 5. Draw $\beta$ from $p(\beta |\Sigma, \lambda, y^{T},S^{T})$ , as described in equation (25).

• • Step 6. Repeat Steps 1–5 for a large number of times.

In step 0, we initialize the VAR parameters at their OLS estimates, that is, $\beta ^{(0)}=\hat{\beta }_{OLS}$ and $\Sigma ^{(0)}=\hat{\Sigma }_{OLS}$ . We set $P^{(0)}$ at its prior mean, $\xi _{0|0}$ at the corresponding ergodic (or unconditional) probabilities of the Markov chain, and $\lambda ^{(0)}$ at their posterior mode.Footnote ⁹ In drawing the unobserved states, we ensure that each regime has a minimum of ten observations.

In producing the results below, we first generate 50,000 draws and discard the first 10,000. From the remaining 40,000 draws, we retain every 20th draw, resulting in 2,000 draws on which our inference is based. According to the relative numerical efficiency (RNE) (Geweke, Reference Geweke and Geweke1992) and Raftery and Lewis (Reference Raftery, Lewis, Raftery and Lewis1992) diagnostics, these draws exhibit good convergence of the Markov chain.

3.1. Data and model specification

In our empirical applications, the frequency-dependent regime-switching VAR model presented in Section 2 is applied to quarterly US data on log GDP (multiplied by 100), inflation rates (measured as percentage changes in the CPI), and the federal funds rate.Footnote ¹⁰ The data are obtained from the Federal Reserve Bank of St. Louis’ FRED database, covering the period from 1959:Q1 to 2023:Q2. To address the issue of outliers during COVID-19, which are associated with extremely high variance that may significantly affect parameter estimates, we include dummy variables in the error covariance matrix for the period from 2020:Q1 to 2020:Q3. This adjustment provides a slightly better fit than the model without the dummies, although the overall results remain largely unchanged.Footnote ¹¹

A few decisions remain regarding the model specifications, such as the number of regimes and lag length. For the lag length, we consider values ranging from $p=2\ \mathrm{to}\ 5$ , covering most of the commonly used options for quarterly VAR models, with a minimum of 2 lags due to prior settings. Regarding the number of regimes, while the trend-cycle decomposition literature often considers three components and thus suggests three frequency regimes (i.e., low, BC, and high) as discussed in Section 2.2, it is not clear a priori whether the third or high-frequency regime is particularly relevant or economically meaningful in the context of our regime-switching VAR model. Therefore, we also explore models with two regimes, combining the high-frequency regime with the BC-frequency regime and setting the prior for the combined regime equal to that of the BC-frequency regime. This results in 8 possible specifications. Based on the log marginal likelihoods, we find that the model with 2 regimes and 4 lags provides the best fit, and this serves as our baseline model.Footnote ¹²

The specific priors are set as follows. For the VAR coefficients in the BC-frequency regime, we follow Planas et al. (Reference Planas, Rossi and Fiorentini2008) and Jarociński and Lenza (Reference Jarociński and Lenza2018) and set $\rho =0.7$ and $\tau =32$ , implying a typical BC lasts 8 years. For the transition probabilities, we set $\alpha _{rr^{\prime}}=8$ if $r=r^{\prime}$ and $\alpha _{rr^{\prime}}=2$ otherwise.Footnote ¹³ This specification ensures that each regime is treated symmetrically and exhibits reasonable persistence, with the continuation probability of $P_{(r,r^{\prime})}=0.8$ . Note that, except for the VAR coefficients, the priors imposed on all other parameters are either the same or symmetric across the regimes.

3.2. Labeling and regime identification

In Bayesian inference for regime-switching models, a practical and important issue is the labeling problem (Frühwirth-Schnatter, Reference Frühwirth-Schnatter2001). Since the likelihood in these models is invariant with respect to permutations of the regime labels, without careful and appropriate identifying restrictions, a model may not be globally identified, which can lead to biased posterior inferences.

To address this issue, a common solution is to impose an identification restriction on the model parameters, such as imposing an order on certain parameters. For models with frequency-domain features, such as ours, another useful approach is to examine the shape of the spectral density. To illustrate this point, Figure 1 plots the spectral densities for each regime implied by these priors, along with a standard (unit root type) Minnesota prior.Footnote ¹⁴

Figure 1.

Spectral densities implied by priors.

Note: The figure presents the spectral densities implied by the regime-specific priors. The vertical lines indicate the frequency band associated with typical business cycles, that is, [ $2\pi /32,2\pi /6]$ .

As shown in the figure, while the spectral densities decrease with frequency (except for the high-frequency regime whose densities are flat), and each regime contains other frequency components beyond its target frequency ranges, each regime appears to have distinctive spectral characteristics in terms of shape. Specifically, the spectral densities are relatively more concentrated within their respective frequency ranges, with varying degrees of decline: the low-frequency regime exhibits the steepest decline, followed by the Minnesota prior, and then the BC-frequency regime. These features allow the construction of simple and reasonable identification criteria based on the properties of spectral densities.

Specifically, once the normalized frequency domain is divided into two parts, $[0, \tilde{\omega }]$ and $[\tilde{\omega },0.5]$ , each representing a relatively low- and high-frequency range, respectively, we can compute the difference in the average log spectral density for regime r between these two ranges, denoted as $D(\overline{\log (SD(r))}=\overline{\log (SD(r))}_{[0, \tilde{\omega }]}-\overline{\log (SD(r))}_{[\tilde{\omega },0.5]}$ . This difference can be roughly interpreted as the “slope” of the density and the identification condition can then be stated in terms of this slope as follows:

(28) \begin{align} D\left(\overline{\log \left(SD\left(1\right)\right)}\right)\gt D\left(\overline{\log \left(SD\left(2\right)\right)}\right). \end{align}

After some experimentation, we find that a cutoff frequency of $\tilde{\omega }=0.083$ (corresponding to a 12-quarter cycle) works well and gives reasonable results.Footnote ¹⁵

3.3. Estimated regimes and their key properties

We now present the main empirical results. Figure 2 shows the smoothed probabilities of the BC-frequency regime from the baseline model, along with the associated spectral densities (the 5%–95% intervals of the estimated log spectrum), while Table 1 contains the estimates of the model parameters. For comparison, the results for the constant parameter VAR (CVAR) are also included in Table 1.Footnote ¹⁶

Table 1.

Estimates for model parameters

Note: The table presents the median estimates for the model parameters (in each frequency regime) along with the 16%–84% credible intervals. $u_{y}$ , $u_{p}$ , and, $u_{i}$ represent the (reduced-form) shocks in the GDP, inflation, and interest rates equations, respectively.

Figure 2.

Estimated regimes and spectral densities.

Note: The left panel of the figure shows the (mean) smoothed probabilities of the BC-frequency regimes. The shaded areas represent NBER recession dates. The right panel shows the 5%–95% intervals of the estimated log spectra for each frequency regime. The vertical lines indicate the frequency band associated with typical business cycles, that is, [ $2\pi /32,2\pi /6]$ .

Figure 2 illustrates clear frequency-dependent regime-switching behavior over time. While the low-frequency regime generally dominates, the period from the 1970s to the mid-1980s is largely characterized by the BC-frequency movements. This pattern of frequency shifts is consistent with the results of other related studies. For example, Aguiar-Conraria et al. (Reference Aguiar-Conraria, Martins and Soares2012) show that, based on wavelet analysis, there are notable changes in the frequency behavior of economic activity, inflation, and interest rates during this period. Similarly, Sargent and Surico (Reference Sargent and Surico2011) show that the low-frequency relationship between inflation and money growth exhibited distinct behavior in the 1970s and early 1980s, while Kliem et al. (Reference Kliem, Kriwoluzky and Sarferaz2016) find a similar pattern in the relationship between public deficits and inflation. Both studies suggest that these low-frequency changes are likely related to shifts in the conduct of monetary policy over time.

While several recession episodes roughly align with the BC-frequency regime, the degree of their association varies across episodes, suggesting that their synchronization is not particularly strong. For example, the 1990-1991 and 2001 recessions are either barely or weakly identified as the BC-frequency regime by the model. Meanwhile, the 1969–1970 recession and the 2007–2009 Global Financial Crisis are identified as the BC-frequency regime, though with some lag.

The spectral densities indicate that, while a large portion of the spectra for both regimes overlaps in the low frequency regions, there is a clear distinction in their slopes. Specifically, the low-frequency regime exhibits a steeper slope and contains relatively more low-frequency components.

The sharp shift in frequency regimes (from BC-frequency to low-frequency) around the mid-1980s is reminiscent of the Great Moderation. As Panel A of Table 1 shows, the BC-frequency regime is indeed associated with higher volatility, consistent with conventional findings. The covariances also show clear differences, with shocks being more highly correlated in the BC-frequency regime. Notably, while the covariance between output and inflation shocks is positive in both the BC-frequency regime and the CVAR, it becomes significantly negative in the low-frequency regime.

Panel B of Table 1 shows that, while the other hyperparameters are roughly similar across the two frequency regimes, there are clear differences in the two specific hyperparameters. Specifically, λ², which controls for shrinkage at longer lags, is much larger in the low-frequency regime, while λ⁴, which is associated with co-persistence, is significantly smaller in the low-frequency regime. These results are intuitive and suggest that: (i) the low-frequency regime is associated with greater persistence, possibly in the form of an integrated process, where tighter shrinkage at longer lags may improve the representation of the dynamics in this regime; and (ii) co-persistence or co-integration properties between variables are more suitable for the low-frequency behavior of the series.

The transition probabilities (Panel C in Table 1) show strong persistence for both regimes, with the low-frequency regime having a higher median continuation probability of 0.955, compared to 0.873 for the BC-frequency regime.

To further illustrate the regime-dependent dynamics implied by the model, Figure 3 presents the impulse responses across the regimes, along with those from the CVAR. The identification is based on a simple Cholesky decomposition, with the ordering of output, inflation, and interest rate (the last being the most exogenous). While it is somewhat unclear to discern the economic implications for some of the responses under this identification scheme, the figure overall reveals sharply contrasting results across the regimes. Specifically, the responses in each regime are characterized by strong persistence and clear cyclical patterns, respectively, with the responses from the CVAR falling between those of the two regimes. Beyond these differences, there are several notable features across the regimes, which we briefly discuss below.

Figure 3.

Impulse responses across frequency regimes.

Note: The figure shows the 16%–84% bands of the impulse responses across frequency regimes (blue solid line), along with those of the constant parameter VAR (red dashed line).

First, in the case of an interest rate shock (which can be viewed as a rudimentary form of a monetary policy shock), both output and inflation increase persistently and significantly in the low-frequency regime. In contrast, in the BC-frequency regime, output decreases significantly for up to six years, while the responses of inflation, aside from an initial short-term increase, remain largely insignificant over a 2- to 10-year horizon. While the literature often attributes the price puzzle, that is, the strong positive response of prices following an interest rate shock, to missing information or omitted variables in VARs, this study suggests that another contributing factor may be the failure to properly extract relevant frequency-domain information about the shock and the associated variable dynamics.

Next, following an inflation shock, while output declines persistently and significantly in the BC-frequency regime (much more so than in the CVAR), its responses in the low-frequency regime are much smaller and become insignificant after the first several quarters. Finally, for an output shock, while the responses of output in the low-frequency regime are strongly persistent (more so than those in the CVAR), the responses in the BC-frequency regime are relatively transitory and become insignificant after a few years. Interestingly, in the low-frequency regime, the short-run responses of inflation to an output shock are negative, although they quickly turn positive after a few quarters, in contrast to the BC-frequency regime and the CVAR, where its short-run responses are significantly positive. This difference in the initial responses reflects the results of the regime-dependent error covariance matrix, as reported in Table 1.

4.1. List of alternative models

This subsection briefly outlines the main features of the alternative models, with detailed descriptions provided in the Appendix. First, to assess the role of regime-dependent heteroskedasticity in the disturbances, we consider regime-switching VAR models with either homoskedastic shocks or stochastic volatility. In the homoskedastic model (denoted as RS-VAR-homo), the error covariances are constrained to be the same in both regimes (i.e., $\Sigma _{{S_{t}}}=\Sigma$ for all $t$ ), while all other features, including the prior means of the VAR coefficients, remain the same as in the baseline model (hereafter denoted as RS-VAR-FD).

For the model with stochastic volatility (denoted as RS-VAR-SV), we follow the approach outlined in the literature (e.g., Primiceri, Reference Primiceri2005). In this model, the (co)variances of the innovations are assumed to (gradually) evolve over time rather than switching between regimes while all other features remain the same as in the baseline model. Specifically, the model relies on the standard decomposition of the (time-varying) error covariance matrix, $\Sigma _{t}=A_{t}^{-1}H_{t}A_{t}^{-1^{\prime}}$ , where $A_{t}$ is a lower triangular matrix and $H_{t}$ is a diagonal matrix. We assume that the evolution of the diagonal elements of $H_{t}$ and the elements of $A_{t}$ follows a (geometric) random walk.

To address the issue of symmetry across the equations associated with the natural conjugate prior, discussed in Section 2.2, we also consider a model where this restrictive structure is relaxed. In specifying this model (denoted as RS-VAR-AP), we largely follow the approach of Carriero et al. (Reference Carriero, Clark and Marcellino2019) and Carriero et al. (Reference Carriero, Chan, Clark and Marcellino2022) and allow each equation to be treated and estimated individually. Thus, in this model, each equation as well as each regime is governed by a set of separate hyperparameters.

We also estimate a conventional regime-switching VAR model with standard time-domain features (denoted as RS-VAR-TD) to identify expansion and recession regimes. While retaining the same essential structure as the RS-VAR-FD, this model differs in several key respects. First, for the VAR coefficients, following the literature (e.g., Sims and Zha, Reference Sims and Zha2006; Sims et al. Reference Sims, Waggoner and Zha2008), we impose a unit root-type Minnesota prior on the VAR coefficients in both regimes. With this symmetric prior, the dynamic properties of each regime are primarily determined by the data. Second, regime identification in this model relies on typical time-domain characteristics of recessions, such as lower output growth and higher volatility. Thus, the associated identification condition (in terms of model parameters) is that, in the recession regime, the intercept in the GDP equation is lower and the determinant of the error covariance matrix is larger.

Diebold et al. (Reference Diebold, Lee, Weinbach, Diebold, Lee and Weinbach1994) and Filardo (Reference Filardo1994) argue that constant transition probabilities in regime-switching models may be too restrictive for many empirical contexts. In the spirit of these studies, we also examine models with time-varying transition probabilities which are specified as functions of lagged output growth, in both the baseline and conventional (time-domain) models (each denoted as RS-VAR-FD-TVTP and RS-VAR-TD-TVTP, respectively). Given that both the transition probabilities and the regime inferences in these models are based on past output dynamics, this approach may provide additional insights into the nature of regime-switching dynamics.

The alternative models are summarized as follows:

• • RS-VAR-homo: Frequency-dependent regime-switching VAR model with homoskedasticity

• • RS-VAR-SV: Frequency-dependent regime-switching VAR model with stochastic volatility

• • RS-VAR-AP: Frequency-dependent regime-switching VAR model with an asymmetric prior

• • RS-VAR-TD: Conventional time-domain regime-switching VAR models with constant transition probabilities

• • RS-VAR-TD-TVTP: Conventional time-domain regime-switching VAR models, both with constant and time-varying transition probabilities

• • RS-VAR-FD-TVTP: Baseline frequency-dependent regime-switching VAR model with time-varying transition probabilities

4.2. Comparison with alternative specifications in variances and priors

The main results for these alternative models are summarized in Figures 4 and 5, which show the smoothed regime probabilities and the associated spectral densities, and in Tables 2 and 3, which show the estimates of selected model parameters.

Figure 4.

Results for alternative models.

Note: Each panel shows the result for the three alternative models. Panel (a) The model with homoskedasticity; panel (b) The model with stochastic volatility; panel (c) the model with asymmetric prior. See notes to Figure 2.

Figure 5.

Results for the conventional (Time-domain) regime-switching model (RS-VAR-TD). Notes: See notes to Figure 2. The figure plots the 5%–95% intervals of the estimated log spectra for the low-frequency and the expansion regime (left panel) and the BC-frequency and the recession regime (right panel).

In the case of the RS-VAR-homo (Figure 4(a)), the number of periods identified as belonging to the BC-frequency regime is significantly smaller compared to the RS-VAR-FD; many periods previously classified as part of the BC-frequency regime in the RS-VAR-FD are now identified as belonging to the low-frequency regime. This is because, unlike in the RS-VAR-FD (where the likelihood is relatively flatter and more dispersed in the BC-frequency regime, while more concentrated in the low-frequency regime), imposing homoskedasticity increases the likelihood of a typical shock or data point favoring the low-frequency regime. In other words, data that might have been identified as belonging to the BC-frequency regime in the RS-VAR-FD are now more likely to be classified as belonging to the low-frequency regime.

In the case of the RS-VAR-SV (Figure 4(b)), although some weak patterns of frequency switching are observed during the 1970s to 1980s, the distinction between regimes is less clear. This result seems to be due to two factors. First, the regime-specific errors in this model are likely to have non-normal and asymmetric distributions. As noted by Hwu and Kim (Reference Hwu and Kim2024), when the distribution of the error term deviates from normality, especially in cases of asymmetry, the accuracy of regime probability inference deteriorates significantly. In addition, the lack of tightly synchronized heteroskedasticity across the series (as shown in Appendix Figure A3) may further complicate clear regime inference.Footnote ¹⁷

Figure 4(c) shows the results for the RS-VAR-AP. Similar to the RS-VAR-homo and the RS-VAR-SV, this model does not provide clear evidence of distinct frequency switching. As with the RS-VAR-SV, the lack of a clear distinction between regimes is likely because the dynamics of each series are not tightly synchronized in terms of their frequency-domain behavior. To illustrate this, we apply the baseline model to each series, with the estimated regimes shown in Appendix Figure A4. Each series exhibits distinct frequency-dependent switching patterns, with common movements observed only during limited periods.

Table 2.

Estimates of parameters for alternative models

Note: The table presents the median estimates for parameters in each frequency regime, along with the 16%–84% credible intervals. See the notes to Table 1.

Table 3.

Estimates of parameters for the RS-VAR-TD

Note: The table presents the median estimates for parameters in each regime (expansion and recession), along with the 16%–84% credible intervals. See the notes to Table 1.

As shown in the figure, in all these models, the spectral densities of the two regimes are not clearly separated and tend to overlap across the entire frequency range, which is consistent with the poor or unclear regime inference.

4.3. Comparison with conventional regime-switching models

We now compare the results of the RS-VAR-TD (shown in Figure 5 and Table 3) with those of the RS-VAR-FD. While the overall results appear similar in terms of the parameter estimates and the estimated latent regimes, a closer inspection reveals several notable and important differences.

First, the identified recession regimes (in the RS-VAR-TD, Panel (a) of Figure 5) are more closely aligned with the NBER reference cycles. For example, while the RS-VAR-FD classifies the periods of the 1990–1991 and 2001 recessions as either in the low-frequency regime or in between the two frequency regimes, respectively, the RS-VAR-TD identifies these periods as the recession regime with much higher probabilities. A similar observation can be made for the early periods of the 1969–1970 and 2007–2009 recessions.

The differences in regime inference between the two models can be attributed to their different approaches. In the RS-VAR-TD, regime identification is essentially based on the first two moments of the data (i.e., mean and volatility), whereas in the RS-VAR-FD, the transition dynamics around recessions play a more important role in regime inference. For example, if the transition into and out of a recession is relatively smooth and entails more low-frequency features, the recession period is more likely to be classified as a low-frequency regime.Footnote ¹⁸

In this regard, comparing the spectral properties of the two models reveals further differences. Panel (b) of Figure 5 compares the spectral densities of both models in each regime (i.e., low-frequency/expansion and BC-frequency/recession, respectively). While the shapes of the spectral densities for the second regime (BC-frequency/recession) are largely similar across the two models (with slightly wider intervals in the RS-VAR-FD), there are notable differences in the spectral densities for the first regime (low-frequency/expansion). Specifically, the spectral densities in the expansion regime are relatively flatter compared to those in the low-frequency regime. This suggests that, although both regimes contain a substantial amount of low-frequency components, the low-frequency regime contains relatively more low-frequency components, implying relatively more persistence, whereas the expansion regime features relatively more high-frequency components, indicating relatively less persistence.Footnote ^{19 ,} Footnote ²⁰

In addition, as shown in Figure 5, although the spectral densities in the RS-VAR-TD do not show clear differences in slope between the two regimes, a substantial portion of the spectra in the recession regime appear to lie above those in the expansion regime. This suggests that, in a rough sense, the spectral densities in the recession regime can be seen as a parallel upward shift from those in the expansion regime. In other words, the differences in the spectral characteristics between the two regimes in the RS-VAR-TD are influenced more by the error covariances than by the VAR coefficients.

Given the several differences between the two types of regime-switching models, an interesting question arises: What would happen if we applied the recession identification criteria (in the time-domain model) to the RS-VAR-FD? Conversely, what would happen if we applied the frequency-domain identification criteria (based on spectral density shapes) to the RS-VAR-TD? Exploring these questions could provide valuable insights into the role of prior specifications and regime identification schemes, potentially revealing further differences between the two types of models.

Figure 6 shows the results of these exercises. When the (time-domain) recession identification is applied to the RS-VAR-FD, the identified regimes are broadly similar to those in the original RS-VAR-FD (Figure 6(a)). This is not surprising, as the BC-frequency regime is also associated with higher volatilities, with the minor differences due to the additional consideration of the intercept term in the GDP equation. In contrast, applying the frequency-domain identification to the RS-VAR-TD results in poor regime inference, making the results difficult to interpret (Figure 6(b)). This result is largely due to the fact that, as discussed earlier, the two regimes in the RS-VAR-TD model are not clearly distinguished in terms of the slope of the spectral densities.Footnote ²¹

Figure 6.

Results from alternative identification.

Note: In each panel, the figure shows the probability of the BC-frequency (panel (a)) or recession regime (panel (b)).

These results suggest that, despite some similarities between the two types of models in terms of the volatilities and persistence of each regime, one cannot directly associate expansion/recession periods with low-/BC-frequency regimes, as each of their spectral properties appears to be different between the two models. Moreover, because of these differences, natural frequency-domain identification cannot simply be applied to identify the recession regime as in conventional time-domain models; careful and model-specific identification is required in each case.

Finally, when the transition probabilities are allowed to vary over time, the estimated regimes remain broadly similar to those obtained with constant transition probabilities (Appendix Figure A7). However, the distinction between regimes becomes clearer and more consistent with the intended regime properties of each model, which is also reflected in a slightly improved model fit, as indicated by the log marginal likelihoods.

The spectral densities for these models show that, in the second regime (BC-frequency/recession), there are notable peaks around 0.1 for the BC-frequency and 0.3 for the recession regime, corresponding to 10-quarter and 4-quarter cycles, respectively. This suggests that time-varying transition probabilities introduce more uncertainty around these frequency ranges, where the associated frequencies are also relatively lower in the frequency-dependent switching model (RS-VAR-FD-TVTP).

Table 4.

Forecast comparison

Notes: At each forecast horizon, the row corresponding to CVAR contains the RMSEs and the average log predictive likelihood (in italics) while other entries are values relative to them. Asterisks (∗) and daggers (†) indicate the significance of the Giacomini and White (Reference Giacomini and White2006) predictive ability test against the CVAR at the ∗ 0.10 and ∗∗ 0.05 levels, and against the RS-VAR-FD at the † 0.10 and †† 0.05 levels, respectively.

4.4. Forecasting performance comparison

Finally, this subsection presents a simple comparison of forecasting performance. Although the main objective of the study is not to examine the predictive ability of the model, evaluating the predictive performance of econometric models has become standard practice, and the forecast evaluation can be thought of as a useful tool for model comparison.

Therefore, we construct a small pool of models for forecasting comparison, including CVAR, RS-VAR-FD, RS-VAR-TD, RS-VAR-TVTP, and RS-VAR-TD-TVTP.Footnote ²² For all models, we use a recursive estimation window to generate pseudo out-of-sample forecasts for $h=1$ , 2, 3, and $4$ -step-ahead horizons. Starting with the estimation sample up to 1983:Q4, the first available forecasts are for 1983:Q4+h. We then proceed by expanding the sample by a quarter and generating forecasts through the end of the sample. In line with conventional practices and the interests of policymakers, h quarter-ahead forecasts for GDP and inflation are the cumulative quarterly rates of change (i.e., the percentage change from period t through t + h). All the forecasts from each model are transformed to be consistent with this cumulative data measure. The exercises involve multi-step forecasting, and we use iterated forecasting to obtain the result, which generally outperforms direct forecasting (Marcellino et al. Reference Marcellino, Stock and Watson2006). In evaluating the forecasting performance of the models, we use the root mean squared error (RMSE) and the predictive likelihood as evaluation metrics. In the latter, we measure the log-predictive score using a Gaussian approximation of the predictive density.

The results of this forecast comparison are summarized in Table 4. We normalize all the results in the table relative to the CVAR. Specifically, the numbers represent the ratio of each model’s RMSE to that of the CVAR and the average of the log predictive likelihoods of a given model minus that of the CVAR, respectively. In addition, to provide a rough indication of the statistical significance of the differences in predictive performance, we report the results of the predictive ability test (Giacomini and White, Reference Giacomini and White2006) for each model against both the CVAR and the RS-VAR-FD.

As shown in the table, forecasting with the RS-VAR-FD leads to substantial improvement over the CVAR in most cases. The improvement is more pronounced for inflation and interest rates, as well as for the density forecast over short horizons, and is often accompanied by strong statistical significance over the CVAR. Compared to the RS-VAR-TD, the overall performance of the RS-VAR-FD is somewhat worse, especially for point forecasts. However, the predictive ability test results show that their performances are not significantly different in most cases and that the RS-VAR-FD remains strong and competitive over short horizons ( $h=1$ and $h=2$ ). Despite their better in-sample fit, the models with time-varying transition probabilities perform slightly worse than their constant probability counterparts in several cases, indicating the possibility of overfitting.