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.
139,516 characters · 15 sections · 0 citation commands
Estimation and Inference for High Dimensional Factor Model with Regime Switching
\baselineskip=20.0pt\ \ \ \setcounter{page}{0}
A great deal of attention has focused on the loading instability issue in high dimensional factor models. For empirical evidences of parameter instability in macroeconomic and financial time series, see for example, Banerjee, Marcellino and Masten (2008), Stock and Watson (2009) and Korobilis (2013). Several procedures are proposed to detect and/or estimate common abrupt breaks in the loadings, including Cheng, Liao and Shorfheide (2016), Baltagi, Kao and Wang (2017, 2021), Bai, Han and Shi (2020), and Ma and Su (2018), to mention a few. Other models of time varying loadings, such as i.i.d./random walk, smooth change, vector autoregression and threshold type, are studied in Bates, Plagborg-Moller, Stock and Watson (2013), Su and Wang (2017), Mikkelsen, Hillebrand and Urga (2019) and Massacci (2017), respectively.
An alternative approach of modeling loading instability is common regime switching. In business cycle analysis, several unobservable factors summarize the comovements of many economic variables and the loadings measure the importance of factors for each economic variable. The importance of each factor may be different depending on fiscal policy (expansionary, contractionary, neutral), or monetary policy (expansionary, contractionary), or the stage of the business cycle (peak, trough, expansion, contraction), hence the loadings may switch synchronously between several states under different scenarios. In stock return analysis, the loadings measure the impact of the factor return on the expected return of each individual stock, hence the loadings may switch synchronously depending on the stock market scenarios (bull versus bear markets, high versus low volatility), see for example Gu (2005) and Guidolin and Timmermann (2008) for related discussions. In bond return analysis, the yields of bonds with different maturities are well captured by the level factor, the slope factor and the curvature factor, see for example Cochrane and Piazzesi (2005) and Diebold and Li (2006). The importance of each factor could be different depending on the stock market volatility, or the stage of the business cycle, or the unemployment rate, hence the loadings may also switch synchronously according to these state variables. In general, large factor models with regime switching in the loadings could also be useful for other topics, such as tracking labor productivity.
There are only a few related results on large factor models with regime switching in the loadings. Liu and Chen (2016) proposes an iterative algorithm for estimating the model parameters and the hidden states based on eigen-decomposition and the Viterbi algorithm, however, the asymptotic properties of the estimated parameters are established only when the true states are known. Considering loadings as general functions of some recurrent states, Pelger and Xiong (2021) develops nonparametric kernel estimator for the loadings and the factors, and establishes the relevant asymptotic theory. However, Pelger and Xiong (2021) requires observable state variables. In general, state variables may be misspecified or unobservable.
This paper proposes maximum (quasi)likelihood estimation for high dimensional factor models with regime switching in the loadings when the state variables are unobservable. This paper also proposes new criteria to consistently determine the number of regimes and the number of factors in each regime. The model parameters are estimated jointly by the EM algorithm, which in the current context only requires calculating principal components iteratively.
More specifically, in the E-step, the probabilities of each regime at each time $t$ are calculated based on the observed data and the parameter values at the current iteration using a recursive algorithm modified from Hamilton (1990), and then the joint (log)likelihood of the observed data and the unobserved states are averaged with respect to the calculated regime probabilities. In the M-step, the estimated loadings for each regime are the principal components of the weighted sample covariance matrix of the observed time series, where the weight on $x_{t}$ (the observed time series at time $t$) equals the probability of that regime at time $t$. Since principal components can be easily calculated even when $N$ (the dimension of time series) is large, our method is very easy to implement.
For the proposed algorithm, this paper establishes the convergence rates of the estimated loading spaces and the estimated factor spaces, the limit distributions of the estimated loadings and the estimated factors, the consistency of the estimated regime probabilities, and the consistency of the estimated transition probability matrix when the true state process is Markovian. Note that asymptotic analysis under the regime switching setup is more difficult than under the structural break setup, because the pattern of regimes for the latter is much simpler.
These asymptotic results are essential in many empirical contexts. First, the limit distributions of the estimated factors allow us to construct confidence intervals for the true factors, which represent economic indices in many applications. The result on the estimated factor spaces implies that if the estimated factors are used in factor-augmented forecasting (or factor-augmented VAR), the forecasting equation (or the VAR equation) would have induced regime switching in the model parameters. Second, for asset management, the estimated loadings of each regime allow us to construct portfolios according to each specific market scenario. For structural dynamic factor analysis, consistently estimated loadings are also crucial for recovering the impulse responses. Third, the consistency of the estimated regime probabilities implies that for each $x_{t}$, we can consistently identify which regime $x_{t}$ belongs to as $N\rightarrow \infty $. For asset management, this allows us to consistently identify the current market scenario. For business cycle analysis, this allows us to consistently date turning points of the business cycle and detect new recessions or expansions, especially when high frequency (weekly, daily) data is utilized.
For cases with small $N$, various methods have been proposed for estimating factor models with regime switching. Kim (1994) proposes approximate Kalman filter for likelihood evaluation and uses nonlinear optimization for likelihood maximization. Kim and Yoo (1995) and Chauvet (1998) apply Kim (1994)'s method to a small number of economic series and obtain recession probabilities and turning points very close to the official NBER dates. Kim (1994) allows for regime switching in both the factor mean and the factor loadings, but when $N$ is large, Kim (1994)'s method would be very time consuming and may have convergence problems\footnote{ This is because the number of parameters grows proportionally to $N$ and the likelihood function is calculated numerically and maximized by nonlinear optimization algorithm.}. Other methods, such as Diebold and Rudebusch (1996) and Kim and Nelson (1998), assume stable loadings and only focus on regime switching in the factor mean. If the loadings are unstable, these methods are not applicable. More importantly, if there is only regime switching in the factor mean, we can not consistently identify each regime even when $N$ is large.
In contrast with Kim (1994), our method is fast and easy to implement even when $N$ is very large. The crucial point behind our EM algorithm is to ignore factor dynamics\footnote{ Factor dynamics are still allowed for the data generating process. } and integrate out the factors in the likelihood function. If factors dynamics are taken into account or factors are treated as parameters in the likelihood function, the estimated loadings would not be the principal components of the weighted sample covariance matrix, and consequently both the algorithm and the asymptotic analysis would become infeasible. On the other hand, the efficiency loss of ignoring factor dynamics is small when $N$ is large.
This paper may also contribute to the literature on dating turning points of the business cycle. Currently there are two main approaches for dating business cycle using multiple time series. The first approach, aggregating then dating, is to date business cycle by focusing on a few highly aggregated time series such as GDP, industrial production and nonfarm employment. The second approach, dating then aggregating, is to date turning points in each disaggregated series and then aggregate these turning points in some appropriate way, see Burns and Mitchell (1946), Harding and Pagan (2006) and Chauvet and Piger (2008). These papers only use a small number of time series. Stock and Watson (2010, 2014) studies this issue using many time series. This paper shows that it is possible to consistently identify turning points if regime switching is synchronous and $N$ is large enough. If $N$ is small, consistency is not possible no matter how large $T$ is. This paper also shows that if $N$ is large, it is possible to consistently detect regime switching right after the turning point with only one observation, thus the speed of detection could be improved significantly. If $N$ is small, we have to wait for enough observations from the new regime.
The rest of the paper is organized as follows. Section (ref) introduces the model setup and the estimation procedures. Section (ref) presents the assumptions and the asymptotic results. Section (ref) proposes criteria for determining the number of regimes and the number of factors in each regime. Section (ref) presents simulation results.\ Section (ref) presents an empirical application of the proposed method to the FRED-MD dataset. Section (ref) concludes. All proofs are relegated to the appendix.
Through out the paper, $(N,T)\rightarrow \infty $\ denotes $N$\ and $T$\ going to infinity jointly, $\delta _{NT}=\min \{\sqrt{N},\sqrt{T}\}$. $ \overset{p}{\rightarrow }$ and $\overset{d}{\rightarrow }$ denotes convergence in probability and convergence in distribution, respectively. For matrix $A$, let $\left\Vert A\right\Vert $, $\left\Vert A\right\Vert _{F} $, $\rho _{\max }(A)$ and $\rho _{\min }(A)$ denote its spectral norm, Frobenius norm, largest eigenvalue and smallest eigenvalue, respectively. Let $P_{A}=A(A^{\prime }A)^{-1}A^{\prime }$ denote the projection matrix and $M_{A}=I-P_{A}$. "w.p.a.1" denotes with probability approaching one.
Consider the following factor model with regime switching: for $i=1,...,N$ and $t=1,...,T,$
where $\lambda _{ji}^{0}$ is an $r_{j}^{0}$ dimensional vector of loadings for regime $j$, $f_{t}^{0}$ is an $r_{z_{t}}^{0}$ dimensional vector of factors, $z_{t}$ is the state variable indicating which regime $x_{it}$ belongs to, and $e_{it}$ is the error term allowed to have cross-sectional and temporal dependence as well as heteroscedasticity. $x_{it}$ is observable and all of the right hand side variables are unobservable. The number of regimes $J^{0}$ and the number of factors in each regime $ r_{j}^{0} $ (could be different across $j$) are fixed as $(N,T)\rightarrow \infty $ and assumed to be known in this section and Section (ref). How to consistently determine $r_{j}^{0}$\ and $J^{0}$\ will be studied in Section (ref).
The factor process $\{f_{t}^{0},t=1,...,T\}$ is allowed to be dynamic, and similar to the principal component estimator (PCE) in Bai (2003) and the maximum likelihood estimator (MLE) in Bai and Li (2012, 2016), factor dynamics are ignored when estimating the model parameters, thus there is no need to model factor dynamics.
For the state process $\{z_{t},t=1,...,T\}$, the asymptotic results in Section (ref) and Section (ref) are valid as long as $\frac{1}{T} \sum\nolimits_{t=1}^{T}1_{z_{t}=j}\overset{p}{\rightarrow }q_{j}^{0}>0$\ for $j=1,...,J^{0}$\ ($q_{j}^{0}=\Pr (z_{t}=j)$\ is the unconditional probability of regime $j$\ and $1_{z_{t}=j}=1$\ if $z_{t}=j$\ and $0$\ otherwise), and Assumptions (ref)-(ref) and (ref)-(ref) in Section (ref) hold conditioning on $\{z_{t},t=1,...,T\}$. Thus $\{z_{t},t=1,...,T\}$\ is allowed to be correlated with $f_{s}^{0}$\ and $e_{is}$\ for all $i$\ and $s$, and we do not need to know the true model of $\{z_{t},t=1,...,T\}$.
In vector form, the model can be written as:
where $\Lambda _{j}^{0}=(\lambda _{j1}^{0},...,\lambda _{jN}^{0})^{\prime }$ , $x_{t}=(x_{1t},...,x_{Nt})^{\prime }$ and $e_{t}=(e_{1t},...,e_{Nt})^{ \prime }$. Let $\Lambda ^{0}=(\Lambda _{1}^{0},...,\Lambda _{J^{0}}^{0})$ and let $E=(e_{1},...,e_{T})^{\prime }$ be the $T\times N$ matrix of errors. When there are no superscripts, $\Lambda _{j}$ and $\Lambda $ denote parameters as variables.
Since the factors are unobservable, regimes are defined in terms of the linear spaces spanned by the loadings. Two regimes are different if their loading spaces are different, and vice versa. More specifically, the identification condition is: for any $j$\ and $k$,
A sufficient condition for ((ref)) is:
Condition ((ref)) requires $\lim\limits_{N\rightarrow \infty }\frac{1}{N} (\Lambda _{j}^{0},\Lambda _{k}^{0})^{\prime }(\Lambda _{j}^{0},\Lambda _{k}^{0})\ $to be full rank for any $j$\ and $k$. Thus $\Lambda _{j}^{0}$\ and $\Lambda _{k}^{0}$\ are not allowed to share some columns, and columns of $\Lambda _{j}^{0}$\ could not be linear combination of $\Lambda _{k}^{0}$ \ and vice versa. An alternative sufficient condition for ((ref)) is:
where $g_{jk}$ is the eigenvector of $\lim\limits_{N\rightarrow \infty } \frac{1}{N}\Lambda _{k}^{0\prime }M_{\Lambda _{j}^{0}}\Lambda _{k}^{0}$ corresponding to nonzero eigenvalue. Condition ((ref)) only requires that the linear spaces spanned by $\Lambda _{j}^{0}$ and $\Lambda _{k}^{0}$ are different. Thus $\Lambda _{k}^{0}$\ and $\Lambda _{j}^{0}$\ are allowed to share some columns, and some columns of $\Lambda _{k}^{0}$ are allowed to be linear combinations of the columns of $\Lambda _{j}^{0}$ and vice versa, but $\Lambda _{k}^{0}$\ is not allowed to be a subset of $\Lambda _{j}^{0}$. For example, if there are two regimes with two factors in each regime and only the loadings of $f_{2t}$ (the second factor) switch across the regimes, then condition ((ref)) requires that $\min_{t}\left\vert f_{2t}\right\vert $\ is nonzero.
Note that condition ((ref)) does not rule out the possibility that any regime $j$ can be further decomposed into multiple regimes. Suppose the true model is $x_{t}=\Lambda _{j}^{0}f_{t}^{0}+e_{t}$\ if $z_{t}=j$, $j=1,2,3$, and $(\Lambda _{1}^{0},\Lambda _{2}^{0},\Lambda _{3}^{0})$\ satisfies condition ((ref)). If we consider $\Lambda _{1}^{0}$\ as the first regime and $(\Lambda _{2}^{0},\Lambda _{3}^{0})$\ as the second regime, the true model can be equivalently written as $x_{t}=\Lambda _{1}^{0}f_{t}^{0}+e_{t}$ \ if $z_{t}=1$, and $x_{t}=(\Lambda _{2}^{0},\Lambda _{3}^{0})f_{t}^{\ast }+e_{t}$\ if $z_{t}=2$\ $or$\ $3$, where $f_{t}^{\ast }=\left( f_{t}^{0\prime },0^{\prime }\right) ^{\prime }$\ if $z_{t}=2$\ and $ f_{t}^{\ast }=\left( 0^{\prime },f_{t}^{0\prime }\right) ^{\prime }$\ if $ z_{t}=3$. The equivalent model also satisfies condition ((ref)). However, while $plim\frac{1}{T}\sum\nolimits_{z_{t}=2\text{ }or\text{ }3}f_{t}^{\ast }f_{t}^{\ast \prime }$\ is positive definite, $plim\frac{1}{T} \sum\nolimits_{z_{t}=2}f_{t}^{\ast }f_{t}^{\ast \prime }$\ and $plim\frac{1}{ T}\sum\nolimits_{z_{t}=3}f_{t}^{\ast }f_{t}^{\ast \prime }$\ are not positive definite. To rule out the possibility that any regime $j$ can be further decomposed, we assume that
where $A_{j}$ denotes any subset of $\{t:z_{t}=j\}$ with cardinality $ \left\vert A_{j}\right\vert $ and $\lim \frac{\left\vert A_{j}\right\vert }{T }>0$. If $\frac{1}{\left\vert A_{j}\right\vert }\sum\nolimits_{t\in A_{j}}f_{t}^{0}f_{t}^{0\prime }$\ is not positive definite as $T\rightarrow \infty $\ for some $A_{j}$, then $A_{j}$\ and $\{t:z_{t}=j,t\notin A_{j}\}$ are considered as two separate regimes.
Consider the following log-likelihood function for Gaussian mixture in covariance:
where $\prod\nolimits_{t=1}^{T}L(x_{t}\left\vert z_{t};\Lambda ,\sigma ^{2}\right. )$ is the density of $(x_{1},...,x_{T})$ conditioning on $ (z_{1},...,z_{T})$, $\Pr (z_{1},...,z_{T})$ is the joint probability of $ (z_{1},...,z_{T})$,
$\Sigma _{j}$ is the covariance matrix of $x_{t}$ for regime $j$, and
The above log-likelihood function avoids estimating the factors. If the factors are estimated jointly with the loadings, we would not have the analytical first order conditions presented below, and consequently the EM algorithm would become infeasible.
Equation ((ref)) is a misspecified log-likelihood function. First, the state process $\{z_{t},t=1,...,T\}$\ is not specified yet, and the probability $\Pr (z_{1},...,z_{T})$\ depends on how we model the state process. Second, similar to the principal component estimator in Stock and Watson (2002) and Bai (2003), equation ((ref)) ignores the cross-sectional and serial dependence and heteroscedasticity of the error term. We may also take into account the heteroscedasticity as Doz, Giannone and Reichlin (2012) and Bai and Li (2012, 2016). With regime switching, the algorithm and the asymptotic analysis would be much more complicated, but the results should be conceptually similar.
Third, the factor dynamics are ignored. As shown in Bai (2003) for PCE and in Bai and Li (2012, 2016) for MLE, when there is no regime switching, the asymptotic properties of the estimated factors and the estimated loadings are robust to the presence of the factor dynamics if both $N$ and $T$ are large. We shall show in Section (ref) that when there is regime switching, the asymptotic results are also robust to the presence of the factor dynamics. More importantly, ignoring the factor dynamics greatly simplifies the computation algorithm for regime switching factor models. As shown below, with factor dynamics ignored, we just need to calculate principal components iteratively. If the factor dynamics are not ignored, Kim (1994)'s method would be very time consuming and may have convergence problems if $N$\ is large\footnote{ When there is no regime switching, as suggested by Doz et al. (2012), large $ N$ factor model with factor dynamics can be calculated by the EM algorithm. However, when there are both regime switching and factor dynamics, the EM algorithm also fails. This is because in the E-step we need to calculate the likelihood for each possible state chain $z_{1},....,z_{T}$ and there are $ (J^{0})^{T}$ possibilities, and in the M-step numerical optimization is still needed.}.
Fourth, equation ((ref)) implicitly assumes that $\mathbb{E}(f_{t}^{0})=0$ and $\mathbb{E}(f_{t}^{0}f_{t}^{0\prime })$ is stable within each regime, and $\mathbb{E}(f_{t}^{0}f_{t}^{0\prime })$ is absorbed into $\Lambda _{j}\Lambda _{j}^{\prime }$ in equation ((ref)). This does not matter, since all results of this paper still hold when $\mathbb{E}(f_{t}^{0})\neq 0$ and $\mathbb{E}(f_{t}^{0}f_{t}^{0\prime })$ is unstable within regime, as long as Assumption (ref) is satisfied.
The parameters $\Lambda $ and $\sigma ^{2}$ are estimated by maximizing $ l(\Lambda ,\sigma ^{2})$. Define $x_{1:t}\equiv (x_{1},...,x_{t})$ and $ z_{1:t}\equiv (z_{1},...,z_{t})$, and let $p_{tj\left\vert T\right. }\equiv \Pr (z_{t}=j\left\vert x_{1:T};\Lambda ,\sigma ^{2}\right. )$ denote the probability of $z_{t}=j$ conditional on $x_{1:T}$. Based on equation ((ref)), it can be easily verified that
where equation ((ref)) follows from
see Chapter 14.3 in Andersen (2003) for the details on calculating these derivatives. Set $\frac{\partial l(\Lambda ,\sigma ^{2})}{\partial \Lambda _{j}}$ to $0$, we have
$S_{j}$ can be considered as sample covariance matrix for $\Sigma _{j}$ based on importance sampling. The weights $p_{tj\left\vert T\right. }/\sum\nolimits_{t=1}^{T}p_{tj\left\vert T\right. }$ depend on the importance of the sample $x_{t}$ for regime $j$, the larger $p_{tj\left\vert T\right. }$ is, the more important $x_{t}$ is for regime $j$.
From equation ((ref)), we have $\Sigma _{j}\Lambda _{j}=\Lambda _{j}(\Lambda _{j}^{\prime }\Lambda _{j}+\sigma ^{2}I_{r_{j}^{0}})$. Left multiply $S_{j}\Sigma _{j}^{-1}$ on both sides, we have $S_{j}\Lambda _{j}=S_{j}\Sigma _{j}^{-1}\Lambda _{j}(\Lambda _{j}^{\prime }\Lambda _{j}+\sigma ^{2}I_{r_{j}^{0}})$. From equation ((ref)), we have $\Lambda _{j}=S_{j}\Sigma _{j}^{-1}\Lambda _{j}$, thus
If $\Lambda _{j}$ is a solution for equation ((ref)) and $\Lambda _{j}^{\ast }$ equals post-multiplying $\Lambda _{j}$ by the eigenvector matrix of $\Lambda _{j}^{\prime }\Lambda _{j}$, then $\Lambda _{j}^{\ast }$ is also a solution for equation ((ref)) and $\Lambda _{j}^{\ast \prime }\Lambda _{j}^{\ast }$ is diagonal. Thus we can directly choose the solution $\Lambda _{j}$ with $\Lambda _{j}^{\prime }\Lambda _{j}$ being diagonal. It follows that the solution $\Lambda _{j}$ is the eigenvectors of $S_{j}$ and $ \Lambda _{j}^{\prime }\Lambda _{j}+\sigma ^{2}I_{r_{j}^{0}}$ is the corresponding eigenvalues. We show in Appendix (ref) that $\sigma ^{2} $ satisfies the following condition:
Note that we do not need to specify the state process $\{z_{1},...,z_{T}\}$\ when deriving the first order conditions ((ref)) and ((ref)), and different models of $\{z_{1},...,z_{T}\}$\ correspond to different ways of calculating $p_{tj\left\vert T\right. }$. In the EM algorithm presented below, we consider $\{z_{1},...,z_{T}\}$\ as a Markov process regardless of what the true process of $\{z_{1},...,z_{T}\}$\ is.
Let $q^{0}=(q_{1}^{0},...,q_{J^{0}}^{0})^{\prime }$ denote the unconditional regime probabilities, $\phi ^{0}=(\phi _{1}^{0},...,\phi _{J^{0}}^{0})^{\prime }$ denote the initial probabilities of $z_{1}$, $Q^{0}$ denote the $(J^{0}\times J^{0})$ matrix of transition probabilities and $ Q_{jk}^{0}$ denote the probability of switching from state $k$ to state $j$. If there are no superscripts, $q$, $Q$ and $\phi $ denote parameters as variables.
For any given $Q$ and $\phi $, at the $h$-th iteration, let $\tilde{\Lambda} ^{(h)}$ denote the estimated loadings, $\tilde{\sigma}^{2(h)}$ denote the estimated variance, and $\Pr (z_{1},...,z_{T}\left\vert x_{1:T};\tilde{\theta }^{(h)}\right. )$ denote the probability of $z_{1:T}$ conditioning on $ x_{1:T}$ and evaluated at $\tilde{\theta}^{(h)}=(\tilde{\Lambda}^{(h)}, \tilde{\sigma}^{2(h)},Q,\phi )$. The EM algorithm maximizes the expectation of the log-likelihood of $(x_{1:T},z_{1:T})$ with respect to $\Pr (z_{1},...,z_{T}\left\vert x_{1:T};\tilde{\theta}^{(h)}\right. )$, i.e.,
Considering $z_{t}$ as a Markov process,\ $\Pr (z_{1},...,z_{T}\left\vert Q,\phi \right. )=\Pr (z_{1}\left\vert \phi \right. )\prod\nolimits_{t=2}^{T}\Pr (z_{t}\left\vert z_{t-1};Q\right. )$. Thus
where $\tilde{p}_{tjk\left\vert T\right. }^{(h)}=\Pr (z_{t}=j,z_{t-1}=k\left\vert x_{1:T};\tilde{\theta}^{(h)}\right. )$ and $ \tilde{p}_{tj\left\vert T\right. }^{(h)}=\Pr (z_{t}=j\left\vert x_{1:T}; \tilde{\theta}^{(h)}\right. )=\sum\nolimits_{k=1}^{J^{0}}\tilde{p} _{tjk\left\vert T\right. }^{(h)}$ are the smoothed probabilities based on $ x_{1:T}$ and $\tilde{\theta}^{(h)}$. Appendix (ref) presents a recursive algorithm for calculating $\tilde{p}_{tjk\left\vert T\right. }^{(h)}$. From equations ((ref)) and ((ref)), we have $\frac{\partial \log L(x_{t}\left\vert z_{t}=j;\Lambda _{j},\sigma ^{2}\right. )}{\partial \Lambda _{j}}=-\Sigma _{j}^{-1}\Lambda _{j}+\Sigma _{j}^{-1}x_{t}x_{t}^{\prime }\Sigma _{j}^{-1}\Lambda _{j}$. Thus
and it follows that
Similar to equation ((ref)), equation ((ref)) implies that
thus the columns of $\tilde{\Lambda}_{j}^{(h+1)}$ are the eigenvectors of $ \tilde{S}_{j}^{(h)}$ and the diagonal elements of $\tilde{\Lambda} _{j}^{(h+1)\prime }\tilde{\Lambda}_{j}^{(h+1)}+\tilde{\sigma} ^{2(h+1)}I_{r_{j}^{0}}$ are the corresponding eigenvalues. To save space, we show in Appendix (ref) that
Choose any $Q$\ and $\phi $\ such that $Q_{jk}>0$ \ for any $j$\ and $k$\ and $\phi _{k}>0$\textit{ \ for all }$k$. \textit{Start from randomly generated initial values of }$ \tilde{\Lambda}^{(0)}$ \textit{and }$\tilde{\sigma}^{2(0)}=1$\textit{. For }$ h=0,1,...,$
(E-step): calculate $\tilde{p}_{tjk\left\vert T\right. }^{(h)}$ \ using the algorithm in Appendix (ref), and calculate $ \tilde{S}_{j}^{(h)}=\sum\nolimits_{t=1}^{T}\tilde{p}_{tj\left\vert T\right. }^{(h)}x_{t}x_{t}^{\prime }/\sum\nolimits_{t=1}^{T}\tilde{p}_{tj\left\vert T\right. }^{(h)}$ with $\tilde{p}_{tj\left\vert T\right. }^{(h)}=\sum\nolimits_{k=1}^{J^{0}}\tilde{p}_{tjk\left\vert T\right. }^{(h)}$ ;
(M-step): given $\tilde{p}_{tjk\left\vert T\right. }^{(h)}$ \ and $\tilde{S}_{j}^{(h)}$, calculate $\tilde{\Lambda} _{j}^{(h+1)} $\ as the eigenvectors of $\tilde{S}_{j}^{(h)}$ \ corresponding to the $r_{j}^{0}$\ largest eigenvalues, and then normalize $\tilde{\Lambda}_{j}^{(h+1)}$\textit{\ such that }$\left\Vert \tilde{\Lambda}_{jl}^{(h+1)}\right\Vert ^{2}+\tilde{\sigma}^{2(h+1)}$\textit{ \ equals the }$l$\textit{-th largest eigenvalue of }$\tilde{S}_{j}^{(h)}$ \textit{\ for }$l=1,...,r_{j}^{0}$\textit{\ and equation ((ref)) is also satisfied, where }$\tilde{\Lambda}_{jl}^{(h+1)}$\textit{\ is the }$l$\textit{ -th column of }$\tilde{\Lambda}_{j}^{(h+1)}$\textit{. Note that the computation of }$\left\Vert \tilde{\Lambda}_{jl}^{(h+1)}\right\Vert ^{2}$ \textit{\ and }$\tilde{\sigma}^{2(h+1)}$\textit{\ requires iteration between equations ((ref)) and ((ref)).}
Iterate the E-step and the M-step until converge. Let $\tilde{ \Lambda}_{j}=(\tilde{\lambda}_{j1},...,\tilde{\lambda}_{jN})^{\prime }$ , $\tilde{\Lambda}=(\tilde{\Lambda}_{1},...,\tilde{\Lambda} _{J^{0}}) $\ and $\tilde{\sigma}^{2}$\ denote the estimated parameters, and let $\tilde{p}_{tj\left\vert T\right. }$\ and $ \tilde{p}_{tjk\left\vert T\right. }$\ denote the smoothed probabilities based on $x_{1:T}$\textit{\ and }$(\tilde{\Lambda},\tilde{ \sigma}^{2},Q,\phi )$\textit{. }
A special case of the above EM algorithm is when we choose $\phi =q$\ and $ Q=q1_{J^{0}}^{\prime }$\ ($1_{J^{0}}$\ denotes the $J^{0}\times 1$\ vector of ones), i.e., we consider $\{z_{1},...,z_{T}\}$\ as an independent process. For this case, the computation of $\tilde{p}_{tj\left\vert T\right. }^{(h)}$\ is simplified because the unsmoothed regime probabilities can be calculated directly by
This case is preferable if we knew the true process of $\{z_{1},...,z_{T}\}$ \ is independent. In practice, since the state process of the business cycle/stock market is highly persistent, smoothed regime probabilities that capture the persistence should perform significantly better, especially when mixed frequency data or ragged edge data (data released at non-synchronized dates) are used. The asymptotic results in Section (ref) and Section (ref) hold for any $Q$\ and $\phi $\ as long as $\phi _{j}>0$\ for any $j$ \ and $Q_{jk}>0$\ for any $j$\ and $k$, i.e., they hold for both the smoothed algorithm and the unsmoothed algorithm.
If the true process of $\{z_{1},...,z_{T}\}$ is Markovian, $Q_{jk}^{0}$ and $ \phi _{k}^{0}$ can be estimated by
We can also plug $\tilde{Q}_{jk}$ and $\tilde{\phi}_{k}$ back in the above EM algorithm and iterate between $(\tilde{\Lambda},\tilde{\sigma}^{2})$ and $ (\tilde{Q},\tilde{\phi})$ until convergence. This is the maximum likelihood estimator when $(Q,\phi )$ is estimated jointly with $(\Lambda ,\sigma ^{2})$ , see Appendix (ref) for details.
The asymptotic results in Section (ref) and Section (ref) also hold as long as $\tilde{\sigma}^{2}$ is bounded and bounded away from zero in probability. Consistency of $\tilde{\sigma}^{2}$ is not needed. We could restrict $\tilde{\sigma}^{2}$ in $[\frac{1}{C^{2}},C^{2}]$ for some large $C$ or simply fix down $\tilde{\sigma}^{2}=1$ to avoid the iteration between $ \tilde{\Lambda}_{j}^{(h+1)}$ and $\tilde{\sigma}^{2(h+1)}$. This only affects the Euclidean norm of the columns of $\tilde{\Lambda}_{j}^{(h+1)}$.
If the factor dynamics are taken into account, the expectation of $f_{t}$ conditioning on $x_{1:t}$ is
which is formidable since we need to calculate $\Pr (z_{1:t}\left\vert x_{1:t};\tilde{\Lambda},\tilde{\sigma}^{2},Q,\phi \right. )$ for each possible $z_{1:t}$, i.e., we need to calculate $(J^{0})^{t}$ probabilities. For large $N$, the benefit of considering factor dynamics is marginal and outweighed by the computational simplicity of ignoring factor dynamics. If the factor dynamics are ignored, the expectation of $f_{t}$ conditioning on $ x_{1:T}$ is
Note that the dimension of $\tilde{\Lambda}_{j}^{\prime }(\tilde{\Lambda}_{j} \tilde{\Lambda}_{j}^{\prime }+\tilde{\sigma}^{2}I_{N})^{-1}x_{t}$\ is different across $j$\ if $r_{j}^{0}$\ is different across $j$. Here and also in the proof of Theorem (ref), when we add two vectors of different dimensions, we implicitly augment the vector of smaller dimension with zeros to make the dimensions of these two vectors equal. Thus $\tilde{f}_{t}$\ is a $\max r_{j}^{0}$\ dimensional vector.
We assume the following conditions hold as $(N,T)\rightarrow \infty $. These conditions are mainly Assumptions A-G in Bai (2003) adapted to the current regime switching setup.
Assumption (ref) corresponds to Assumption A in Bai (2003). Assumption (ref)(1) rules out the possibility that for regime $j$, the subsample $\{t:z_{t}=j\}$ can be further decomposed into multiple regimes, see the discussion in Section (ref). The factor process is allowed to be dynamic such that $C(L)f_{t}=\epsilon _{t}$. Assumption (ref)(2) assumes that the factors have bounded moments.
Assumption (ref)(1) corresponds to Assumption B in Bai (2003). Assumption (ref)(1) ensures that each factor has a nontrivial contribution within each regime, and $\left\Vert \lambda _{ji}^{0}\right\Vert $ is assumed to be uniformly bounded over $i$. Assumption (ref)(2) is the identification condition for determining which regime each $x_{t}$ belongs to, see Section (ref) for details on the implication of this condition.
Assumption (ref) is modified slightly from Assumption C in Bai (2003). The error term is allowed to have limited cross-sectional and serial dependence as well as heteroscedasticity.
The asymptotic results in Section (ref) and Section (ref) are valid as long as Assumption (ref) holds and the other assumptions in this section hold conditioning on $\{z_{t},t=1,...,T\}$. Thus the state process $\{z_{t},t=1,...,T\}$\ is allowed to be non-Markovian and correlated with $f_{s}^{0}$\ and $e_{is}$\ for all $i$\ and $s$. Knowledge of the true state process is not needed.
Assumption (ref) is modified slightly from Assumption D in Bai (2003). Assumption (ref)(1) assumes that the errors are weakly correlated across $i$ for each $t$. When $\beta =2$, Assumption (ref) (1) is implied by Assumptions (ref)(1), (ref)(1) and (ref)(2). Assumption (ref)(2) assumes that the errors are weakly correlated across $t$ for each $i$. Assumption (ref)(2) is implied by Assumptions (ref)(2), (ref)(1) and (ref)(4) if we further assume the factors are nonrandom or independent with the errors.
Assumption (ref) corresponds to Assumption G in Bai (2003). With Assumption (ref), the loadings and the factors are identifiable up to a rotation. For identification of the loading space and the factor space, Assumption (ref) is not needed.
Assumption (ref) corresponds to Assumption F in Bai (2003). Part (3) and part (4) are just central limit theorems and will be used for deriving the limit distributions of the estimated factors and loadings.
Theorem (ref) shows that the estimated loading space is consistent without observing the state variable $z_{t}$. Note that the estimated loadings $\tilde{\Lambda}_{j}$\ and the estimated regime probabilities $ \tilde{p}_{tj\left\vert T\right. }$\ depend on each other, but the standard technique in Bai (2003) for analyzing $\tilde{\Lambda}_{j}$\ is applicable only when $\tilde{p}_{tj\left\vert T\right. }=1_{z_{t}=j}$. This is the first technical difficulty we encounter in going from one regime to multiple regimes. The crucial point for Theorem (ref) is that if the linear spaces spanned by $\tilde{\Lambda}_{j}$\ and $\Lambda _{j}^{0}$\ differ too much, as long as $\min \phi _{j}>0$\ and $\min Q_{jk}>0$, the likelihood of $ \tilde{\Lambda}$\ would be smaller than the likelihood of $\Lambda ^{0}$\ uniformly over all possible $\{z_{1},...,z_{T}\}$, i.e.,
This crucial point is due to large $N$, see the Appendix for the formal proof. Based on Theorem (ref), we show that the estimated regime probabilities are consistent.
Note that $\eta $ could be large but it is fixed as $(N,T)\rightarrow \infty $. $\alpha $ and $\beta $ could also be large as long as Assumptions (ref)(2), (ref)(1) and (ref)(1) are satisfied. Theorem (ref) shows that $\tilde{p}_{tj\left\vert T\right. }$ is consistent as $ N\rightarrow \infty $ and is uniformly consistent if $T$ is relatively small compared to $N$. The proof utilizes the exponential likelihood ratio.
Theorem (ref) implies that we can consistently identify which regime $ x_{t}$ belongs to for all $t$, if there is common regime switching in the loadings and the dimension of $x_{t}$ tends to infinity. Theorem (ref) also implies that we can consistently detect regime switching right after the turning point with only one observation, so that we do not need to wait for many observations of the time series from the new regime. This could improve the speed of detection of new turning points, especially when high frequency data is used.
An interesting special case is when the proposed algorithm is applied to factor models with common breaks in the loadings. Various methods are proposed recently for estimating the break points, Theorem (ref) implies that we can also consistently estimate the break points using the proposed EM algorithm.
If the true states $z_{t}$ were known, asymptotic properties of the estimated loadings and factors are straightforward. Based on Theorem (ref), we shall show that using estimated regime probabilities does not affect the asymptotic results. Define $W_{jNT}=\frac{1}{N}(\tilde{\Lambda} _{j}^{\prime }\tilde{\Lambda}_{j}+\tilde{\sigma}^{2}I_{r_{j}^{0}})(\frac{1}{T }\sum\nolimits_{t=1}^{T}\tilde{p}_{tj\left\vert T\right. })$ and $H_{j}= \frac{\sum\nolimits_{t=1}^{T}f_{t}^{0}f_{t}^{0\prime }1_{z_{t}=j}}{T}\frac{ \Lambda _{j}^{0\prime }\tilde{\Lambda}_{j}}{N}W_{jNT}^{-1}$, then we have:
Proposition (ref) is an important auxiliary result, and part (1) and part (2) corresponds to Lemma A.3 and Proposition 1 in Bai (2003), respectively. Lemma A.3 in Bai (2003) is based on the fact that the estimated factors are $\sqrt{T}$\ times the eigenvectors corresponding to the $r$\ largest eigenvalues\footnote{$r$ denotes the number of factors in Bai (2003).} of $XX^{\prime }$\ and consequently $\tilde{\Lambda}^{\prime } \tilde{\Lambda}$\ is a diagonal matrix consisting of the $r$\ largest eigenvalues of $\frac{1}{T}\sum\nolimits_{t=1}^{T}x_{t}x_{t}^{\prime }$. However, here the first order condition ((ref))\ only tells us the columns of $\tilde{\Lambda}_{j}$\ are the eigenvectors of $S_{j}$\ and $ \tilde{\Lambda}_{j}^{\prime }\tilde{\Lambda}_{j}+\tilde{\sigma} ^{2}I_{r_{j}^{0}}$\ are the corresponding eigenvalues. Condition ((ref) ) does not tells us whether these eigenvalues are the $r_{j}^{0}$\ largest eigenvalues of $S_{j}$\ or not.\ This is the second technical difficulty we encounter in going from one regime to multiple regimes. Our proof strategy of Proposition (ref) utilizes Theorem (ref) and is totally different from Bai (2003)'s proof for his Proposition 1.
Theorem (ref) establishes the convergence rate of the estimated loading space for each regime. This could help us study the effect of using estimated loadings on subsequent applications. For example, if the estimated loadings are used to construct portfolios, Theorem (ref) could help us calculate how the estimation error contained in $\tilde{\Lambda}_{j}$ would affect the performance of these portfolios.
Theorem (ref) shows that for each $j$ and $i$, $\tilde{\lambda}_{ji}$ has a limiting normal distribution. This allows us to construct confidence interval for the estimated loadings. Also note that the rotation matrix $ H_{j}$ is different for different regime.
Theorem (ref)(2) shows that the limit distribution of $\tilde{f}_{t}$ is mixed normal, since the rotation matrix $H_{z_{t}}^{-1}$ and the asymptotic variance depend on the state variable $z_{t}$. Theorem (ref)(1) establishes the convergence rate of the estimated factor space. Note that if $\{\tilde{f}_{t},t=1,...T\}$ is used as proxies for the true factors in factor-augmented forecasting (or factor-augmented VAR), the forecasting equation (or the VAR equation) would have induced regime switching in the model parameters, because $H_{z_{t}}^{-1}$ depends on $ z_{t} $. For illustration, consider the following $h$-period ahead forecasting model using factors and some other observable variables $W_{t}$: $y_{t+h}=a^{\prime }f_{t}^{0}+b^{\prime }W_{t}+u_{t+h}$. If $\tilde{f}_{t}$ is used as proxies for $f_{t}^{0}$, the model can be written as
The first term on the right hand side is asymptotically negligible. It is easy to see that the coefficient $a^{\prime }H_{z_{t}}$ depends on $z_{t}$ and this need to be taken into account when we estimate the forecasting equation. Finally, we show that the estimated transition probability matrix is also consistent when $\{z_{1},...,z_{T}\}$\ is a Markov process.
Given the number of factors $(r_{1},...,r_{J})$ and the number of regimes $J$ , let $(\tilde{\Lambda}_{1,r_{1}},...,\tilde{\Lambda}_{J,r_{J}})$ be the solution for maximizing the log-likelihood $l(\Lambda _{1,r_{1}},...,\Lambda _{J,r_{J}},\sigma ^{2},Q,\phi )$. Here we use $\Lambda _{j,r_{j}}$ to emphasize that $\Lambda _{j,r_{j}}$ is of dimension $N\times r_{j}$. The criterion we propose for model selection is:
where $g(N,T)$ is a penalty function depending on both $N$ and $T$, and $ b(\cdot )$ is a positive and decreasing function with $b(1)=1$, e.g., $ b(r_{j})=\frac{1}{r_{j}}$. For each $J$, the numbers of factors are estimated by
and then the number of regimes is estimated by
where $\bar{r}$ is the maximal number of factors in each regime and $\bar{J}$ is the maximal number of regimes. In the following theorem we show that $( \tilde{r}_{1},...,\tilde{r}_{J})$ and $\tilde{J}$ are consistent.
Note that the condition $\lim\limits_{N\rightarrow \infty }\frac{1}{N} \Lambda _{k}^{0\prime }M_{\Lambda _{j}^{0}}\Lambda _{k}^{0}\neq 0$ allows $ \Lambda _{k}^{0}$ and $\Lambda _{j}^{0}$ to share some columns, i.e., Theorem (ref) holds for the case where the loadings of some (but not all) factors remain the same across different regimes.
The basic idea behind Theorem (ref) is similar to Theorem 2 of Bai and Ng (2002), i.e., add a penalty term that converges to zero but slowly enough so that underparameterized models and overparameterized models will not be chosen. Here the penalty $(g(N,T))^{b(r_{j})}$ converges to zero because $ g(N,T)\rightarrow 0$ and $b(r_{j})$ is positive, and $\delta _{NT}(g(N,T))^{b(r_{j})}\rightarrow \infty $ because $\delta _{NT}g(N,T)\rightarrow \infty $ and $b(r_{j})$ is a decreasing function of $ r_{j}$ with $b(1)=1$. Compared to Bai and Ng (2002), the difference and difficulty here is that the number of regimes is unknown and the number of factors in each regime may be different. For example, suppose the true model is $(r_{1}=2,r_{2}=1,J=2)$ and the two columns in $\Lambda _{1}^{0}$ are linearly independent with $\Lambda _{2}^{0}$. This model can be equivalently written as
i.e., there is only one regime and there are three factors in this regime. The difference between the log-likelihood of the true model $ (r_{1}=2,r_{2}=1,J=2)$ and the log-likelihood of the equivalent model $ (r_{1}=3,J=1)$ is negligible and clearly Bai and Ng (2002) is not applicable to this example.
Our solution is to add a penalty term for each regime and let the penalty term of different regime have different asymptotic order, so that overestimating the number of factors in one regime can not be compensated by underestimating the number of factors in another regime. For example, the penalty for the equivalent model is $(g(N,T))^{b(3)}$ while the penalty for the true model is $(g(N,T))^{b(2)}+(g(N,T))^{b(1)}$. Since $\frac{ (g(N,T))^{b(3)}}{(g(N,T))^{b(2)}+(g(N,T))^{b(1)}}\rightarrow \infty $ as $ (N,T)\rightarrow \infty $, the true model would be chosen with probability approaching one as $(N,T)\rightarrow \infty $. The formal proof of Theorem (ref) is provided in the Appendix.
Our method can also be used to consistently determine the number of factors and the number of breaks for factor models with multiple common breaks in the loadings. If we replace $l(\tilde{\Lambda}_{1,r_{1}},...,\tilde{\Lambda} _{J,r_{J}},\sigma ^{2},Q,\phi )$ in expression ((ref)) by minus the minimum of the least squares over all possible break points and calculate $( \tilde{r}_{1},...,\tilde{r}_{J},\tilde{J})$ as expressions ((ref))-((ref)), then it is not difficult to prove that we still have $\Pr (\tilde{J} =J^{0}$ and $\tilde{r}_{j}=r_{j}^{0}$ for all $j)\rightarrow 1$ as $ (N,T)\rightarrow \infty $. As we discussed in the Introduction, recently the literature on the factor loading instability issues developed quite a lot, but as far as we know, there are very few (if any) consistent model selection procedures that allow $r_{j}$ to be different across $j$ and allow $\Lambda _{k}^{0}$ and $\Lambda _{j}^{0}$ to share some columns.
In this section, we perform simulations to confirm the theoretical results and examine the finite sample performance of our methods under various empirically relevant scenarios.
The data is generated as follows:
i.e., we consider two regimes. For the factors and the loadings, we consider four data generating processes (DGP) as listed below:
DGP 1: There are two factors in both regimes and the loadings of both factors have regime switching.
DGP 2: There are two factors in both regimes and only the loadings of the second factor have regime switching.
DGP 3: There is one factor in both regimes and its loadings have regime switching.
For DGP1 - DGP3, the factors are generated as follows:
$\epsilon _{t,p}$\ is i.i.d. $N(0,1)$, and $f_{1,p}^{0}$\ is i.i.d. $N(0, \frac{1}{1-\rho ^{2}})$ so that the distributions of the factors are stationary. Serial correlation of the factors is controlled by the scalar $ \rho $.
DGP 4: The loadings are generated in the same way as DGP2. $ f_{t,1}^{0}$\ is generated as i.i.d. $N(0,1)$\ and $ f_{t,2}^{0}$\ is generated as uniform $(0.5,1.5)$.
The errors are generated as follows:
where $v_{t}=(v_{1,t},...,v_{N,t})^{\prime }$\ is i.i.d. $N(0,\Omega )$ for $ t=2,...,T$ and $(e_{1,1},...,e_{N,1})^{\prime }$\ is $N(0,\frac{1}{1-\zeta ^{2}}\Omega )$\ so that the distributions of the errors are stationary. Serial correlation of the errors is controlled by the scalar $\zeta $. For $ \Omega $, we set $\Omega _{ij}=\xi ^{\left\vert i-j\right\vert }$ for some $ \xi $ between 0 and 1, thus cross-sectional dependence of the errors is controlled by $\xi $. In addition, the processes $\{\epsilon _{t,p}\}$\ and $ \{v_{it}\}$\ are mutually independent for all $p$ and $i$.
The loadings are generated as follows: For DGP1, both $\lambda _{1i}^{0}$ and $\lambda _{2i}^{0}$\ are generated as i.i.d. $N(0,\frac{1-\rho ^{2}}{ 1-\zeta ^{2}}\frac{2R^{2}}{1-R^{2}}I_{2})$ across $i$, and $\lambda _{1i}^{0} $ and $\lambda _{2i}^{0}$ are also independent with each other. For DGP2, $\lambda _{1i}^{0}$ and the second element of $\lambda _{2i}^{0}$ are generated as i.i.d. $N(0,\frac{1-\rho ^{2}}{1-\zeta ^{2}}\frac{2R^{2}}{ 1-R^{2}}I_{3})$ across $i$. For DGP3, both $\lambda _{1i}^{0}$ and $\lambda _{2i}^{0}$\ are generated as i.i.d. $N(0,\frac{1-\rho ^{2}}{1-\zeta ^{2}} \frac{R^{2}}{1-R^{2}})$ across $i$, and $\lambda _{1i}^{0}$ and $\lambda _{2i}^{0}$ are also independent with each other. All loadings are independent of the factors and the errors. The variance $\frac{1-\rho ^{2}}{ 1-\zeta ^{2}}\frac{2R^{2}}{1-R^{2}}$ guarantees that the regression R-square\ of each series $i$\ is equal to $R^{2}$, this controls the signal-noise ratio. Following the literature, we set $R^{2}=0.5$.
For the state process $\{z_{t},t=1,...,T\}$, we consider four cases as listed below:
Regime Pattern 1: US business cycle 1945Q2-2020Q1
Regime Pattern 2: single common break at $t=T/2$
Regime Pattern 3: two common breaks at $t=T/3$ and $t=2T/3$, and the loadings switch back after the second break
Regime Pattern 4: a randomly generated Markov process
Regime pattern 1 is based on the US business cycle from 1945 Quarter 2 to 2020 Quarter 1, as determined by the NBER business cycle dating committee. There are 75 years (300 quarters) in total, thus we have $T=300$. For $ t=1,...,300$, $z_{t}=1$ if the US economy at time $t$ is in expansion and $ z_{t}=2$ if the US economy at time $t$ is in recession. The transition probabilities of the state process calibrated to the US business cycle is $ Q_{11}^{0}=0.95$ and $Q_{22}^{0}=0.72$ (average duration of expansion is $ 1/(1-Q_{11}^{0})=20$ and average duration of recession is $ 1/(1-Q_{22}^{0})\approx 3.5$).
Regime patterns 2 and 3 correspond to the case where loadings have single common break\ and multiple common breaks, respectively. Regime patterns 3 is especially interesting since the case where there are multiple breaks and the loadings switch back to their original values after the second break is rarely studied in the literature. Various methods are proposed in the literature recently for estimating the break points, here we perform simulations for regime patterns 2 and 3 to evaluate the finite sample performance of our method when it is applied to these interesting cases.
Regime pattern 4 is a Markov process randomly generated with transition probabilities $Q_{11}^{0}=0.95$ and $Q_{22}^{0}=0.72$, and $ \{z_{t},t=1,...,T\}$ is independent with $f_{s}^{0}$ and $e_{is}$ for all $i$ and $s$. Regime patterns 1-3 are prespecified and are not necessarily Markov processes, thus here we consider regime pattern 4 to evaluate the performance of our method when applied to a Markov state process.
We study both the unsmoothed algorithm and the smoothed algorithm.\ The key difference is that in the E-step, the former uses unsmoothed regime probabilities while the latter uses smoothed regime probabilities. Both algorithms start from randomly generated initial values of the loadings and iterate between the E-step and the M-step until convergence. To search for the global maximum of the likelihood function, we generate initial values randomly for many times and take the one with the largest likelihood. For other parameters, we set $\sigma ^{2}=1$, $q_{j}=0.5$\ for $j=1,2$, $\phi _{k}=0.5$ for $k=1,2$, $Q_{11}=0.95$ and $Q_{22}=0.72$. $Q_{11}$ and $Q_{22}$ are calibrated to regime pattern 1. Once we get the estimated regime probabilities and the estimated loadings, $\tilde{Q}_{11}$ and $\tilde{Q} _{22}$ are estimated by equation ((ref)), and the factors are estimated by equation ((ref)).
Figure 1 displays the smoothed probabilities of regime 2 for DGP 1 with $ (N,T)=(100,300)$ and $(\rho ,\zeta ,\xi )=(0,0,0)$. Subfigures 1-4 of Figure 1 correspond to regime patterns 1-4, respectively. It is easy to see that in all subfigures when the true regime is regime 1, the smoothed probabilities stay at zero with only a few short and mild spikes. At the beginning of each shaded region, the smoothed probabilities increase to one instantly, and at the end of each shaded region, the smoothed probabilities instantly decrease to zero. Figure 2 displays the unsmoothed probabilities of regime 2 for DGP1 under the four regime patterns with $(N,T)=(100,300)$ and $(\rho ,\zeta ,\xi )=(0,0,0)$. The estimated probabilities still stay at zero when it's regime 1 and instantly increase to one (decrease to zero) when there is regime switching, but compared to Figure 1, Figure 2 shows more and sharper spikes (upward or downward). These spikes are false positives in detecting regime switching. Figure 3 and Figure 4 display the smoothed and the unsmoothed probabilities of regime 2 for DGP2, respectively. The performance of the estimated probabilities deteriorates since for DGP2 only one factor has regime switching in its loadings. Overall, Figures 1-4 confirm the theoretical results that turning points (break points) can be identified consistently if $N$ is large.
Comparing Figure 2 to Figure 1 and Figure 4 to Figure 3, it is obvious that the smoothed probabilities performs much better than the unsmoothed probabilities. Many false positives in Figure 3 and Figure 4 are eliminated by the smoother. This is because for each $t$, regimes at $t-1$ and $t+1$ contains information for detecting the regime at period $t$. Comparing subfigures 2-3 to subfigures 1 and 4 in Figures 1-4, we can see that the performance of the estimated probabilities under regime patterns 2-3 is better than the performance under patterns 1 and 4. This is also because regimes at the neighborhood periods provide information for the current regime. Roughly speaking, the performance is better when the regime pattern is relatively simple. In addition, we can also see that the performance under regime pattern 1 is slightly better than the performance under regime pattern 4. This is because the subsample size of regime 2 under pattern 4 is larger than the subsample size under pattern 1 (72 vs 45). In general, we find that to guarantee good performance, the subsample size for each regime should be not less than 40.
Figure 5 focuses on regime pattern 1 and displays the estimated probabilities of regime 2 for DGP1 and DGP2 with $N=200$. Comparing the subfigure 3 and subfigure 4 of Figure 5 to subfigure 1 of Figure 3 and subfigure 1 of Figure 4, it is easy to see that $N=200$ improves the performance of the estimated probabilities. Figure 6 also focuses on regime pattern 1 and displays the smoothed probabilities of regime 2 for DGP1 and DGP2 with $(\rho ,\zeta ,\xi )=(0.5,0,0)$ or $(0,0.5,0.5)$. Comparing to subfigure 1 of Figure 1 and subfigure 1 of Figure 3, it seems that the value of $(\rho ,\zeta ,\xi )$ does not affect the performance too much if they were far away from 1.
Figure 7 displays the smoothed and unsmoothed probabilities for regime pattern 4 (a randomly generated Markov process) under DGP1 with $ (N,T)=(50,500)$\ or $(N,T)=(500,50)$. The smoothed probabilities still perform well even under such extreme case, especially when $T=50$, the subsample size of regime 2 is just 13. However, the unsmoothed probabilities in subfigures 3-4 deteriorate obviously, compared to subfigure 4 of Figure 2.
Figure 8 displays the smoothed and unsmoothed probabilities under DGP4 for regime patterns 1 and 2 with $(N,T)=(100,300)$. DGP4 modifies DGP2 so that the second factor stays away from zero. Comparing subfigures 1-2 to subfigures 1-2 of Figure 3 and subfigures 3-4 to subfigures 1-2 of Figure 4, we can see that the performance improvement is quite significant. This is because under DGP2, the second factor $f_{t,2}^{0}$ is likely to be close to zero and it is difficult to identify the regime of $x_{t}$\ when $ f_{t,2}^{0} $\ is close to zero. Thus Figures 3-4 reflect more of the identification problem when the factors equal zeros.
Finally, to access the adequacy of the asymptotic distributions of the estimated loadings and factors in approximating their finite sample counterparts, we display in Figures 9-12 the histograms of the standardized estimated factors for $t=T/2$ and the standardized estimated loadings for $ i=N/2$ under DGP3. The number of simulations is 1000. The histograms are normalized to be a density function and the standard normal density curve is overlaid on them for comparison. It is easy to see that in all subfigures of Figures 9-12, the standard normal density curve provides good approximation to the normalized histograms. The histograms of the estimated factors in Figure 9 are slightly fat-tailed because of bad initial values. Comparing the four rows in each of Figures 9-12, we can see that the estimated loadings and factors using the smoothed algorithm perform better than using the unsmoothed algorithm, $(\rho ,\zeta ,\xi )=(0.5,0.5,0.5)$ does not matter too much, and $N=200$ significantly improves the performance.
The number of initial value trials also significantly affect the performance. We find that for regime pattern 1, normally 5 trials are enough, but to guarantee good performance in all of 1000 replications, 30 trials are needed. For regime pattern 4, normally 2 trials are enough and 15 trials are needed to guarantee good performance in all replications. For regime patterns 2-3, 5 trials are enough for all replications. In general, more trials are needed when the regime pattern is complex and the subsample size is small.
In addition, we also present in Table (ref) the average $R^{2}$ of the estimated loadings of regime 1 and regime 2 projecting on the true loadings, the average $R^{2}$ of the estimated factors, and the average absolute error of the estimated transition probabilities. It is easy to see that in Table 1, $R_{l1}^{2}$ and $R_{l2}^{2}$ are always close to one. $R_{Hf}^{2}$ is always close to one but $R_{f}^{2}$ is much smaller than $R_{Hf}^{2}$. This is because $R_{Hf}^{2}$ considers the regime specific rotation matrix, as shown in Theorem (ref)(1). In summary, results in Figures 1-12 and Table 1 lend strong support to the theoretical results and illustrate the usefulness of the proposed EM algorithms.
In this section we apply the proposed method to detect turning points of US business cycle from 02/1980 to 01/2023 in real-time using the FRED-MD (Federal Reserve Economic Data - Monthly Data) data set. The FRED database is maintained by the Research division of the Federal Reserve Bank of St. Louis, and is publicly accessible and updated in real-time. The 02/2023 vintage of the FRED-MD data set contains 128 unbalanced monthly time series from 01/1959 to 01/2023, including eight groups (output and income, labor market, housing, consumption and inventories, money and credit, prices, stock market). After removing those series with missing values and data transformation\footnote{ See the Appendix of McCracken and Ng (2016) for the details of data description and transformation.}, we have 106 balanced monthly series ranging from 03/1959 to 01/2023. Finally, the data is demeaned and standardized.
For each month from 02/1980 to 01/2023 (516 months in total), we use the data from 03/1959 to that month for calculating the probability of recession of that month, i.e., we behave as if we were standing at that month\footnote{ For simplicity, we do not use the vintage data of that month. Compared to the vintage data, the data we use contains revision in some series if more accurate observations were available after that month, but previous studies on business cycle dating show that data revisions have little effects on the results.}. More specifically, we apply the EM algorithm in Section (ref) to the data from 03/1959 to the previous month to estimate the model parameters\footnote{ US business cycle from 03/1959 to the previous month as determined by NBER is used as the initial values for probabilities.}, and then use the estimated parameters and the data from 03/1959 to that month to calculate the filtered probability of recession for that month. Since the data of that month is available at the end of that month or the beginning of the next month, new recession or expansion starting from the beginning of that month could only be detected with at least one month delay. \
To convert the recession probability of each month into a binary variable that indicates the state of the economy in that month, we compare the estimated recession probability to a prespecified threshold. More specifically, if the previous turning point is a trough and the recession probability of month $t$ exceeds 0.8 for the first time after the previous turning point, month $t$ would be considered as a new turning point from expansion to recession. Similarly, if the previous turning point is a peak and the recession probability of month $t$ falls below 0.2 for the first time after the previous turning point, month $t$ would be considered as a new turning point from recession to expansion. For robustness check, we also consider $(0.9,0.1)$ as the threshold, the results are quite similar.
We consider the turning points determined by the NBER BCDC (business cycle dating committee) as the benchmark for comparison and we mainly focus on the accuracy and speed of the proposed method in detecting turning points. The proposed method is applied to both the whole panel and a subset of the whole panel which consists of only the first 50 series among all 106 series. The results of using only the first 50 series are better. We conjecture that this is mainly because not all 106 series had regime switching in the factor loadings at each turning point determined by the NBER BCDC\footnote{ The NBER BCDC mainly focuses on four series, (1) non-farm payroll employment, (2) industrial production, (3) real manufacturing and trade sales, and (4) real personal income excluding transfer payments.}, or some series had regime switching in their loadings at time periods that are different from the NBER BCDC turning points. Thus we may further improve the performance of the proposed method by selecting series that are most relevant to and synchronous with the business cycle. A careful selection is out of the scope of this paper.
Table (ref) presents the real-time results of 02/1980-02/2020 using the first 50 series. The number of factors in each regime is set to be six. mm/yyyy in the second and the seventh row indicate the starting month of each recession and expansion. The row corresponds to "NBER BCDC", "Chauvet Piger" and "This paper" shows the number of months it takes the NBER BCDC, Chauvet and Piger (2008) and this paper to detect each recession and expansion, respectively. For example, the recession starting from the beginning of February 1980 would be detected by the NBER BCDC at the beginning of June 1980, by Chauvet and Piger (2008) at the beginning of August 1980, and by this paper at the beginning of May 1980, respectively. Overall, it is easy to see that this paper detects turning points much faster than NBER BCDC and slightly faster than Chauvet and Piger (2008). On average, this paper detects recessions with 6.25 months delay and expansions with 5.4 months delay, NBER BCDC detects recessions with 7.4 months delay and expansions with 14.8 months delay, and Chauvet and Piger (2008) detects recessions with 8.6 months delay and expansions with 6.2 months delay.
We also detect two recessions after the Covid-19 pandemic, one from 03/2020 to 08/2020 and the second from 02/2021 to 05/2021, so we would have detected the 03/2020-08/2020 recession in 04/2020 because the data for 03/2020 is available with one month delay. This is quite interesting, given that our method detects recessions with 6.25 months delay on average during the period 02/1980-02/2020.
Table (ref) shows that using more series could improve the speed of turning points detection. However, using more series could also bring in false positives (turning points detected by the proposed method using many series but not detected by NBER BCDC), because the extra series may not be synchronous with the NBER business cycle. Here we detect eight false recessions: 09/1983-11/1983, 10/1986-02/1987, 07/1989-10/1989, 01/1993-02/1993, 01/1995-03/1995, 08/1998, 05/2000-08/2000, 06/2010-10/2010, and one false expansion: 02/1982. While these false positives should not be ignored, most of them only last for a very short periods and would have little effect on macroeconomic policy. Overall, our results illustrate the potential of using a large number of series and factor models with common loading switching for real-time detection of the business cycle turning points.
To get rid of those false recessions, one possible solution is to select time series that are synchronous with the NBER business cycle, and another promising solution is to extend our results to the case where regime switching in the loadings is approximately synchronous rather than exactly synchronous. In fact, Stock and Watson (2014) mainly focuses on how to combine different turning points of many individual series (determined by the Bry-Boschan algorithm) into a single point. \
The exposure of economic time series to common factors may switch depending on state variables such as fiscal policy, monetary policy, business cycle stage, stock market volatility, technology and so on. For consistent estimation of the factor structure, it is crucial to take into account such regime switching phenomena. This paper considers maximum likelihood estimation for large factor models with common regime switching in the loadings and proposes EM algorithm for computation, which is easy to implement and runs fast even when $N$ is large. Convergence rates and limit distributions of the estimated loadings and the estimated factors are established under the approximate factor model setup. This paper also shows that when $N$ is large, regime switching can be identified consistently and only one observation after the switching point is needed. This allows us to detect regime switching at very early times. Monte Carlo simulations confirm the theoretical results and good performance of our method. An application to the FRED-MD dataset demonstrates the potential of using many time series with our method for detection of the business cycle turning points.
Some related topics are worth further study. First, it would be interesting to see the performance of the portfolio constructed using regime specific loadings, and how the identified regime is related to exogenous variables such as market volatility and money growth. Second, our results imply that the forecasting equation would have induced regime switching if the estimated factors are used for forecasting, so we want to know whether it indeed matters. Finally, a selection of time series that are most synchronous with or related to business cycle could improve the speed and accuracy of our method for turning points detection, so we would like to see how much we can achieve after careful selection.
We are deeply indebted to the Editor, Serena Ng, the Associate Editor, and two referees for very useful comments and suggestions which have helped to improve and develop the paper further. We wish to thank Daniele Massacci for some useful discussion on a previous version of the paper. The usual disclaimer applies. We acknowledge financial support from the Centre for Econometric Analysis, Bayes Business School, London (UK), and the Fundamental Research Funds for the Central Universities, Peking University, Beijing (China).
\setcounter{page}{1}