Modelling Large Dimensional Datasets with Markov Switching Factor Models

Modelling Large Dimensional Datasets with Markov Switching Factor Models

\thispagestyle{empty}

\footnotetext[1]{Universit\`a di Bologna, Department of Economics, \url{[email removed]}.}

\footnotetext[2]{King's College London, King's Business School, \url{[email removed]}. Corresponding author.

The paper greatly benefited from comments from conference participants at the 10th Italian Congress of Econometrics and Empirical Economics, the 28th International Panel Data Conference, the EEA-ESEM Barcelona 2023, the 2023 NBER-NSF Time Series Conference, and the 16th Annual SoFiE Meeting. Any errors and omissions are the authors' own responsibility only.}

Introduction

This paper develops a comprehensive approach for the analysis of large dimensional models exhibiting an approximate factor structure, in which the loadings are subject to regime shifts driven by a first order latent Markov process. We label these large dimensional Markov Switching factor models.

Since the works of hamilton1989new, and diebold78measuring, and inspired by the seminal paper of GQ1973, Markov switching models have been widely used in the empirical analysis of macroeconomic and financial time series data: Hamilton_2016_Chapter gives an overview from a macroeconomic perspective, and Doz_Ferrara_Pionnier_2020_WP present recent evidence of their usefulness for turning-point detection and macroeconomic forecasting; Guidolin_2011_Chapter, and AT2012, provide a comprehensive survey in relation to financial markets; see also Qu_Zhuo_2021_ReStat and references therein for more recent advances. However, to the very best of our knowledge, the existing literature has focused on small dimensional Markov switching models, which are not applicable to high dimensional cross-sections. We aim at filling a gap in the literature by studying Markov switching models as applied to large panels.

There now exists strong empirical evidence that macroecononomic and financial variables exhibit an approximate factor structure, as stressed in giannone_lenza_primiceri_2021. This nature of the data naturally leads to approximate latent factor specifications as a tool to model time series comovement in large dimensional cross-sections. For example, following the seminal contribution of chamberlainrothschild83, static approximate factor representations have been considered in Connor_Korajczyk_1986_JFE to develop measures of portfolio performance, and in stockwatson02JASA,stockwatson02JBES to forecast large macroeconomic panels and to build indexes of macroeconomic activity. The full inferential theory is developed by bai03. Settings allowing for dynamic factor representations have been also extensively studied: see FHLZ17 and references therein. A broad overview of large factor models is provided in stockwatson16. To the very best of our knowledge, the vast majority of existing contributions has looked at the linear setting. However, this may not be flexible enough to accommodate the discrete regimes typically observed in macroeconomic and financial series.

A number of contributions have extended linear static factor models to allow for discrete shifts in the loadings by assuming that these shifts are driven by an observable state variable. A first and growing stream of literature assumes that this state variable is a deterministic time index, which leads to a factor model with structural instability in the loadings: see Breitung_Eickmeier_2011_JoE, corradi2014testing, baltagi2016estimation, cheng2016shrinkage, BCF18, Barigozzi_Trapani_2020_SPA, duan2022quasi, among others, and Bai_Han_2016_FEC for a survey of the literature. The presence of structural breaks implies that regime changes are not recurrent and are related to events such as technological changes or shifts in monetary policy regimes. Alternatively, the states could be driven by the realisation of an observable stationary variable with respect to a reference value, in which case a threshold factor model would arise: see Massacci_2017_JoE,Massacci_2023_JFEcon. Under this set up, regimes are recurrent and associated to cyclical events such as business and financial cycles. Smoothly varying loadings are considered in motta2011locally and Pelger_Xiong_2022_JBES. Finally, Chen_Chen_Chen_2023_WP follow Su_Wang_2017_JoE and propose a time-varying matrix factor model with smooth changes in the loadings driven by a time index.

In this paper, we are interested in large dimensional factor models in relation to recurrent regime changes. A major drawback of threshold factor models is that they require a priori identification of the state variable. This may lead to model misspecification and unreliable empirical findings should the wrong state variable be employed to identify the regimes. In order to overcome this problem, we resort to the two-state Markov switching model of GQ1973 with a latent state variable, and we extend it to allow for an underlying large dimensional factor structure. Within this setting, we make the following major methodological contributions: we propose an algorithm to estimate the conditional state probabilities, as well as the loadings and the factors; and we derive the asymptotic properties of the estimators for loadings and factors. Remarkably, our results do not require knowledge of the true number of factors in any regime, and they are robust to the number of factors being unknown and estimated. This is an important aspect of our paper. Estimating the number of factors is challenging in a linear setting, as evidenced by the high number of relevant contributions: baing02, ABC10 and ahnhorenstein13, develop model selection criteria; Kapetanios_2010_JBES, onatski10, and Trapani_2018_JASA, propose inferential procedures. Dealing with an unknown number of factors clearly becomes even more engaging in the presence of regimes driven by a latent state variable and it therefore is an important contribution of our paper.

To the very best of our knowledge, the literature on large dimensional Markov Switching factor models is still in its infancy. However, two existing contributions are important to discuss. First, Liu_Chen_2016_SS study a model similar to ours, but their definition of common factors differs from ours in that they consider factors that are pervasive along the time dimension rather than along the cross-sectional dimension. As a consequence, their idiosyncratic components are assumed to be white noise. Second, Urga_Wang_2022_WP study a set up similar to ours, with some important differences: they assume a priori knowledge of the number of factors; they consider a model with serially homoskedastic idiosyncratic components. In addition, the Maximum Likelihood estimation approach of Urga_Wang_2022_WP adapts the EM algorithm by RT82 and baili12 to the case of Gaussian mixtures, where the weights are given by the probability of the latent variables to be in a given regime. Furthermore, the fact that the proposed EM algorithm is just an approximation to Maximum Likelihood estimation is however not accounted for when deriving the asymptotic properties of the considered estimators, in other words no formal proof that such algorithm is a contraction towards the Maximum Likelihood estimator is given.

Our approach is as follows. We introduce an algorithm to estimate factors, loadings, and transition probabilities, which extends to high dimensional factor models the state-space approach advanced in hamilton1989new and kim1994 to handle low dimensional Markov switching autoregressive models. In particular, we generalize the Baum-Lindgren-Hamilton-Kim filter and smoother, the original version of which was proposed to estimate Markov-switching VAR models: for example, see the reviews by Guidolin_2011_Chapter, krolzig2013markov, Hamilton_2016_Chapter, and guidolin2018essentials. An important feature of our approach is that it provides closed form expressions for all estimators. Even more remarkably, we not require a priori knowledge of the number of factors in each regime, which is instead needed by Urga_Wang_2022_WP.

We obtain our theoretical results by exploiting the well known property that a factor model with neglected discrete regime changes admits an equivalent representation with a higher number of factors: for example, see the discussions in Breitung_Eickmeier_2011_JoE, BCF18, and duan2022quasi, in the case of structural breaks; and Massacci_2023_JFEcon for threshold factor models. We use this property to estimate the latent factors by means of Principal Component Analysis (PCA) as applied to the linear representation. We then input these estimated factors into our algorithm, which allows us to recover the loadings and the transition probabilities. We then derive the asymptotic properties of the estimator for the loadings: we prove the asymptotic normality; we characterise the bias, which is induced both by the well known identification problem, and by the incomplete information related to the underlying data generated process. We also study the asymptotic properties of the estimated factors, which are obtained by projecting the data onto the estimated loadings. We corroborate our theoretical results through a comprehensive set of Monte Carlo experiments, which confirm the good finite sample properties of the estimation procedure we propose.

Finally, we assess the empirical validity of our model through three applications to large U.S. datasets of stock returns, macroeconomic variables, and inflation indexes. Markov switching models have been widely used to capture the cyclical behaviour of small-dimensional portfolios of financial assets: see Guidolin_2011_Chapter, and AT2012, and references therein. We apply our Markov switching factor model to a large dimensional portfolio of financial assets: the results show that the regimes described by the model closely follow U.S. business cycle dynamics, and complement the findings in Massacci_Sarno_Trapani_2021_WP, who identify the regimes based on an observable state variable. We then consider a large set of U.S. macroeconomic variables, and we use them to identify turning points in the U.S. business cycle in the spirit of Burns_Mitchell_1946_measuring: through appropriate metrics, we show that our model performs very well also on this respect. Finally, building upon the recent contribution of AL, we illustrate how our model may be employed to identify regimes in a large set of inflation indexes. Overall, these results confirm the usefulness of our theoretical framework to conduct empirical analysis.

The rest of the paper is organised as follows. Section (ref) introduces the two-state model. Section (ref) describes the estimation algorithm. Section (ref) derives the asymptotic theory. Section (ref) presents two further results related to estimation of the number of factors and to underspecification of the number of regimes. Section (ref) deals with the issue of unobserved heterogeneity. Section (ref) discusses the problem of testing for regime changes. Section (ref) runs a comprehensive set of Monte Carlo experiments. Section (ref) presents the empirical applications. Finally, Section (ref) concludes. Details about the estimation algorithm are given in Appendix (ref). Mathematical derivations are collected in Appendices (ref) and (ref). Additional Monte Carlo and empirical results are to be found in Appendices (ref) and (ref), respectively.

Notation

We denote as $\otimes$ the Kronecker product, with $\odot$ the element-wise (Hadamard) product, and with $\oslash$ the element-wise ratio. For a vector $\bm v=(v_1\cdots v_m)'$ we denote its Euclidean norm as $\Vert \bm v\Vert=\sqrt{\sum_{i=1}^m v_i^2}$. For a matrix $\mathbf C$ we denote the spectral norm as $\Vert\mathbf C\Vert=\sqrt{\mu_1(\mathbf C\mathbf C^\prime)}$, where $\mu_1(\mathbf C\mathbf C^\prime)$ indicates the largest eigenvalue of $\mathbf C\mathbf C^\prime$. If $\text{rk}(\mathbf C)=r<\infty$, then, we sometimes use the same notation $\Vert\mathbf C\Vert$ to denote also the Frobenius norm $\Vert\mathbf C\Vert_F=\sqrt{\text{tr}(\mathbf C\mathbf C^\prime)}$. Indeed, $\Vert\mathbf C\Vert_F\le \sqrt r \Vert\mathbf C\Vert$ and since it is always true that $\Vert\mathbf C\Vert\le \Vert\mathbf C\Vert_F$, then, bounding the Frobenius or the spectral norm is asymptotically equivalent.

For a scalar discrete random variable $Z$, the notation $\mathsf{P}( Z= z)$ is its probability mass function computed using the true value of the parameters. For random variables $\mathbf Y$ and $\mathbf W$ the notations $\mathsf{E}[\mathbf Y]$ and $\mathsf{E}[\mathbf Y|\mathbf W]$ are the expectation and conditional expectation given $\mathbf W$, respectively, computed with respect to the true distributions $F_Y(\mathbf y)$ and $F_{Y|W}(\mathbf y|\mathbf W)$ which in turn are computed using the true value of the parameters. If, in place of the true value of the parameters, we use an estimate of the parameters, say $\widehat{\theta}$, then we adopt the notations $\mathsf{P}_{\widehat{\theta}}( Z= z)$, $\mathsf{E}_{\widehat{\theta}}[\mathbf Y]$, and $\mathsf{E}_{\widehat{\theta}}[\mathbf Y|\mathbf W]$, respectively.

Finally, we let $\mathbf I_m$ be the identity matrix of dimension $m$, $\bm\iota_m$ an $m$-dimensional vector of ones, and $\mathbf 0$ any matrix or vector of zeros whose dimensions depend on the context.

Markov switching factor model

Setup

We study a two-state large dimensional Markov switching factor model. Formally, we consider

We assume that the elements of the $N\times1$ vector process of observable dependent variables $\{\mathbf{x}_{t}\}$ have zero mean, and we consider the more general case in which they are allowed to have mean different from zero in Section (ref); $\{\mathbf{f}_{jt}\}$ is the $r_{j}\times1$ vector process of latent factors such that $r_j$ is fixed and $r_{j}\ll N$, for $j=1,2$; $\bm\Lambda_j$ is the $N\times r_{j}$ matrix of factor loadings with rows equal to $\bm\lambda_{ji}^\prime$, for $i=1,\ldots, N$ and $j=1,2$; $\{\mathbf{e}_{t}\}$ is the $N\times1$ vector process of idiosyncratic components with innovations $\bm{\nu}_{t}\sim\left(\mathbf{0},\mathbf{I}_{N}\right)$. Note that we allow the elements of $\{\mathbf e_t\}$ to be both serially and cross-sectionally weakly correlated, and we refer to Section (ref) for the specific assumptions. It is also important to point out that the number of factors $r_j$ within each state is allowed to be unknown.

The model in (ref) and (ref) explicitly allows for two regimes: the case in which the number of states is actually underspecified is dealt with in Section (ref). Also, the number of factors $r_1$ and $r_2$ is allowed to change between the regimes: in this, our approach is more general than in Liu_Chen_2016_SS, who assume that $r_1=r_2$ and the dimension of the factor space is a priori the same between the two regimes.

As it is standard in the literature, we assume that $s_{t}$ follows a discrete-state, homogeneous, irreducible and ergodic, first-order Markov chain such that

with matrix of transition probabilities

Defining the $2\times1$ vector of state indicators

allows us to write the transition equation

where $\{\mathbf{v}_{t}\}$ is a discrete-valued zero mean martingale difference sequence whose elements sum to zero. Because, $\Vert \mathbf P\Vert <1$, $\{s_t\}$ follows an ergodic Markov chain, thus, there exists a stationary vector of probabilities $\bar {\bm\xi}$ satisfying: \[ \bar {\bm\xi}=\mathbf P'\bar {\bm\xi}. \] Hence, the elements of $\bar {\bm\xi}$ are long-run or unconditional state probabilities. In particular, we have $\bar {\bm\xi}=\mathsf{E}[\bm\xi_t]$, such that

where $0<\mathsf{P}(s_t=j)<1$, for $j=1,2$, by Assumption (ref) in Section (ref) below, which makes the Markov chain irreducible. In particular, (ref) and (ref) are related by guidolin2018essentials

Finally, unlike the low-dimensional model of diebold78measuring, we do not specify the factor dynamics. In particular, diebold78measuring allow for regime-specific factor mean, whereas the loadings do not vary: in this setting, the variance of the dependent variables remains constant over time. On the other hand, the large-dimensional model in (ref) and (ref) allows for regime-specific covariance matrix of $\mathbf{x}$: this is relevant for modelling both macroeconomic variables and financial returns, as stressed in McConnell_Perez-Quiros_2000_AER, and Perez_Quiros_Timmermann_2000_JF,Perez_Quiros_Timmermann_2001_JoE, respectively. We exploit this feature in the empirical analysis in Section (ref), where we use the model in (ref) and (ref) to study large U.S. datasets of stock returns, macroeconomic variables, and inflation indexes. On the other hand, we explain in Section (ref) how we can deal with datasets displaying regime-specific individual effects.

State space representation

Let the $\left(r_{1}+r_{2}\right)\times1$ vector process $\{\mathbf{g}_t\}$ be defined as

Let $\mathbf{B}_{1}=[\bm\Lambda _{1}\ \mathbf{0}]$ and $\mathbf{B}_{2}=[\mathbf{0}\ \bm\Lambda_{2}]$, where $\mathbf{B}_{1}$ and $\mathbf{B}_{2}$ are $N\times \left(r_{1}+r_{2}\right)$ matrices. The model in (ref), (ref) and (ref) admits the equivalent state space representation\footnote{Note that $\bm{\xi }_{t}\otimes \mathbf{g}_{t}= [\mathbf f_{1t}'~\mathbf 0_{}~\mathbf f_{2t}'~\mathbf 0_{}]'$.}

Under standard assumptions, the term $\left( \mathbf{B}_{1}~\mathbf{B}_{2}\right) \left( \bm{ \xi }_{t}\otimes \mathbf{g}_{t}\right) $ is identifiable up to a relabelling of the states. This means that the indices of the states can be permuted without changing the law governing the process for $\mathbf{x}_{t}$: on this, see Section 3 in LEROUX1992127. Also note that, even for given $\bm{\xi }_{t}$, identification of $\textbf{B}_{1}$ and $\textbf{B}_{2}$, and therefore of the elements of $\textbf{g}_{t}$, is in general possible only up to an invertible linear transformation (see bai03).

Linear representation

The model in (ref) admits the same equivalent linear representation as a model with either one change point or a single threshold effect: see BCF18, and Massacci_2017_JoE, respectively. It can then be rewritten as the $r_1+r_2$ linear factor model

where $\mathbf{A}=\left[\mathbf{\Lambda}_{1} ~ \mathbf{\Lambda}_{2}\right]$. Therefore, large dimensional factor models with two discrete regimes, be them modelled through a permanent structural change, or through cyclical threshold or Markov switching dynamics, admit the same equivalent linear representation. Then $\mathbf{A}$ and $\mathbf{g}_{t}$ may be estimated by standard Principal Component Analysis (PCA) stockwatson02JASA,stockwatson02JBES,bai03. Since PCA gives, as $N,T\to\infty$, consistent estimators of the factors up to premultiplication by an invertible matrix (see bai03), for ease of exposition we first consider estimation of the model in (ref) by treating $\mathbf g_t$ as known. We then briefly review the implementation of PCA and its effect on the estimation of the model in Section (ref).

Log-likelihood

Following the approaches by DGRqml, BLqml, and baili16, all developed for QML estimation of linear factor models, we consider a misspecified Gaussian quasi-likelihood of an exact factor model with white noise idiosyncratic components. This implies that the idiosyncratic components are treated as if they were cross-sectionally and serially uncorrelated. This approach is adopted also by Urga_Wang_2022_WP in the case of Markov switching factor models. It is important to stress that we are not assuming that the idiosyncratic components are uncorrelated, as we are just considering likelihood estimation of a misspecified model. Furthermore, in the linear case, baili16 and BLqml, show that such misspecifications are asymptotically negligible as $N,T\to\infty$.

The parameters of interest are then partitioned as

so that the vector of parameters of interest, denoted as $\mathbf{q}$, is defined as

Notice that we estimate only the diagonal elements of $\bm\Sigma _{e1}$ and $\bm\Sigma _{e2}$ in (ref). Let $\mathbf{X}=\left( \mathbf{x} _{1}^{\prime },\ldots ,\mathbf{x}_{T}^{\prime }\right) ^{\prime }$, $\bm{\mathcal G}=\left( \mathbf{g}_{1}^{\prime },\ldots ,\mathbf{g}_{T}^{\prime }\right) ^{\prime }$, where $\mathbf{X}$ is an $NT\times 1$ vector, $\bm{\mathcal G}$ is an $(r_1+r_2)T\times 1$ vector. These are $T$-dimensional realizations of the stochastic processes $\{\mathbf x_t\}$ and $\{\mathbf g_t\}$, respectively. Moreover, let $\bm{X}_{v}$ be the $\sigma$-algebra generated by the random variables $\{\mathbf{x}_{t}\}_{t=1}^v$, for $v=1,\ldots ,T$; in a similar way, define $\bm{G}_{v}$ as the $\sigma$-algebra generated by the random variables $\{\mathbf{g}_{t}\}_{t=1}^v$, for $v=1,\ldots ,T$. And for simplicity we write $\bm X\equiv \bm X_T$ and $\bm G\equiv \bm G_T$.

The likelihood function, denoted by $f\left( \mathbf{X};\mathbf{q}\right)$, can be decomposed as

in the last step we account for the fact that $f\left(\bm{\mathcal G};\mathbf{q}\right)\equiv f\left(\bm{\mathcal G}\right)$, since it does not depend on the parameters of our model, as we do not specify any dynamic model for the process $\{\mathbf g_t\}$.

Furthermore, following krolzig2013markov, we have

Here, to avoid heavier notation, we use the same notation $\{\bm\xi_t\}_{t=1}^T$ both for a generic $T$ dimensional realization of the process $\{\bm\xi_t\}$ and for the $\sigma$-algebra generated by the random variables $\{\bm\xi_t\}_{t=1}^T$. Notice that the sum is over $2^T$ possible values since, given a realization for $\{\xi_{1t}\}_{t=1}^T$, the realizations of $\{\xi_{2t}\}_{t=1}^T$ are given by $\xi_{2t}=1-\xi_{1t}$ for all $t$.

Given that we treat the idiosyncratic components as if they were uncorrelated, and using the Markov property of $\{\bm\xi_t\}$, up to omitted constant terms we have

where $\mathbf{\Sigma} _{et}=\left(\text{diag}(\mathbf{\Sigma} _{e1})~\text{diag}(\mathbf{\Sigma} _{e2})\right) \left( \bm{\xi } _{t}\otimes \mathbf{I}_{N}\right)$. Note that in this case the likelihood (ref) is not Gaussian; rather, it is a mixture of Gaussian distributions. Finally, again by the Markov property of $\{\bm\xi_t\}$, we can write

Estimation

In this section, we assume that the data generating process is characterised by two regimes as in the model in (ref) and (ref). In Section (ref) we study the case in which the model is underspecified and the data generating process exhibits a higher number of regimes. We also assume that the dimension of the vector $\mathbf{g}_{t}$ in (ref) is known. Should this not be the case, the dimension of $\mathbf{g}_{t}$ can be determined using information criteria such as those proposed in baing02, ABC10, and ahnhorenstein13, or inferential techniques such as those developed in onatski10 and Trapani_2018_JASA. This issue is discussed also in Section (ref).

In what follows, Section (ref) defines the steps of the proposed Expectation Maximization (EM) algorithm. Section (ref) describes the Baum-Lindgren-Hamilton-Kim filter and smoother. Section (ref) details the estimator for the factor space. Section (ref) discusses the estimator for the parameters. Section (ref) deals with initialization and convergence of the algorithm.

EM algorithm

The algorithm outlined in this section is a generalization of the procedure described by krolzig2013markov. The EM algorithm is made of two steps repeated at each iteration $k\ge 0$. The E step involves taking the expected value of the log-likelihood derived from (ref) conditional on $\bm{X}$ given an estimate of the parameters $\widehat{\mathbf{q}}^{\left(k\right)}$, namely

The M step solves the constrained maximization problem with respect to $\mathbf{q}=\left[ \bm\varphi^\prime,\bm\rho^\prime\right] ^{\prime }$, that is

where the constraints ensure that probabilities add up to one. In principle, in the M step we should also account for the term $\mathsf{E}_{\mathbf{\widehat{q}}^{\left( k\right) }}\left[ \log f \left( \bm{\mathcal G}\right)\left\vert\bm X\right. \right]$, which however in our context does not depend on any parameter.

It is well known that the iteration of these steps produces a series of increasing log-likelihoods. Indeed, $\mathsf{E}_{\mathbf{\widehat{q}}^{\left( k\right) }}\left[ \log f\left( \bm{\mathcal G}\left\vert \bm{X};\mathbf{q}\right. \right) \left\vert \bm X\right.\right]$ does not contribute to the convergence of the EM algorithm (see DLR77, and wu83). Moreover, if the maximum is identified and unique, then the EM algorithm will eventually lead to the Maximum Likelihood estimator of $\mathbf q$. As shown below, the solution of the M step can be computed explicitly using the expressions given in (ref) and (ref). This solution is unique and in closed form. Therefore, no identification issue arises due to multiple maxima, or related to the existence of such maxima.

Baum-Lindgren-Hamilton-Kim filter and smoother

From (ref) and (ref), in order to compute the expected likelihood in the E step we need to compute $\mathsf{E}_{\widehat{\mathbf q}^{(k)}}[\bm\xi_t|\bm X]$, $\mathsf{E}_{\widehat{\mathbf q}^{(k)}}[\bm\xi_t\otimes \mathbf g_t|\bm X]$, and $\mathsf{E}_{\widehat{\mathbf q}^{(k)}}[(\bm\xi_t\otimes \mathbf g_t)(\bm\xi_t\otimes \mathbf g_t)'|\bm X]=\mathsf{E}_{\widehat{\mathbf q}^{(k)}}[(\mathbf I_2\otimes \mathbf g_t\mathbf g_t')|\bm X]$.

We start by considering the case in which both $\{\mathbf g_t\}_{t=1}^T$ is observed and the true value of the parameters $\mathbf q$ is known, while we postpone the discussion of the estimation of the factors to Section (ref). Then, for the E step we just need to compute $\mathsf{E}[\bm\xi_t|\bm X]$, since in this case $\bm\xi_t$ and $\mathbf g_t$ are independent for all $t$. This is accomplished by means of a generalization the Baum-Lindgren-Hamilton-Kim filter and smoother explained in detail in Appendix (ref). It is an iterative procedure through which we first compute the sequences of conditional one-step-ahead predicted probabilities $\{\bm{\xi }_{t\left\vert t-1\right. }\}_{t=1}^{T}$, such that $\bm{\xi }_{t\left\vert t-1\right .}=\mathsf{E} \left[ \bm{\xi }_{t}\left\vert \bm{X}_{t-1}\right. \right]$, and filtered probabilities $\{\bm{\xi }_{t\left\vert t\right. }\}_{t=1}^{T}$ such that $\bm\xi_{t|t}=\mathsf{E}[\bm\xi_t|\bm X_t]$. Second, by means of those sequences, we compute the sequence of smoothed probabilities $\{\bm{\xi }_{t\left\vert T\right. }\}_{t=1}^{T}$ such that $\bm\xi_{t|T}=\mathsf{E}[\bm\xi_t|\bm X]$.

The final recursions for the filtered probabilities are given by (e.g., see krolzig2013markov, and hamilton1989new)

where

The filter can be started by setting either $\bm\xi_{0|0}=\left[ 1~0\right]^{\prime }$, or, equivalently, $\bm\xi_{0|0}=\left[ 0~1\right]^{\prime }$.

The final recursions for the smoothed probabilities are given by (e.g., see krolzig2013markov, and kim1994)

This backward recursion is initiated at $\bm\xi_{T|T}$, which is the last iteration of the filter in (ref).

The above description of the Baum-Lindgren-Hamilton-Kim filter and smoother assumes that $\mathbf q$ and $\mathbf g_t$ are observed. However, in practice both need to be estimated. This is discussed in the next two Sections (ref) and (ref) below.

Estimating the factor space

In order to estimate the factors $\mathbf g_t$, and their dimension $r_1+r_2$, we exploit the fact that the Markov switching factor model in (ref) is observationally equivalent to a linear factor model with $r_1+r_2$ common factors $\mathbf g_t$ and factor loadings $\mathbf A$: see Section (ref) and, in particular, equation (ref). The number of factors in (ref) can be estimated using methods already available in the literature: for example, see baing02, onatski10, ahnhorenstein13, and Trapani_2018_JASA. The factors $\mathbf g_t$ can be estimated by PCA as follows. First, the estimator $\widehat{\mathbf A}$ of the loadings matrix $\mathbf A$ is obtained as $\sqrt N$ times the normalized eigenvectors corresponding to the $r_1+r_2$ largest eigenvalues of the sample $N\times N$ covariance matrix $T^{-1}\sum_{t=1}^T \mathbf x_t\mathbf x_t'$. Second, the factors are estimated by linear projection of the data $\mathbf x_t$ onto the estimated loadings:

This is the same approach followed by stockwatson02JASA. It is also the dual approach of the one adopted by bai03. Consistency of $\widehat{\mathbf A}$ and $\widehat{\mathbf g}_t$ follow from Lemma (ref) and Lemma (ref)(a) in Appendix (ref), respectively. Note that the steps described in this section do not require knowing the latent state indicator $\bm\xi_t$, and they can be carried out independently. Because of these results, $\bm\xi_t$ and $\widehat{\mathbf g}_t$ can also be treated as independent for all $t$. As a consequence, the Baum-Lindgren-Hamilton-Kim filter described in Section (ref) can be implemented by just replacing the true factors $\mathbf g_t$ with their estimator $\widehat{\mathbf g}_t$ defined in (ref).

Estimating the parameters

At each iteration $k\ge 0$ of the EM algorithm, the filtered and smoothed probabilities, given in (ref) and (ref), respectively, and the smoothed cross-probabilities given in (ref), are computed using an estimator $\widehat{\mathbf q}^{(k)}$ of the parameters and an estimator $\widehat{\mathbf g}_t$ of the factors. Hereafter, we denote as $\bm\xi_{t|t}^{(k)}$, $\bm\xi_{t|T}^{(k)}$, and $\bm\xi_{t,t-1|T}^{(k)}$ such estimators. This defines the E step.

In the M step we have to solve the constrained maximization problem in (ref). Here we just give the final results, while we refer to Appendix (ref) for their derivation. The estimates of the loadings $\mathbf B_j$, $j=1,2$, are given by

and, consistently with the fact that we use a mis-specified likelihood with uncorrelated idiosyncratic components, we set

where $\mathbf{\widehat{b}}_{ji}^{(k+1)\prime}$ is the $i$th row of $\mathbf{\widehat{B}}_{j}^{(k+1)}$. Concerning the estimates of $\bm\rho$, which are subject to the adding up condition,

By letting $k^*$ be the last iteration of the EM algorithm, we define our final estimator of the parameters as $\widehat{\mathbf q}\equiv \widehat{\mathbf q}^{(k^*+1)}$, as given by (ref), (ref), and (ref). The final estimator of $\bm\xi_t$ is defined as $\widehat{\bm\xi}_{t|T}\equiv \bm\xi_{t|T}^{(k^*+1)}$, i.e., obtained by running one last time the Baum-Lindgren-Hamilton-Kim filter using the final estimates of the parameters.

Initialization and convergence of the EM algorithm

To start the algorithm we need initial estimators $\widehat{\mathbf q}^{(0)}$ for the parameters. Specifically, we set $\widehat{\mathbf B}_1^{(0)}=\widehat{\mathbf B}_2^{(0)}=\widehat{\mathbf A}$, as defined in Section (ref). Then, given also $\widehat{\mathbf g}_t$ as in (ref), let $\widehat{\mathbf e}_t=\mathbf x_t- \widehat{\mathbf A}\widehat{\mathbf g}_t$, and we set $\widehat{\bm\Sigma}_{e1}^{(0)}=\widehat{\bm\Sigma}_{e2}^{(0)}=\text{diag}\left(T^{-1}\sum_{t=1}^T \widehat{\mathbf e}_t\widehat{\mathbf e}_t'\right)$. Finally, we set $$ \widehat{\mathbf P}^{(0)}= \left(

\right), $$ where $\omega_1,\omega_2\in(0,0.5)$ and $\omega_1>\omega_2$. This initialization implicitly identifies state 1 as the most probable one, i.e., it is the state with largest unconditional probability as defined in (ref).

We say that the EM algorithm converged at iterations $k^*$, where $k^*$ is the first value of $k$ such that: \[ \frac{\left\vert\log f\left( \mathbf{X}\left\vert \bm{G};\widehat{\bm\varphi}^{(k)},\widehat{\bm\rho}^{(k)}\right. \right)- \log f\left( \mathbf{X}\left\vert \bm{G};\widehat{\bm\varphi}^{(k-1)},\widehat{\bm\rho}^{(k-1)}\right. \right) \right\vert} {\frac 12 \left\{\vert\log f\left( \mathbf{X}\left\vert \bm{G};\widehat{\bm\varphi}^{(k)},\widehat{\bm\rho}^{(k)}\right. \right)+ \log f\left( \mathbf{X}\left\vert \bm{G};\widehat{\bm\varphi}^{(k-1)},\widehat{\bm\rho}^{(k-1)}\right. \right) \right\}} < \epsilon, \] for some a priori chosen threshold $\epsilon>0$.

Asymptotic theory

In what follows, Section (ref) states the assumptions, whereas Section (ref) presents the asymptotic properties of the estimators.

Assumptions

For ease of reference, let us write (ref) and (ref) in scalar notation as

We consider the following set of assumptions, which generalizes to our framework the settings in bai03 and Massacci_2017_JoE.

Assumption (ref) restricts the factor processes $\left\{\mathbf{f}_{jt}\right\}$, for $j=1,2$, so that appropriate moments exist. The sequence $\{h_{kt}\}_{t=1}^T$ can be random or deterministic, and it is introduced to account for the fact that we estimate the expected value of $\xi_{jt}$, and not its actual value. Assumption (ref) implies that $0<\mathsf{P}\left[s_t=j\right]<1$, for $j=1,2$, thus ruling out the possibility that any of the states is absorbing, as discussed in Section (ref). It also implies that for $j=1,2$, as $T\rightarrow \infty $,

where $\mathbf{\Sigma }_{\mathbf{f}j}$ is positive definite and

In particular, note that (ref) allows the covariance matrix of $\mathbf{f}_j$ to be state-dependent, as advocated in Massacci_2023_JFEcon. It is also easy to see that if $j\ne k$, then for all $T\in\mathbb N$

According to Assumption (ref), loadings are nonstochastic and factors have a nonnegligible effect on the variance of $\{\mathbf{x}_{t}\}$ within each regime. In particular, part (b) implies that at least one common factor is present within each regime. The condition in part (d) ensures that the regimes are identified and it is analogous to the alternative hypothesis in the test for change in loadings developed in Pelger_Xiong_2022_JBES. This condition is trivially satisfied if $r_1 \neq r_2$, since the number of factors changes between regimes; if instead $r_1 = r_2$, then part (d) rules out the possibility that the columns of $\mathbf{\Lambda }_{1}$ are a linear combination of the columns of $\mathbf{\Lambda }_{2}$, in which case the regimes cannot be separately identified. From Assumption (ref) it also follows that, as $N\to\infty$,

and

Part (b) of Assumption (ref) controls the amount of cross-sectional correlation we can allow for. It implies the usual assumption for approximate factor models of nondiagonal idiosyncratic covariances $\bm\Sigma_{ej}$, $j=1,2$. Note that the sequence $\left\{h_{kt}\right\} _{t=1}^{T}$ has the same role as in Assumption (ref), which we refer to for further comments. Part (b) of Assumption (ref) also implies

and hence $N^{-1/2}\Vert \mathbb I(s_t=j) \mathbf e_t\Vert=O_p(1)$ for $j=1,2$, and for all $t\in\mathbb Z$. Part (c) of Assumption (ref) limits time dependence, and it is guaranteed together with part (a) if we assume finite 8th order cumulants for the bivariate process $\{(e_{it},e_{lt})\}$. Notice that the constant $M$ in the three parts of the assumption does not have to be the same one.

Assumption (ref) limits the degree of dependence between factors, state variable $s_t$, and idiosyncratic components.

Assumption (ref) guarantees a unique limit for $N^{-1}\mathbf{A}^{\prime }\mathbf{\widehat{A}}$, as stated in Lemma (ref) in Appendix (ref). By assuming distinct eigenvalues, we can uniquely identify the space spanned by the eigenvectors, which are linear combinations of the columns of $\mathbf A$. Notice that $\bm\Sigma_{\mathbf g}$ is block diagonal because of (ref).

Assumptions (ref) to (ref) are sufficient to prove the consistency of the estimators we propose. In order to derive their asymptotic distributions, we further introduce the following Assumptions (ref) and (ref).

Parts (a) and (b) of Assumption (ref) are suitable moment bounds, whereas parts (c) and (d) are central limit theorems.

Assumption (ref) imposes standard restrictions on the convergence rates.

Define the $\left(r_{1}+r_{2}\right) \times \left(r_{1}+r_{2}\right)$ matrix $\mathbf{\widehat{H}}$ as

where $\mathbf{G}=\left( \mathbf{g}_{1},\ldots ,\mathbf{g}_{T}\right)$ and $\mathbf{\widehat{V}}$ is the $\left(r_{1}+r_{2}\right) \times \left(r_{1}+r_{2}\right)$ diagonal matrix containing the first $r_{1}+r_{2}$ eigenvalues of $\mathbf{\widehat{\Sigma}}_{\mathbf{x}}=\left( NT\right) ^{-1}\sum\nolimits_{t=1}^{T}\mathbf{x}_{t}\mathbf{x}_{t}^{\prime }$ sorted in decreasing order. In Lemma (ref) we prove that

where $\mathbf{V}$ is the $\left(r_1+r_2\right)\times\left(r_1+r_2\right)$ diagonal matrix of the first $\left(r_1+r_2\right)$ eigenvalues of $\mathbf{\Sigma }_{\mathbf{g}}^{1\left/ 2\right. }\mathbf{\Sigma }_{\mathbf{A}}\mathbf{\Sigma }_{\mathbf{g}}^{1\left/ 2\right. }$ in decreasing order, and $\mathbf{\Psi}$ is the corresponding matrix of eigenvectors such that $\mathbf{\Psi}^{\prime}\mathbf{\Psi}=\mathbf{I}_{r_1+r_2}$. Likewise define $\mathbf{Q}_{j}=p\lim_{N,T\rightarrow \infty }N^{-1}\mathbf{\Lambda}_{j}^{\prime }\mathbf{\widehat{A}}$, for $j=1,2$, which is an $r_j\times(r_1+r_2)$ matrix such that $\mathbf{Q}=\left[\mathbf{Q}_{1}^{\prime}\ \mathbf{Q}_{2}^{\prime}\right]^{\prime}$. Thus, by Lemma (ref) we have

where $\mathbf{\Psi}_{j}$ is the $r_j\times\left(r_1+r_2\right)$ matrix such that $\mathbf{\Psi}=\left[\mathbf{\Psi}_{1}^{\prime}\ \mathbf{\Psi}_{2}^{\prime}\right]^{\prime}$. Therefore, because of (ref), (ref), and by Lemma (ref) according to which $\widehat{\mathbf V}\overset{p}{\rightarrow } \mathbf V$,

Asymptotic results

For $j=1,2$, let $\mathbf{\widehat{B}}_{j}=\mathbf{\widehat{B}}_{j}^{(k^{*}+1)}$, where $k^{*}$ is the last iteration of the EM algorithm as defined in Section (ref). For given $j=1,2$ and $i=1,\ldots,N$, let $\mathbf{\widehat{b}}_{ji}$ be the estimator for $\mathbf{b}_{ji}$ such that $\mathbf{\widehat{B}}_{j}=[\mathbf{\widehat{b}}_{j1},\ldots,\mathbf{\widehat{b}}_{jN}]^{\prime}$ and $\mathbf{B}_{j}=[\mathbf{b}_{j1},\ldots,\mathbf{b}_{jN}]^{\prime}$. The following theorem states the asymptotic distribution of $\mathbf{\widehat{b}}_{ji}$.

Theorem (ref) shows that the estimator $\mathbf{\widehat{b}}_{k_{1}i}$ for $\mathbf{b}_{k_{1}i}$ is subject to two sources of bias. The first is standard and it is induced by the usual indeterminacy due to the latency of both factors and loadings, and it is captured by the invertible matrix $\mathbf{\widehat{H}}$ defined in ((ref)) (see bai03). If we assume $T^{-1}\sum_{t=1}^T \mathbf g_t\mathbf g_t^\prime=\mathbf I_{r_1+r_2}$, then $\widehat{\mathbf H}$ becomes a rotation, namely an orthogonal matrix. However, additional restrictions on the loadings are necessary to reduce $\widehat{\mathbf H}$ to the identity: for a discussion on identification of factors see inter alia baing13. The second source of bias is induced by $\mathbf{\widehat{I}}_{\mathbf{\widehat{\xi}}k_{1}}$ defined in (ref), which depends on the probability of the state being asymptotically correctly estimated. If the unconditional probability of being in state $k_1$ were correctly estimated with probability one, that is, if $\widehat{\xi}_{k_1,t|T}\stackrel{p}{\to} \mathbb I(s_t=k_1)$, as $N,T\to\infty$, then $\mathbf{\widehat{I}}_{\mathbf{\widehat{\xi}}k_{1}}\stackrel{p}{\to}\mathbf{I}_{r_1+r_2}$ and $\mathbf{\widehat{b}}_{k_{1}i}$ would consistently estimate a linear transformation of $\mathbf{b}_{k_{1}i}$.

Therefore, $\mathbf{\widehat{b}}_{k_{1}i}$ estimates a linear transformations of $\mathbf{b}_{k_{1}i}$ and $\mathbf{b}_{k_{2}i}$, with weights determined by $\mathbf{\widehat{I}}_{\mathbf{\widehat{\xi}}k_{1}}$ and $( \mathbf{I}_{r_1+r_2}-\mathbf{\widehat{I}}_{\mathbf{\widehat{\xi}}k_{1}})$, respectively. This second source of bias is due to the fact that the process $s_t$ is latent, and it is specific to Markov switching models. As such, it does not affect threshold or structural break models, in which the state is identified with probability one.

Theorem (ref) has implications for the estimation of the regime specific loadings $\bm\Lambda_j$, $j=1,2$. To see this, let $\widehat{\mathbf R}_{k}= \widehat{\mathbf H} \widehat{\mathbf I}_{\widehat\xi k}$, for $k=1,2$, and consider the partition

where $\widehat{\mathbf R}_{ k, j\ell}$, $k,j,\ell=1,2$ and $\widehat{\mathbf H}_{ j\ell}$, $j,\ell=1,2$, are $r_j\times r_\ell$. Then, from Theorem (ref), for any given $i=1,\ldots,N$, as $N,T\to\infty$, we obtain

and

This means that $r_1+r_2$ columns of $\widehat{\mathbf B}_j$, $j=1,2$, estimate two different linear transformations of the columns of $[\bm\Lambda_1 \,\bm\Lambda_2]$. We can distinguish two cases. On the one hand, if $r_1=r_2=r$, as assumed for example in Liu_Chen_2016_SS, there is no need to know the true values of $r_1$ and $r_2$ to get consistent estimates of the space spanned by the true loadings in the two different regimes. Indeed, in this case $\mathbf B_1$ and $\mathbf B_2$ have an even number of columns, equal to $2r$, and from the first line of ((ref)) and ((ref)) we see that we can consider the first half of the columns of either $\widehat{\mathbf B}_1$ or $\widehat{\mathbf B}_2$ as an estimator of a linear transformation of $\bm\Lambda_1$ and the second half of the columns of either $\widehat{\mathbf B}_1$ or $\widehat{\mathbf B}_2$ as an estimator of a linear transformation of $\bm\Lambda_2$. Hence, we can define the following estimators of the loadings:

or

where $\widehat{\mathbf b}_{ji,1:r}$ denotes the first $r$ elements of $\widehat{\mathbf b}_{ji}$, and $\widehat{\mathbf b}_{ji,r+1:2r}$ denotes the second $r$ elements of $\widehat{\mathbf b}_{ji}$, for $j=1,2$ and $i=1,\ldots,N$. The property of these estimators are formalized in the following corollary, which is a direct consequence of Theorem (ref), and of ((ref)) and ((ref)).

This corollary has some interesting implications. If we strengthen Assumption (ref)(c) to add the identification constraint $\bm\Sigma_{\bm\Lambda_{12}}=\mathbf 0$, which is natural given Asssumption (ref)(d), then it is immediate to see that $\widehat{\mathbf H}_{12}\stackrel{p}{\to} \mathbf 0$ and $\widehat{\mathbf H}_{21}\stackrel{p}{\to} \mathbf 0$, as $N,T\to\infty$, in other words $\widehat{\mathbf H}\stackrel{p}{\to} \mathbf H$ which is now a block-diagonal matrix (see (ref) and recall that $\bm\Sigma_{\mathbf g}$ is block-diagonal by construction). It follows that if the unconditional probability of being in a given state were correctly estimated with probability one, so that, as $N,T\to\infty$, we had $\mathbf{\widehat{I}}_{\mathbf{\widehat{\xi}}k_{1}}\stackrel{p}{\to}\mathbf{I}_{r_1+r_2}$, then, as $N,T\to\infty$, for $k=1,2$ we have $\widehat{\mathbf R}_{k}\stackrel{p}{\to} \mathbf H$, which implies $\widehat{\bm{\lambda }}_{ki}^{\prime }\stackrel{p}{\to}\bm{\lambda }_{ki}^{\prime }\widehat{\mathbf{H}}_{kk}$, while $\widetilde{\bm{\lambda }}_{ki}^{\prime }\stackrel{p}{\to}\mathbf 0$. These results, which allow for a clear separation of $\bm\Lambda_1$ and $\bm \Lambda_2$, hold only under the restrictive assumption $\bm\Sigma_{\bm\Lambda_{12}}=\mathbf 0$. However, in general it is not possible to verify such condition and the two sets of estimators $\widehat{\bm{\lambda }}_{1i}^{\prime }$ and $\widehat{\bm{\lambda }}_{2i}^{\prime }$ or $\widetilde{\bm{\lambda }}_{1i}^{\prime }$ and $\widetilde{\bm{\lambda }}_{2i}^{\prime }$ will estimate consistently only a linear combination of the true loadings in both regimes.

On the other hand, if $r_1\ne r_2$, we need consistent estimators of $r_1$ and $r_2$ in order to be able to isolate the first $r_1$ columns of $\widehat{\mathbf B}_1$ and the last $r_2$ columns of $\widehat{\mathbf B}_2$, respectively. Therefore, if we only know that $r_1\ne r_2$ without knowing their true values, then we can consistently estimate a linear transformation of the columns of $\mathbf B_j$, but nothing can be said about $\bm\Lambda_j$, $j=1,2$.

Theorem (ref) describes the asymptotic properties of the estimator for the factor loadings $\mathbf{\widehat{B}}_{1}$ and $\mathbf{\widehat{B}}_{2}$. Complementary results can be obtained with respect to the estimated factors associated to the loading matrices $\mathbf{\widehat{B}}_{1}$ and $\mathbf{\widehat{B}}_{2}$. Formally, the true factors that correspond to $\mathbf{B}_{1}$ and $\mathbf{B}_{2}$ are $\xi _{1t}\mathbf{g}_{t}$ and $\xi _{2t}\mathbf{g}_{t}$, respectively, and their estimators are $\widehat{\xi}_{1,t\left\vert T\right. }\mathbf{\widehat{g}}_{t}$ and $\widehat{\xi}_{2,t\left\vert T\right. }\mathbf{\widehat{g}}_{t}$, respectively. The following theorem states the asymptotic distribution of these estimators.

In general, $\widehat{\mathbf{I}}_{\widehat{\mathbf{\xi }}j}\ne \mathbf I_{r_1+r_2}$ and so also $\mathbf{I}_{\mathbf{\xi }j}\ne \mathbf I_{r_1+r_2}$. Then, because of Theorem (ref), the estimator $\widehat{\mathbf b}_{ji}$ is biased and it is straightforward to see that the asymptotic covariance in Theorem (ref) is positive definite. Note that if we know that $r_1=r_2=r$ holds, then we can build consistent estimators for linear combinations of ${\mathbf f}_{jt}$, $j=1,2$, by simply regressing $\mathbf x_t$ onto the estimators $\widehat{\bm\Lambda}_j$ or $\widetilde{\bm\Lambda}_j$ which are defined in (ref) and (ref), respectively, and, as shown in Corollary (ref), are consistent for linear transformation of ${\bm\Lambda}_j$. Formally, this means we can build the sequence of factor estimators by running the cross-sectional regressions

or

If the unconditional probability of being in a given state is correctly estimated then $\widehat{\mathbf{I}}_{\widehat{\mathbf{\xi }}j}\stackrel{p}{\to} \mathbf I_{r_1+r_2}$ as $N,T\to\infty$, and Theorem (ref) is redundant: in this case, asymptotic normality of (ref) and of (ref) follows from arguments analogous to those in bai03. In the more general case we are considering, the asymptotic distribution of $\widehat{\mathbf f}_{jt}$ is stated in the following theorem (an analogous result holds for $\widetilde{\mathbf f}_{jt}$ and it is omitted for brevity).

According to Theorem (ref), $\widehat{\mathbf{f}}_{jt}$ estimates the space spanned by either $\mathbf{f}_{jt}$ or $\mathbf{f}_{kt}$, for $j,k=1,2$, with $j \neq k$, depending on which the true underlying regime is in period $t$.

On the number of factors and regimes

This section deals with two further issues related to the model in (ref) and (ref). Section (ref) studies estimation of the number of factors within each regime. Section (ref) discusses the consequences of an underspecified model.

Estimating the number of factors within each regime

Theorems (ref) and (ref) rely on the factor estimator $\mathbf{\widehat{g}}_{t}$ obtained from the equivalent linear representation in (ref). This estimator does not embed any information related to the likelihood of observing a regime $j$ at a given point in time $t$, for $j=1,2$ and $t\in\mathbb Z$. We now study the property of the estimator for the dimension of the factor space that is obtained when such information is accounted for. In particular, we are interested in separately identifying the number of factors within each regime, namely $r_{1}$ and $r_{2}$, given the dimension $r_1+r_2$ of the factor space of the equivalent linear representation in (ref). Note that under Assumption (ref)(b), at least one factor is present in each regime, which means that $r_1\geq1$ and $r_2\geq1$. Our framework is then more general than Liu_Chen_2016_SS and Urga_Wang_2022_WP: in the former $r_1=r_2$, and the two regimes have the same number of factors; the latter assumes that $r_1$ and $r_2$ are both known and do no have to be estimated. We do not impose any restriction on $r_1$ and $r_2$, except that $r_1\geq1$ and $r_2\geq1$, as required in Assumption (ref)(b). This is the natural extension of the linear set up, and it is aligned to Assumption B in baing02.

Formally, for $j=1,2$, we consider the regime-specific covariance matrix

where $0 < \sum\nolimits_{t=1}^{T}\widehat{\xi}_{jt\left\vert T\right. } < T$. The matrix $\mathbf{\widehat{\Sigma}}_{\widehat{\xi},\mathbf{x}j}$ includes information about the regimes through the estimated sequence $\{\widehat{\xi}_{jt\left\vert T\right. }\}_{t=1}^{T}$. Define the $r_{j}\times 1$ vectors

and the $r_{j}\times T$ matrices

For $1 \leq p \leq \bar{p}$, with $\bar{p} < \infty$, let $\mathbf{\widehat{V}}_{\widehat{\xi},j}^{\left( p\right) }$ be the $p \times p$ diagonal matrix containing the first $p$ eigenvalues of $\mathbf{\widehat{\Sigma}}_{\widehat{\xi},\mathbf{x}j}$ in decreasing order. Finally, let $\mathbf{\widehat{\Lambda}}_{\widehat{\xi},j}^{\left(p\right) }=[\bm{\widehat{\lambda}}_{\widehat{\xi},j1}^{\left( p\right) },\ldots,\bm{\widehat{\lambda}}_{\widehat{\xi},jN}^{\left( p\right) }]^{\prime}$ be the $N \times p$ matrix estimator for $\bm{\Lambda}_{j}$, which is obtained as $\sqrt N$ times the normalized eigenvectors corresponding to the $p$ largest eigenvalues of the $N\times N$ sample covariance matrix $\mathbf{\widehat{\Sigma}}_{\widehat{\xi},\mathbf{x}j}$ in (ref). The following theorem characterises the mean square convergence of $\bm{\widehat{\lambda}}_{\widehat{\xi},ji}^{\left( p\right) }$ for a given value of $p$.

Theorem (ref) extends Theorem 1 in baing02 and Theorem 3.4 in Massacci_2017_JoE to the case of the Markov switching factor model in (ref) and (ref). For $j,k=1,2$ with $j\neq k$, the theorem shows that $\bm{\widehat{\lambda}}_{\widehat{\xi},ji}^{\left( p\right) }$ estimates a linear combination of the vector $\left(\bm{\lambda }^{\prime}_{ji},\bm{\lambda }^{\prime}_{ki}\right)^{\prime}$ and not just of $\bm{\lambda }_{ji}$. It implies that the dimension of the estimated underlying factor space is $r_{1}+r_{2}$ even when the available information about the regimes is accounted for. Imperfect knowledge of the regimes therefore leads to an enlarged factor space: this makes our setting analogous to large dimensional change point factor models, as previously discussed in Section (ref). This complements what proved in Breitung_Eickmeier_2011_JoE, and corradi2014testing, who show that model misspecification in the form of omitted discrete regime shifts leads to an inflated number of factors. More generally, Theorem (ref) implies that, without further assumptions on the number of factors within each regime, it is not possible to separately estimate $r_1$ and $r_2$ even when the dimension $r_1+r_2$ of the equivalent linear representation in (ref) has been accurately estimated.

As in Liu_Chen_2016_SS, we now make the additional assumption that $r_1=r_2$, which means that the number of factors is equal across regimes. If the estimated number of factors in the equivalent linear representation in (ref) is an even number, we can recover the number of factors within each regime, as this is equal to $r_1=r_2=\left(r_1+r_2\right)/2$. On the other hand, if the estimated number of factors in the linear representation in (ref) is an odd number, an additional third regime might actually be neglected, as discussed in Section (ref) below.

Finally, under the assumption that both $r_1$ and $r_2$ are known as in Urga_Wang_2022_WP, the number of factors is known in both regimes and does not have to be estimated.

The case of an underspecified number of regimes

Up to know we have a priori assumed that the data are generated according to the model with two regimes in (ref) and (ref). This is consistent with existing empirical studies employing Markov switching models: for example, see diebold78measuring. However, in some cases the underlying data generating process of the dependent variables of interest displays a higher number of regimes: for example, Guidolin_Timmermann_JAE_2006 show that the joint distribution of stock and bond returns requires a four-state model. Therefore, the two-regime specification in (ref) and (ref) leads to model misspecification in case the joint distribution of the dependent variables $\mathbf{x}_{t}$ is characterised by a higher number of regimes.

We now study the case in which the model is underspecified and the data are generated by a process with a number of regimes that is finite and greater than two.

Since the number of regimes is finite, without loss of generality we consider the model with three regimes

and let

Suppose that only two regimes are accounted for. Given a natural ordering of the regimes, this means that we have to consider two cases, namely: $\left(a\right)$ $s_t=1$ and $s_t \neq 1$; $\left(b\right)$ $s_t=3$ and $s_t \neq 3$. The model in $\left( \ref{eq: reg_0}\right)$ admits the following two equivalent two-regime representations

where the loadings are defined as $\mathbf{B}_{1}^{\left( 1\right) }=\left( \mathbf{\Lambda }_{1}~\mathbf{0}~ \mathbf{0}\right)$, $\mathbf{B}_{2}^{\left( 1\right) }=\left( \mathbf{0}~ \mathbf{\Lambda }_{2}~\mathbf{\Lambda }_{3}\right)$, $\mathbf{B}_{1}^{\left( 3\right) }=\left( \mathbf{\Lambda }_{1}~\mathbf{ \Lambda }_{2}~\mathbf{0}\right)$, $\mathbf{B}_{2}^{\left( 3\right) }=\left( \mathbf{0}~\mathbf{0}~\mathbf{\Lambda }_{3}\right)$, the latent state process is defined as

the idiosyncratic covariance matrices are defined as $\bm{\Sigma}_{e1}^{\left( 1\right) }=( \bm{\Sigma} _{e1}~\mathbf{0}~\mathbf{0})$, $\bm{\Sigma} _{e2}^{\left( 1\right) }=( \mathbf{0}~\bm{\Sigma} _{e2}~\bm{\Sigma} _{e3})$, $\bm{\Sigma}_{e1}^{\left( 3\right)}=( \bm{\Sigma} _{e1}~\bm{\Sigma} _{e2}~\mathbf{0})$, $\bm{\Sigma} _{e2}^{\left( 3\right) }=( \mathbf{0}~ \mathbf{0~}\bm{\Sigma} _{e3})$, and the transition probabilities are equal to

For $j=1,3$, define the vector of parameters $\mathbf{q}^{\left(j\right)}=\left[ \bm\varphi^{\left(j\right) \prime},\bm\rho^{\left(j\right)\prime}\right] ^{\prime }$, where

Let $\left(NT\right)^{-1}\log f\left( \mathbf{X};\mathbf{q }^{\left( j\right) } \right) $ be the normalised log-likelihood function of $ \left( \ref{eq: reg_1}\right) $. Assume that

In a likelihood sense, the condition in ((ref)) captures a larger regime shift for $j=1$ than for $j=3$. Further, let $\widehat{\mathbf{q}}$ be the generic maximum likelihood estimator for the parameter of an underspecified model that allows for only two regimes when in fact the data generating process is given by (ref).

We proceed by contradiction, see also Appendix (ref) for more details. If $\widehat{\mathbf{q}}$ were an estimator for $\mathbf{q}^{\left( 3\right) }$, then

which leads to a contradiction since $\left(NT\right)^{-1}\log f\left( \mathbf{X};\widehat{\mathbf{q}}\right)$ is the estimated log-likelihood function. On the other hand, if $\widehat{\mathbf{q}}$ were an estimator for $\mathbf{q}^{\left( 1\right) }$, then

Therefore, when one regime is neglected, the maximum likelihood estimator estimates the regimes that maximise the likelihood according to the inequality in ((ref)). Provided that a sufficient number of iterations is done, the EM algorithm proposed in Section (ref) delivers an estimator that is close enough to the maximum likelihood estimator, such that the inequality in ((ref)) is preserved: see MR93,MR94. Therefore, the EM algorithm delivers the estimator for the underspecified representation that is associated to the highest likelihood. This also implies that when running the filter with just two regimes the estimated state $\widehat{\xi}_{1,t|T}$ is still correctly estimating the conditional expectation of the indicator related to the most likely regime, i.e., $\mathsf E[\mathbb I(s_t=1)|\bm X]$.

This result is consistent with the homologous finding in bai_1997, and Bai_Perron_1998_Econometrica, in relation to regression models with structural instability. Therefore, our result is the potential starting point for an inferential procedure on the number of regimes in large dimensional Markov switching factor models. It is also important to note that any neglected regime will be accounted for by an enlarged factor space, as discussed in Section (ref).

Unobserved heterogeneity

The model in $\left(\ref{eq:model}\right)$ assumes no individual effects. However, these may be important when modelling macroeconomic series as in diebold78measuring. In our set up, individual effects can be introduced by extending baili12,baili16 and considering

where $\bm{\alpha}_{j}=\left(\alpha_{j1},\ldots,\alpha_{jN}\right)^{\prime}$, for $j=1,2$, and $\alpha_{ji}$ captures the individual effect of cross-sectional unit $i$ within regime $j$. The vectors $\bm{\alpha}_{1}$ and $\bm{\alpha}_{2}$ introduce unobserved heterogeneity. If the state variable driving the regimes were observable, the resulting identification problem could be solved by expressing the model in terms of deviations of $\mathbf{x}_{t}$ from the conditional means within each regime: on this, see Massacci_Sarno_Trapani_2021_WP. However, since the state variable $s_{t}$ in $\left(\ref{eq:model_unobs_hetero}\right)$ is latent, this strategy no longer is applicable since the state is not observable with probability one. For this reason, we express the model in terms of the deviation of $\mathbf{x}_{t}$ from the unconditional mean.

Formally, consider the $N\times 1$ vector of centred variables $\mathbf{y}_{t}$ defined as

where $d_{jt}=\mathbb{I}\left( s_{t}=j\right) -\mathsf{E}\left[ \mathbb{I} \left( s_{t}=j\right) \right]$, $j=1,2$. If $\bm{\alpha }_{1}=\bm{\alpha }_{2}$, $\mathbf{x}_{t}$ has the same expected value in both regimes, and $\mathbf{y}_{t}=\mathbf{\Lambda }_{1}\mathbf{f}_{1t}\mathbb{I}\left( s_{t}=1\right) +\mathbf{\Lambda }_{2}\mathbf{f}_{2t}\mathbb{I}\left( s_{t}=2\right) +\mathbf{e}_{t}$. In the more general case in which $\bm{\alpha }_{1}\neq \bm{\alpha }_{2}$, unconditional demeaning leads to a larger factor space of dimension $r_{1}+r_{2}+2$. The additional two factors $d_{1t}$ and $d_{2t}$ take only two values, namely $d_{jt}= -\mathsf{E}\left[ \mathbb{I}\left( s_{t}=j\right) \right]$ or $d_{jt}= 1-\mathsf{E}\left[ \mathbb{I}\left( s_{t}=j\right) \right]$, depending on whether $\mathbb{I}\left( s_{t}=j\right)=0$ or $\mathbb{I}\left( s_{t}=j\right)=1$, respectively, for $j=1,2$. In this case, the equivalent linear representation in $\left(\ref{eq:linear_model}\right)$ holds with $\mathbf{g}_{t}=\left[ d_{1t},\mathbb{I}\left( s_{t}=1\right) \mathbf{f }_{1t}^{\prime },d_{2t},\mathbb{I}\left( s_{t}=2\right) \mathbf{f} _{2t}^{\prime }\right]$ and $\mathbf{A}=\left[ \bm{\alpha }_{1},\mathbf{\Lambda }_{1},\bm{\alpha }_{2},\mathbf{\Lambda }_{2}\right]$. The measurement equation in $\left(\ref{eq:state_space_mes}\right)$ of the state space representation remains valid with $\mathbf{B}_{1}=\left[ \bm{\alpha }_{1},\mathbf{\Lambda }_{1},\bm{\alpha }_{2},\mathbf{0}\right]$ and $\mathbf{B}_{2}=\left[ \bm{\alpha }_{1},\mathbf{0},\bm{\alpha }_{2},\mathbf{\Lambda }_{2}\right]$. Therefore, the tools developed in this paper can be applied to the sample counterpart of $\mathbf{y}_{t}$, namely to $\mathbf{\widehat{y}}_{t}=\mathbf{x}_{t}-\left(T^{-1}\sum_{t=1}^{T}\mathbf{x}_{t}\right)$, which consistently estimates $\mathbf{y}_{t}$ as $T\rightarrow \infty $. Corollary (ref) holds accordingly with respect to $\left(\alpha_{1i},\bm{\lambda}^{\prime}_{1i}\right)^{\prime}$ and $\left(\alpha_{2i},\bm{\lambda}^{\prime}_{2i}\right)^{\prime}$ instead of with respect to $\bm{\lambda}_{1i}$ and $\bm{\lambda}_{2i}$ only, respectively, for $i=1,\ldots,N$.

Detecting regime changes

The model in (ref) and (ref) a priori assumes the existence of two regimes. However, in practice Markov switching dynamics should be detected with suitable statistical tools. The development of rigorous inference goes beyond the purpose of this paper. In what follows, we give an overview of the relevant literature, which we use to discuss a possible starting point to run inference on the number of regimes in large dimensional Markov switching factor models.

First of all, it is however important to note that the Monte Carlo experiments in Section (ref) show that, when we fit the model in (ref) and (ref) to a linear factor model with just one regime (which means a model with no regime change), the algorithm detailed in Section (ref) assigns probability almost equal to unity to one state and therefore does not require any inferential procedure on the number of regimes. We refer to Appendix (ref) and the related Tables (ref) and (ref) for all relevant details.

As discussed in Qu_Zhuo_2021_ReStat, there exist three approaches to detect Markov regime switching in low dimensional models. A first one involves testing parameter homogeneity against heterogeneity: this is done in Carrasco_Hu_Ploberger_2014_Econometrica, who develop a class of tests for parameter constancy in random coefficient models; the power of these tests may however be limited, as they detect parameter heterogeneity of general form and are not specific to Markov switching models. A second approach, put forward in Hamilton_1996_JoE, proposes specification tests in Markov switching models: if the null hypothesis of correct model specification is rejected, as a solution one may include additional regimes; however, also this approach may suffer from low power, as it detects model misspecification of unknown form. Finally, a third approach proposes likelihood ratio based tests for the null hypothesis of a given number of regimes against the alternative of a higher number of regimes: this is followed in Hansen1992 and Qu_Zhuo_2021_ReStat, and it needs to account for the problem highlighted in Davies1977,Davies1987 as the additional transition probabilities are identified only under the alternative.

The above mentioned contributions are valid for low dimensional models. They are not directly applicable to large dimensional factor models, as these require imposing a number of restrictions on the loadings that goes to infinity as $N\rightarrow\infty$. This problem has been addressed when the variable driving the state is observable. Chen_Dolado_Gonzalo_2014_JoE, and Han_Inoue_2015_ET, test for a break in the loadings by testing for a change in the covariance matrix of the estimated factors. This approach, also used in Massacci_2017_JoE in threshold factor models, is valid provided that the covariance matrix of the true factors is stable over time. However, this may not be realistic in practice, as discussed in Chen_Dolado_Gonzalo_2014_JoE. Massacci_2023_JFEcon develops an inferential procedure for threshold factor models that is robust to factor heteroskedasticity. However, these solutions are not directly applicable to large dimensional Markov switching factor models, since the state variable is latent rather than observable.

Given the above discussion, a possible strategy to conduct inference on the number of regimes in large dimensional Markov switching factor models is to merge the tests available for low dimensional models with those in use for large dimensional factor models with observable state variable. This is a complex problem that goes beyond the purpose of this paper and will be addressed in future research.

Monte Carlo

We set $N=\{100,200\}$ and $T=\{250, 500, 750, 1000\}$. At each time period $t=1,\ldots, T$, we simulate the $N\times 1$ vector of data $\mathbf x_t$ according to (ref) and (ref). This requires to simulate the latent state $\bm \xi_t$, the loadings $\bm\Lambda_1$ and $\bm\Lambda_2$, the factors $\mathbf f_{1t}$ and $\mathbf f_{2t}$, and the idiosyncratic components $\mathbf e_t$.

We simulate the latent state $\bm \xi_t$ according to (ref), with $\mathbf P$ having entries $p_{11}=0.9$ and $p_{22}=0.7$, so that $p_{12}=0.1$ and $p_{21}=0.3$. This configuration corresponds to the unconditional probabilities to be equal to $\mathsf{P}(s_t=1)=\mathsf{E}[\xi_{1t}]=\frac{1-p_{22}}{2-p_{11}-p_{22}}=0.75$ and $\mathsf{P}(s_t=2)=\mathsf{E}[\xi_{2t}]=\frac{1-p_{11}}{2-p_{11}-p_{22}}=0.25$. Then, we generate the innovations $\mathbf v_t$ of the VAR in (ref) as follows: at each given $t$ we generate $u_t\sim \mathcal U[0,1]$ and

We set the number of factors in each state to $r_j=r=\{1,2\}$, $j=1,2$. The common component is generated according to model (ref). Let $\chi_{it}=\bm\lambda_{1i}^\prime \mathbf{f}_{1t}\mathbb I(s_t=1)+\bm\lambda_{2i}^\prime \mathbf{f}_{2t}\mathbb I(s_t=2)$, $i=1,\ldots, N$, $t=1,\ldots, T$. The $r$ entries of $\bm\lambda_{1i}$ and $\bm\lambda_{2i}$ are generated from a $\mathcal N(1,1)$ distribution. The matrices $\bm\Lambda_1$ and $\bm\Lambda_2$ are then transformed in such a way that $\bm\Lambda_1^\prime \bm\Lambda_1$ and $\bm\Lambda_2^\prime \bm\Lambda_2$ are diagonal matrices. The factors are such that $\mathbf f_{jt}=\mathbf f_t$, $j=1,2$, and satisfy $T^{-1}\sum_{t=1}^T\mathbf f_t\mathbf f_t^\prime =\mathbf I_r$, where each component of $\mathbf f_t$ is such that $f_{kt}=\rho_f f_{k,t-1}+z_{kt}$, $k=1,\ldots,r$, with $\rho_f=\{0,0.7\}$ and $z_{kt}\sim \mathcal N(0,1)$.

The idiosyncratic components are generated according to (ref), where $\mathbf \Sigma_{je}=\mathbf \Sigma_{je,a}+\mathbf \Sigma_{je,b}$, $j=1,2$, with $\mathbf \Sigma_{je,a}$ diagonal and $\mathbf \Sigma_{je,b}$ banded. Specifically, the entries of $\mathbf \Sigma_{1e,a}$ are generated from a $\mathcal U[0.25,1.25]$ and those of $\mathbf \Sigma_{2e,a}$ are generated from a $\mathcal U[0.75,1.75]$, while $\mathbf \Sigma_{1e,b}$ is a Toeplitz matrix with $\tau^{k}$ on the $k$th diagonal for $k=1,2$ and zero elsewhere, and, finally $\mathbf \Sigma_{2e,b}$ is a Toeplitz matrix with $\tau^{k-1}$ on the $k$th diagonal for $k=1,2,3$ and zero elsewhere. We set $\tau=\{0,0.5\}$. Moreover, each component of $\bm\nu_t$ is such that $\nu_{it}=\rho_{i}\nu_{i,t-1}+\omega_{it}$, $i=1,\ldots, N$, $t=1,\ldots,T$, with $\rho_i=\{0,\rho\}$ and $\rho\sim\mathcal U[0,0.5]$. Finally, we set the average noise-to-signal ratio across all $N$ simulated time series to be $N^{-1}\sum_{i=1}^N\frac{\sum_{t=1}^T e_{it}^2}{\sum_{t=1}^T \chi_{it}^2}=0.5$.

We simulate the model above 100 times for different values of $r$, $\rho_f$, $\tau$, and $\rho$. The EM is run allowing for at most 100 iterations and using a convergence threshold equal to $10^{-6}$. We initialize the algorithm using PCA as described in Section (ref). Since the states are identified only up to a permutation at each iteration of the algorithm we assign label 1 to the state with the highest estimated unconditional probability.\footnote{Note that the initialization such that $\omega_1=\omega_2=0.5$ is not empirically feasible, as it leads to no convergence of the EM algorithm. We conjecture that this has to do with the relabelling issue discussed in Section (ref), since for $\omega_1=\omega_2=0.5$ both states are equally likely.}

Results are collected in Tables (ref)-(ref) and are organised as follows: $\left(i\right)$ $r=1$, $\rho_f=0$, $\tau=0$, $\rho=0$ in Table (ref); $\left(ii\right)$ $r=1$, $\rho_f=0.7$, $\tau=0.5$, $\rho=0.5$ in Table (ref); $\left(iii\right)$ $r=2$, $\rho_f=0$, $\tau=0$, $\rho=0$ in Table (ref); $\left(iv\right)$ $r=2$, $\rho_f=0.7$, $\tau=0.5$, $\rho=0.5$ in Table (ref).

The first four columns of Tables (ref)-(ref) report the mean and, between brackets, the corresponding standard deviation over all replications of the estimated diagonal entries of the transition matrix $\widehat p_{jj}$, $j=1,2$, of the unconditional probabilities $\mathsf{P}(s_t=j)$, estimated as $\bar{\widehat{\xi}}_{j,t|T}=T^{-1}\sum_{t=1}^T \widehat{\xi}_{j,t|T}$, $j=1,2$.

Since the loadings are not identified, in the fifth column of Tables (ref)-(ref) we report the multiple $R^2$ coefficient obtained from regressing the columns of $\widehat{\mathbf B}_1$ onto the columns of $\mathbf B_1^*=\mathbf B_1 \mathbf{\widehat{I}}_{\mathbf{\widehat{\xi}}1}+ \mathbf B_2 (\mathbf I_{2r}-\mathbf{\widehat{I}}_{\mathbf{\widehat{\xi}}1})$, thus correcting for the bias described in Theorem (ref). Namely, we compute \[ R^2_{B^*} = \frac{\text{tr}\left\{ \left(\mathbf B_1^{*\prime} \widehat{\mathbf B}_1\right) \left(\widehat{\mathbf B}_1^{\prime} \widehat{\mathbf B}_1\right)^{-1} \left(\widehat{\mathbf B}_1^\prime\mathbf B_1^{*}\right) \right\}} {\text{tr}\left({\mathbf B}_1^{*\prime}{\mathbf B}_1^{*}\right)}. \] The closer this number is to one, the closer is the space spanned by the columns of $\widehat{\mathbf B}_1$ to the space spanned by the columns of ${\mathbf B}_1^*$ (see DGRqml).

In the sixth column of Tables (ref)-(ref) we report the MSE of the estimated common components defined as \[ \text{MSE}(\chi) = \frac{\sum_{i=1}^N\sum_{t=1}^T (\widehat \chi_{it}-\chi_{it})^2}{\sum_{i=1}^N\sum_{t=1}^T \chi_{it}^2}, \] where $\widehat \chi_{it} = \left( \widehat{\mathbf{b}}_{1i}~\widehat{\mathbf{b}}_{2i}\right)^\prime \left( \widehat{\bm{ \xi }}_{t}\otimes \widehat{\mathbf{g}}_{t}\right)$.

In the last column of Tables (ref)-(ref) we report the average number of iterations needed for the EM algorithm to converge.

The results in Tables (ref)-(ref) confirm the empirical validity of the estimation procedure detailed in Section (ref). In all four scenarios, as $N$ and $T$ increase the estimators $\widehat p_{11}$, $\widehat p_{22}$, $\bar{\widehat{\xi}}_{t|T,1}$ and $\bar{\widehat{\xi}}_{t|T,2}$ all converge to the true values of the corresponding parameters. In addition, $R^2_{B^*}$ and MSE($\chi$) are very to $1.00$ and $0.00$, respectively. Finally, note that the average number of iterations declines almost monotonically as $N$ and $T$ increase.

So far, the considered data generating process studies the performance of the proposed EM algorithm when in the model in (ref)-(ref) the loadings and idiosyncratic covariances are regime specific but the factors and their number do not change. We then consider three more scenarios which we briefly describe here while we refer to Appendix (ref) for details on the data generating process and simulation results.

First, we consider the same data generating process as the one considered in this section, but when setting a different number of factors in each regime, specifically, we set $r_1=3$ and $r_2=1$. We the run our EM algorithm initialized by means of PCA using $r_1+r_2=4$ factors. Results show that we correctly estimate the conditional and unconditional probabilities, as well as we correctly retrieve the loadings space (see Tables (ref) and (ref)).

Second, we set $r=r_j=1$, $j=1$, and we let only the autocorrelation of the factors be regime specific, while the loadings and idiosyncratic covariances are constant. In this case the EM algorithm wrongly overestimates the probability of being in the regime with highest simulated probability, thus it does not find evidence of a Markov switching dynamics, but it correctly retrieves the constant loadings space as the PCA estimator would do. Indeed, PCA is known to deliver consistent estimates of the loadings space even when the factors dynamics is piecewise constant BCF18,duan2022quasi (see Tables (ref) and (ref)).

Last, we simulate data from a linear factor model with $r=2$ factors, i.e., when no change is present, but then we fit on the same data our Markov switching model as if there were two regimes. The EM algorithm correctly assigns 97% probability to one regime at all time periods, i.e., as if there were just one regime (see Tables (ref) and (ref)).

Overall, our Monte Carlo findings provide evidence in support of the estimation algorithm proposed in Section (ref).

Empirical analysis

In this section we show how the methodological framework we propose can be used to model three different large U.S. datasets involving stock returns, macroeconomic time series, and inflation indexes. This is done in Sections (ref), (ref), and (ref), respectively. For each application, the estimated factors $\widehat{\mathbf f}_{jt}$, as defined in (ref) for $j=1,2$, are shown in Appendix (ref).

Stock returns

This application relates to a vast literature that models stock return dynamics using Markov switching specifications. Perez_Quiros_Timmermann_2000_JF,Perez_Quiros_Timmermann_2001_JoE document business cycle asymmetries in U.S. stock returns using decile-sorted portfolios. Ang_Bekaert_2002_RFS, and Guidolin_Timmermann_RFS_2008, study portfolio allocation in international equity markets under regime switching. In a multi asset setting, Guidolin_Timmermann_JAE_2006 describe the joint distribution of equity and bonds under regime switching. Guidolin_2011_Chapter, and AT2012, provide a review of the literature. We contribute to this literature by characterizing stock return dynamics using a Markov switching model in a large dimensional setting. To the very best of our knowledge, we are the first to do so.

The vector of observable dependent variables $\mathbf{x}_{t}$ in $\left(\ref{eq:model}\right)$ is made of monthly value weighted returns in excess of the risk-free rate from the $N=49$ industry portfolios kindly made publicly available on Kenneth French website.\footnote{See \url{https://mba.tuck.dartmouth.edu/pages/faculty/ken.french/data_library.html}.} Consistently with the discussion in Section (ref), the unconditional mean of $\mathbf{x}_t$ is equal to $\mathbf{0}$, which means that the returns have been demeaned along the time series dimension over the whole sample period. To obtain a balanced panel, the sample runs from July 1969 through December 2021, a total of $T=630$ time periods.

Using the eigenvalue ratio criterion of ahnhorenstein13 as applied to the equivalent linear representation in ((ref)), we find that the dimension of the vector $\mathbf{g}_t$ is equal to $r_1+r_2=2$ common factors. As commonly assumed in the related literature (see AT2012), we let the number of regimes be equal to two. Therefore, there is one common factor in each regime, so $r_1=r_2=r=1$. Based on this result, we apply the algorithm detailed in Section (ref). We stress that, in this case, it is crucial to allow for heteroskedastic idiosyncratic components, namely $\bm\Sigma_{e1}\ne \bm\Sigma_{e2}$ as assumed in the general model specification in (ref), since the idiosyncratic components on average account for about 35% of the total variation in the data. Given this set up, the EM algorithm converges in 22 iterations.

The realisation of the estimator $\widehat{\mathbf P}$ for the matrix of conditional probabilities $\mathbf P$ in ((ref)) is

The estimated unconditional probability for regime $j$ is equal to the sample average $\bar{\widehat{\xi}}_{j|T}=T^{-1}\sum_{t=1}^{T}\widehat{\xi}_{j,t|T}$, for $j=1,2$. It follows that $\bar{\widehat{\xi}}_{1|T}=0.8044$ and $\bar{\widehat{\xi}}_{2|T}=0.1956$.\footnote{The analytical formulas of the unconditional probabilities in (ref) give $\bar{\widehat{\xi}}_{1|T}=0.8081$ and $\bar{\widehat{\xi}}_{2|T}=0.1919$.} Therefore, regime $j=1$ is approximately four times more frequent than regime $j=2$. This lead us to label $\widehat{\xi}_{2,t|T}$ as the probability of a recession, since expansions occur more often than recessions.

Figure (ref) plots the sequences of estimates $\widehat{\xi}_{1,t|T}$ and $\widehat{\xi}_{2,t|T}$, for $t=1,\ldots,T$. In order to provide economic understanding of the regimes described by the model, we define the estimated recession indicator $\widehat{REC}_t$ as being equal to one if $\widehat{\xi}_{2,t|T}\geq 0.5$ and to zero otherwise. Formally, this means that $\widehat{REC}_t=\mathbb{I}\left(\widehat{\xi}_{2,t|T}\geq 0.5\right)$. Note that $\widehat{REC}_t$ has correlation equal to $0.99$ with $\widehat{\xi}_{2,t|T}$, which suggests that the underlying states are precisely estimated. We then follow harding2006synchronization and compute the degree of concordance between the estimated recession indicator and the NBER recession indicator, denoted as $REC_t$.\footnote{The NBER recession indicator is publicly available at \url{https://fred.stlouisfed.org/series/USREC}.} The degree of concordance is given by

For the dataset of stock returns we consider, we have $DoC=0.8048$. We also compute the probabilities of misclassification, which are given by $FP=T^{-1}\sum_{t=1}^T \widehat{REC}_t\, (1-REC_t)$ (namely, the frequency of false positives) and $FN=T^{-1}\sum_{t=1}^T (1-\widehat{REC}_t)\, REC_t$ (namely, the frequency of false negatives). We obtain $FP=0.1286$ and $FN=0.0667$. Therefore, the state $j=1$ is related to periods of economic expansions, whereas the state $j=2$ is more likely to occur during recessionary phases. Our model therefore captures regime changes in equity markets related to business cycle dynamics.

We then turn to the estimated factors. Since $r_1=r_2$, the estimators for $\bm\Lambda_j$, for $j=1,2$, are readily available from (ref) or (ref). Next, by projecting the data onto the estimated loadings weighted by the probability of being in a given state, we obtain the estimated scalar factors $\widehat{{f}}_{jt}$ and $\widetilde{{f}}_{jt}$, for $j=1,2$ and $t=1,\ldots,T$, as given in (ref) and (ref), respectively.

Table (ref) displays the correlations between the estimated latent factors and the six observable factors considered in FF2016RFS, namely: the value-weighted return on the market portfolio in excess of the one-month Treasury bill rate ($RM_t$); size ($SMB_t$); value ($HML_t$); profitability ($RMW_t$); investment ($CMA_t$); momentum ($MOM_t$). These correlations are computed both over the whole sample period, as well as within regimes. These in turn are defined in two ways: through the NBER recession indicator $REC_{t}$ (Panel A); through the predicted NBER recession indicator $\widehat{REC}_t$ previously defined (Panel B). The results in Table (ref) show that, over the whole sample period, $\widehat{{f}}_{1t}$ is strongly correlated with $RM_t$, and reasonably correlated with $SMB_t$, $HML_t$ and $CMA_t$. The estimate $\widehat{{f}}_{2t}$ is correlated with $MOM_t$. A similar picture comes from $\widetilde{{f}}_{1t}$ and $\widetilde{{f}}_{2t}$. When we compute the correlations during NBER expansions and recessions, additional findings arise (Panel A). On one hand, in expansionary periods, the correlations between $\widehat{{f}}_{1t}$ and $\widetilde{{f}}_{1t}$, and $RM_t$, $SMB_t$, $HML_t$ and $CMA_t$, are similar to those computed over the whole sample period. On the other hand, $\widehat{{f}}_{2t}$ and $\widetilde{{f}}_{2t}$ display sizeable correlations in recession with $SMB_{t}$ and $HML_{t}$, as well as with $MOM_{t}$. The homologous correlations calculated for the regime $j=2$ identified by the model are generally of lower magnitude, with the exception of those related to $MOM_{t}$ (Panel B). This confirms that $f_{2t}$ is a factor that drives the cross-section of equity returns during macroeconomic recessionary periods. Whereas a linear factor model would not be able to uncover this feature, our model can detect these asymmetric dynamics. This shows the empirical usefulness of our framework to model large dimensional portfolios of financial assets.

Macroeconomic time series

We now apply our methodology to a large set of macroeconomic variables to measure the probability of recessions and expansions in the U.S. economy. This relates our work to a large literature on business cycle dating, which goes back to the pioneering work of Burns_Mitchell_1946_measuring: see Romer_Romer_2020_WP for a recent discussion of the topic. We follow hamilton1989new, diebold78measuring, and Chauvet_1998_IER, in employing a Markov switching approach. In the spirit of Stock_Watson_2014_JoE, we use a large set of time series data to estimate recession and expansion probabilities. Finally, we study the ability of our model in dating turning points both using the full-sample and in real-time in a spirit similar to Chauvet_Piger_2008.

Formally, the vector of observable dependent variables $\mathbf{x}_{t}$ in $\left(\ref{eq:model}\right)$ is made of the monthly macroeconomic dataset FRED-MD described by mccracken2016fred formed of $N=126$ times series covering both the real and nominal sectors of the U.S. economy and including also labor market indicators, and financial variables.\footnote{See \url{https://research.stlouisfed.org/econ/mccracken/fred-databases/}.} The data is transformed to stationarity and missing values are imputed by means of the routines made available by mccracken2016fred, which produce a balanced panel, with a sample running from April 1959 through March 2024, for a total of $T=780$ time periods.

Using the information criterion of baing02 as applied to the equivalent linear representation in ((ref)), we find that the dimension of the vector $\mathbf{g}_t$ is equal to $r_1+r_2=8$ common factors. As commonly assumed in the literature Romer_Romer_2020_WP, we consider two regimes. Therefore, under the assumption that the number of factor is the same across states, there are four common factors in each regime, namely $r_1=r_2=r=4$. We then apply the algorithm detailed in Section (ref). We further impose homoskedastic idiosyncratic components, namely $\bm\Sigma_{e1}=\bm\Sigma_{e2}$. This is because, in the dataset in use, idiosyncratic components are often negligible, explaining on average less than 10% of the total variation of real variables BN06.\footnote{Results with heteroskedastic idiosyncratic components are similar and available upon request.} In this set up, the EM algorithm converges in 12 iterations.

The estimate of the matrix of conditional probabilities $\mathbf P$ in ((ref)) is equal to

The estimated unconditional probabilities are $\bar{\widehat{\xi}}_{1|T}=0.8354$ and $\bar{\widehat{\xi}}_{2|T}=0.1646$.\footnote{The analytical formulas in (ref) give unconditional probabilities equal to $\bar{\widehat{\xi}}_{1|T}=0.8362$ and $\bar{\widehat{\xi}}_{2|T}=0.1638$.} In this sample, the unconditional probability of a recession, as measured by the NBER recession indicator, is 0.1218. Therefore, we can identify regime $j=2$ as the recession regime.

Figure (ref) plots the sequences of estimates $\widehat{\xi}_{1,t|T}$ and $\widehat{\xi}_{2,t|T}$, for $t=1,\ldots,T$. The two most recent main recessions, which are due to the Great Financial Crisis (2007-2009) and the Covid19 pandemic (2020-2021), are well captured. To quantify the performance of our model, we once again follow harding2006synchronization and compute the degree of concordance $DoC$ in (ref) between the estimated recession indicator $\widehat{REC}_t$ defined as in Section (ref), and the NBER recession indicator. We obtain $DoC=0.7718$, with frequency of false positives and false negatives equal to $FP=0.1333$ and $FN=0.0949$, respectively. All these measures show the goodness of our method to ex-post dating business cycle turning points.

Turning to real-time dating of turning points, for each month, starting from February 1980 up to March 2024, we re-estimate our model from April 1959 up to that month and compute the filtered probability of recession, $\widehat{\xi}_{1,t|t}$ as given in (ref), for the last observation in the considered sample. So our first prediction is for February 1980. This is the same approach as Urga_Wang_2022_WP with two main differences. First, our indicator of recessions is very stable meaning that most of the times our indicator is equal either 0 or 1 and a thresholding procedure is seldom needed. Second, we do not use a sub-set of the $N$ series but include all of them. In Table (ref), we report the time delay of our method in detecting turning points as defined by the NBER recession indicator $REC_t$. We compare our results with those reported by Urga_Wang_2022_WP. A negative delay means that we anticipate the turning point. Our method predicts well the starting of recessions sometimes with a smaller delay than its competitors, while it tends to underestimate their duration, thus anticipating the end of recessions and resulting in a negative delay in predicting expansions.

Inflation indexes

In the last application, we consider a panel of $N=142 $ U.S. disaggregated Personal Consumption Expenditure (PCE) price monthly inflation rates from February 1959 to December 2023, for a total of $T=779$ time periods. The dataset is built as described in AL, who analyze the same data by means of a time-varying linear dynamic factor model allowing for both short and long memory dynamics. They show evidence of a structural change in the mid/end-1980s or even mid-1990s, depending on the size of the moving window considered; using the hallinliska07 information criterion, they find evidence of one factor before and after the change-point.

In Section (ref) we discussed that the model in (ref) admits the same equivalent linear representation as a model with one change point. We then apply the algorithm detailed in Section (ref) with two regimes and one common factor in each regime, namely $r_1=r_2=r=1$. Note that, in this application, it is crucial to allow for heteroskedastic idiosyncratic components, namely with $\bm\Sigma_{e1}\ne \bm\Sigma_{e2}$, as assumed in the general specification of our model in (ref): in this case, idiosyncratic components on average account for about 80% of the total variation in the data. The EM algorithm converges in 10 iterations.

The estimate of the matrix of conditional probabilities $\mathbf P$ in ((ref)) is equal to

The estimated unconditional probabilities are $\bar{\widehat{\xi}}_{1|T}=0.3770$ and $\bar{\widehat{\xi}}_{2|T}=0.6230$.\footnote{The analytical formulas in (ref) give $\bar{\widehat{\xi}}_{1|T}=0.4154$ and $\bar{\widehat{\xi}}_{2|T}=0.5846$.} By just looking at these numbers, it may seem hard to interpret the two regimes. However, by plotting $\widehat{\xi}_{1,t|T}$ and $\widehat{\xi}_{2,t|T}$ as in Figure (ref), we immediately see that, from March 1996 onwards, regime $j=2$ occurs with probability one in all time periods. Therefore, this regime can be identified with the most recent part of the sample. On the other hand, in the first part of the sample regime $j=1$ is often the most likely to occur. This finding is consistent with the results in AL: they show that the first part of the sample, in which regime $j=1$ is more likely to happen, is characterized by periods of high volatility and long memory, namely by persistent dynamics; conversely, the second part of the sample, which corresponds to regime $j=2$, is characterized by low volatility and short memory, namely by fast mean reversion. More generally, this shows that our model can also be used as a starting point to model stochastic breaks in large dimensional factor models, in the spirit of Chib_1998_JoE.

Concluding remarks

This paper develops estimation and inferential theory for high dimensional factor models with discrete regime changes in the loadings driven by a latent first order Markov process. Our estimator employs a EM algorithm based on a modified version of the Baum-Lindgren-Hamilton-Kim filter and smoother. Remarkably, the estimator does not need knowledge of the number of factors in either states. It only requires the true number of factors in the equivalent linear representation, which can be estimated using existing techniques. We derive convergence rates and asymptotic distributions of the estimators for factors and loadings, and we show their good finite sample performance through an extensive set of Monte Carlo experiments. Finally, we empirically validate our methodology through three applications to large U.S. datasets of stock returns, macroeconomic variables, and inflation indexes.

Our work can be extended along several dimensions. Two are worth mentioning. Our model allows for two regimes and the case of multiple states to capture richer dynamics is worth exploring. The challenging task of making inference on the number of regimes is also worth considering. These extensions are part of our ongoing research agenda and will be studied in future work.

{ {{.2cm} }}