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.
88,374 characters · 9 sections · 93 citation commands
Estimation of large approximate dynamic matrix factor models based on the EM algorithm and Kalman filtering
Matrix-variate time series data are becoming increasingly popular in economics and finance. For example, when forecasting regional specific economic activity chernis2020three, investigating the dynamics of international trade flows chen2022modeling, measuring of financial connectedness billio2021matrix. This has stimulated the development of high-dimensional methods to analyze matrix time series data, including matrix autoregressive models chen2021autoregressive,hsu2021matrix,billio2023bayesian, matrix panel regression models kapetanios2021estimation, and matrix factor models wang2019factor,yu2021projected,chen2021statistical,xu2024quasi,yu2024dynamic.
In this paper, we study a matrix factor model for a $p_1\times p_2$ zero-mean matrix-valued stationary process $\left\{\mathbf{Y}_t\right\}$, with latent factors following a Matrix Autoregressive (MAR) model of order $P$, i.e., for $t\in \mathbb Z$,
In (ref), $\mathbf{F}_t$ is a $k_1 \times k_2$ matrix of latent factors with $k_1,k_2<\min(p_1,p_2)$, $\mathbf{R}$ and $\mathbf{C}$ are $p_1\times k_1$ and $p_2\times k_2$ matrices of unknown row and column loadings, $\mathbf{E}_t$ is a $p_1\times p_2$ matrix of idiosyncratic components with $p_1\times p_1$ row covariance matrix $\mathbf{H}$ and $p_2\times p_2$ column covariance matrix $\mathbf{K}$. In (ref), $\mathbf{A}_\ell$ and $\mathbf{B}_\ell$, $\ell=1,\dots,P$, are $k_1\times k_1$ and $k_2\times k_2$ matrices of autoregressive parameters, and $\textbf{U}_t$ is a $k_1\times k_2$ matrix of innovations from a matrix-variate distribution with $k_1\times k_1$ row covariance matrix $\textbf{P}$ and $k_2\times k_2$ column covariance matrix $\textbf{Q}$. The processes $\left\{\mathbf{F}_t\right\}$ and $\left\{\mathbf{E}_t\right\}$ are assumed to be uncorrelated (at all leads and lags).
In our setting the idiosyncratic components are allowed to be correlated both across rows and columns, i.e., in general $\mathbf{H}$ and $\mathbf{K}$ are allowed to be full matrices, and we say that the factor model is {\it approximate}. Furthermore, the model is {\it dynamic} since the factors are autocorrelated as specified by the MAR in (ref), and, moreover, we also allow the idiosyncratic components to be autocorrelated, although no explicit model for their dynamics is introduced. Therefore, we call a model defined by (ref)-(ref) an approximate dynamic matrix factor model (DMFM). It combines the matrix factor model, as formulated by yu2021projected and chen2021statistical, and the MAR proposed by chen2021autoregressive. A DMFM has also been considered by yu2024dynamic (see Section (ref) for a detailed comparison with our work).
In this paper, we propose a new estimator of the factor loading matrices and factor matrices of the DMFM, implemented via the Expectation Maximization (EM) algorithm jointly with the Kalman smoother. We prove consistency of the spaces spanned by the estimated loadings and by the factors as $\min(p_1,p_2,T)\to\infty$.
We argue that accounting for factor dynamics via the Kalman smoother, thus considering joint estimation of all parameters and the factors, is particularly convenient as it allows the user to impose {\it a priori} restrictions on the models' parameters and/or dynamics, construct counterfactual scenarios, conditional forecasts, obtain now-casts, and deal with missing values due to different sampling frequencies or plain unavailability of the data (see, e.g., the applications in banbura2014maximum and banbura2015conditional in the case of vector time series). Furthermore, we also show that, thanks to the use of the Kalman filter, this approach is also particularly convenient to handle the case in which the data is driven by common stochastic trends, i.e., when (some of) the factors are $I(1)$ barigozzi2023measuring.
Similarly to the vector case, the DMFM can be identified only in the limit $p_1,p_2\to\infty$ due to its approximate structure. That is to say that the numbers of factors $k_1$ and $k_2$ can be consistently estimated only when both dimensions grow large. This is what allows one to disentangle the factor driven component from the idiosyncratic one. However, this forces us to work in a high-dimensional setting. This makes joint Maximum Likelihood estimation of all the parameters and the factors in (ref)-(ref) a hard if not unfeasible task due to the large number of parameters we need to estimate, which is $O( (p_1^2+p_2^2)T)$ (all autocovariances of the factors and idiosyncratic components), and due to the lack of a closed form solution.
The estimation approach we consider has two main features which allow us to solve both problems. First, it is based on a mis-specified likelihood where the idiosyncratic components are treated as if they were uncorrelated. This reduces the number of parameters to be estimated to $O( p_1+p_2)$. Second, it is an iterative approach where, in a first step, for given parameters we estimate the factors via the Kalman smoother, and, in a second step, for given factors we estimate all parameters by maximizing the expected likelihood conditional on the factors. This allows us to derive a closed form expression for all estimators.
Our approach is the generalization of the approach proposed by doz2012quasi for the vector case. However, such generalization is non-trivial, indeed, in the present matrix time series setting, we need, at each iteration of the EM algorithm, to jointly estimate the two matrices of loadings, $\mathbf R$ and $\mathbf C$, which depend on each other (the same goes for the row and column idiosyncratic covariances, $\mathbf H$ and $\mathbf K$, the MAR coefficients, $\mathbf A$ and $\mathbf B$, and the MAR innovation covariances, $\mathbf P$ and $\mathbf Q$). Respecting this bilinear structure requires modifying the algorithm accordingly and makes the derivation of the asymptotic properties more challenging.
Finally, we show the potential of the proposed approach through two applications. First, we analyze a matrix times series containing various volatility proxies for many stocks. Since not all proxies are available for all stocks, we show how to adapt the EM algorithm to deal with missing values and we then produce volatility forecasts for all stocks. Second, we analyze a matrix of time series of real macroeconomic variables of various Euro Area countries, which are clearly driven by few common trends.
The rest of the paper is organized as follows. Section (ref) discusses related works. Section (ref) presents the estimator obtained via the EM algorithm. Section (ref) presents the assumptions and the consistency results. Sections (ref) and (ref) explain how to extend the EM algorithm in presence of missing data and/or common stochastic trends. Section (ref) studies the finite sample properties of the EM algorithm through Monte Carlo simulations. Section (ref) presents two real data applications on variance proxies of financial assets and on macroeconomic indicators of the Euro Area. Appendix (ref) contains all notation, as well as relevant results on matrix operations. Appendix (ref) contains details on the EM updates. Appendix (ref) contains all proofs; Appendix (ref) explains how to identify $I(1)$ and $I(0)$ factors in the case of $I(1)$ data; Appendix (ref) contains additional simulation results.
There exist many works considering estimation only of the matrix factor model in (ref), thus without explicitly accounting for the factors' dynamics. First, wang2019factor introduce the class of large matrix factor models under the assumption of serially uncorrelated idiosyncratic components, and propose to estimate the loadings by means of eigenvectors of a long-run covariance matrix (see also chen2020constrained). Second, yu2021projected and chen2021statistical extend this approach to the case of possibly autocorrelated idiosyncratic components, and propose two different generalizations to the matrix setting of the Principal Component (PC) estimators typically used in the vector case. Both these work consider also methods for determining the number of factors (see also he2023one, and han2022rank, for alternative methods). In a similar setting, gao2023denoising consider estimating in the case of idiosyncratic components containing weak signals. Last, yuan2023two and xu2024quasi consider QML estimation of two different specifications of a matrix factor model.
To the best of our knowledge only yu2024dynamic consider a DMFM as specified by (ref)-(ref). However, our work differs in several aspects. First, we consider joint estimation of factors and parameters of the model, while they consider a two-step approach where first the loadings and the factors are estimated and then a MAR is estimated on the factors. Second, we allow the idiosyncratic components to be serially correlated, while they impose a different factor structure with time independent idiosyncratic components. Third, we derive the asymptotic properties of the factors estimated via the Kalman smoother, while they do not study the asymptotic properties of such estimator, although entertaining the possibility of retrieving the factors via filtering. As a last difference, we also study our estimation approach in presence of arbitrary patterns of missing data or stochastic trends.
Our work is also related to three other strands of the literature. First, the idea of considering a misspecified likelihood in factor analysis to make its maximization more treatable dates back to tipping1999probabilistic who, in a vector context, treated the idiosyncratic components as i.i.d.. This idea was then extended by doz2012quasi and bai2016maximum to the case of high-dimensional vectors of time series having serially and cross-sectionally correlated idiosyncratic components. In particular, doz2012quasi explicitly model the factors dynamics.
Second, there exist many factor model approaches for handling missing values in high-dimensional vector time series. On the one hand, banbura2014maximum propose an EM-based approach which we generalize to the matrix setting in this paper. On the other hand, there are a few approaches based on various modifications of standard PC analysis, see, e.g., the recent works by xiong2023large and cahan2023factor. Finally, cen2024tensor consider a PC based approach for the tensor case, which includes the matrix case.
Third, in the case of $I(1)$ vector time series, estimation of factor models via PC has been studied in a few works either under the assumption of stationary idiosyncratic components, which can be serially uncorrelated zhang2019identifying or autocorrelated bai2004estimating, or when allowing for $I(1)$ idiosyncratic components bai2004panic,barigozzi2021large. Recently, chen2025inference considered estimation via PC methods for matrix time series with $I(1)$ and $I(0)$ factors and stationary idiosyncratic components.
\paragraph{The log-likelihood.} Let consider a DMFM as defined in (ref)-(ref), and without loss of generality assume that the MAR is of order $P=1$. For a $p_1\times p_2$ matrix-valued covariance stationary process $\left\{\mathbf{Y}_t\right\}$ our data generating process is then given by:
where $\mathbf{R}$ is a $p_1\times k_1$ matrix of row loadings, $\mathbf{C}$ is a $p_2\times k_2$ matrix of column, $\mathbf{F}_t$ is a $k_1\times k_2$ matrix of latent factor, $\mathbf{E}_t$ is a $p_1\times p_2$ matrix of idiosyncratic components with covariances $\mathbf{H}$ and $\mathbf{K}$, $\mathbf A$ and $\mathbf B$ are both $k_1\times k_2$ matrices of MAR coefficients, and $\textbf{U}_t$ is a $k_1\times k_2$ matrix of innovations with covariances $\mathbf{P}$ and $\mathbf{Q}$. As usual in factor models, for simplicity and without loss of generality, we assume $\mathbb{E}\left[\mathbf{F}_t\right]=\mathbf 0_{k_1,k_2}$ and $\mathbb{E}\left[\mathbf{E}_t\right]=\mathbf 0_{p_1,p_2}$. Therefore, model (ref)-(ref) implies that $\mathbb{E}\left[\mathbf{Y}_t\right]=\mathbf 0_{p_1,p_2}$ in other words, we implicitly assume for simplicity to be working with centered data.
Denote as $\mathsf{y}_t=\textrm{vec}\left(\mathbf{Y}_t\right)$, $\mathsf{f}_t=\textrm{vec}\left(\mathbf{F}_t\right)$, $\mathsf{e}_t=\textrm{vec}\left(\mathbf{E}_t\right)$ and $\mathsf{u}_t=\textrm{vec}\left(\mathbf{U}_t\right)$, the vectorized versions of the matrices $\mathbf{Y}_t$, $\mathbf{F}_t$, $\mathbf{E}_{t}$ and $\mathbf{U}_t$, respectively. Then, consider a sample of $T$ observations, and let $\mathsf{Y}_T=\left(\mathsf{y}'_1 \cdots \mathsf{y}'_T\right)'$ and $\mathsf{E}_T=\left(\mathsf{e}'_1 \cdots \mathsf{e}'_T\right)'$ be $(p_1p_2T)$-dimensional vectors containing all the observations and idiosyncratic components, respectively, and let $\mathsf{F}_T=\left(\mathsf{f}'_1 \cdots \mathsf{f}'_T\right)'$ be the $(k_1k_2T)$-dimensional vector of factors. Let $\boldsymbol{\Omega}^{\mathsf{Y}}_T=\mathbb{E}\left[\mathsf{Y}_T\mathsf{Y}'_T\right]$, $\boldsymbol{\Omega}^{\mathsf{E}}_T=\mathbb{E}\left[\mathsf{E}_T\mathsf{E}'_T\right]$, $\boldsymbol{\Omega}^{\mathsf{F}}_T=\mathbb{E}\left[\mathsf{F}_T\mathsf{F}'_T\right]$ be covariance matrices containing all the cross-sectional row and column covariances and all the autocovariances up to lag $(T-1)$. Notice that $\boldsymbol{\Omega}^{\mathsf{F}}_T$ is fully characterized by the matrices of MAR parameter $\mathbf{A}$, $\mathbf{B}$, $\mathbf{P}$ and $\mathbf{Q}$, thus, hereafter, we denote it as $\boldsymbol{\Omega}^{\mathsf{F}}_{T(\mathbf A,\mathbf B,\mathbf P,\mathbf Q)}$.
It follows that the DMFM is fully characterized by the covariance matrix of $\mathsf{Y}_T$, which must be such that $ \boldsymbol{\Omega}^{\mathsf{Y}}_T = (\mathbb{I}_T\otimes \mathbf{C} \otimes \mathbf{R})\boldsymbol{\Omega}^{\mathsf{F}}_{T(\mathbf A,\mathbf B,\mathbf P,\mathbf Q)} (\mathbb{I}_T\otimes \mathbf{C} \otimes \mathbf{R})' + \boldsymbol{\Omega}^{\mathsf{E}}_T. $ In an approximate DMFM as the one we consider, $\boldsymbol{\Omega}^{\mathsf{E}}_T$ is allowed to be a full-matrix, but this implies that it has $\frac{p_1p_2T(p_1p_2T+1)}2$ entries to be estimated while we have only $p_1p_2T $ observations. This makes Maximum Likelihood estimation unfeasible.
A solution consists in considering a misspecified likelihood where the idiosyncratic components $\mathbf{E}_t$ are treated as if they were serially and cross-sectionally uncorrelated, i.e., when we replace $\boldsymbol{\Omega}^{\mathsf{E}}_T$ with $\mathbb{I}_T\otimes \mathrm{dg}\left(\mathbf{K}\right) \otimes \mathrm{dg}\left(\mathbf{H}\right)$. This is the approach followed in the vector factor model case, by, e.g., bai2016maximum and doz2012quasi. Under this misspecification, the vector of parameters to be estimated reduces to $ \boldsymbol{\theta}=\left(\textrm{vec}\left(\mathbf{R}\right)',\textrm{vec}\left(\mathbf{C}\right)',\textrm{vec}\left(\mathrm{dg}\left(\mathbf{H}\right)\right)',\textrm{vec}\left(\mathrm{dg}\left(\mathbf{K}\right)\right)',\textrm{vec}\left(\mathbf{A}\right)',\textrm{vec}\left(\mathbf{B}\right)',\textrm{vec}\left(\mathbf{P}\right)',\textrm{vec}\left(\mathbf{Q}\right)'\right)', $ which has dimension now growing as $p_1+p_2$, thus it can be estimated using $p_1p_2T $ observations. Consequently, we consider a misspecified, or quasi, log-likelihood given by:
Due to the introduced misspecifications, we say that the maximizer of (ref) is a QML estimator.
In principle, QML estimation of $\boldsymbol{\theta}$ can be performed by writing the model in vectorized form and maximizing the prediction error decomposition of the Gaussian likelihood obtained from the Kalman filter (see e.g. section 7.2 in durbin2012time). This approach is applicable when $p_1$ and $p_2$ are relatively small, but it becomes quickly unfeasible in larger settings due to a lack of closed form solution. Furthermore, by vectorizing the data we lose the bilinear structure of the model. We resort to the EM algorithm instead.
\paragraph{EM algorithm.} The EM algorithm is an iterative procedure proposed by dempster1977maximum to maximize the log-likelihood in problems where missing or latent observations make the likelihood intractable. This procedure works in two steps: given a set of parameter values, the E-step computes the expectation of the log-likelihood conditional on the observed data, thus “filling” the missing observations; the M-step maximizes the expected log-likelihood with respect to the model parameters. These two steps are iterated until a convergence criterion is satisfied. As the factor process $\{\mathbf{F}_t\}$ is unobserved in our setting, the EM algorithm is a suitable option to perform QML estimation.
In general, wu1983convergence proves that, when considering a Gaussian quasi-likelihood the EM algorithm converges to one of its maxima. As such, doz2012quasi consider their EM approach as Quasi Maximum Likelihood (QML) estimation of a dynamic vector factor model, a conjecture then proved by barigozzi2019quasi. For this reasons in this section we refer to our estimator as a QML estimator. However, to formally prove that the estimator defined below is effectively achieving QML estimation we would need more assumptions on the distribution of the data and the identification of the loadings space, which, in this paper, we refrain to make. For this reason, here we do not prove such equivalence but we limit to notice that, by construction, the considered log-likelihood is effectively increasing at each iteration (see the numerical results in Section (ref)).
\paragraph{Kalman smoother.} For any iteration $n\geq 0$ of the EM algorithm, and given an estimator of the parameters $\widehat{\boldsymbol{\theta}}^{(n)}$, we run the Kalman smoother on a vectorized version of the DMFM (ref)-(ref). This gives as an estimator of $\mathsf{f}_t=\textrm{vec}\left(\mathbf{F}_t\right)$ the linear projection $\mathsf{f}^{(n)}_{t|T} = \mathrm{Proj}_{\widehat{\boldsymbol{\theta}}^{(n)}}\left[\mathsf{f}_t|\mathsf{Y}_{T}\right]$ and the associated MSE, denoted as $\mathbf{\Pi}^{(n)}_{t|T}$. Moreover, by considering the Kalman smoother for the augmented state vector $(\mathsf{f}_t^\prime \;\; \mathsf{f}_{t-1}^\prime)^\prime$, we denote the top-left $k_1k_2\times k_1k_2$ block of the associated $2k_1k_2\times 2k_1k_2$ MSE as $\mathbf{\Delta}^{(n)}_{t|T}$ (see, e.g., Section 4.4 in durbin2012time, for the explicit expressions of these quantities). In particular, to run the Kalman smoother we need first to run the Kalman filter, which, in turn, requires an estimate of the inverse idiosyncratic covariance matrix. This is a hard taks in high-dimensions, but here, consistently with the misspecified log-likelihood (ref), we always consider an estimator of the misspecified diagonal covariance matrix $\text{dg}(\mathbf K)\otimes \text{dg}(\mathbf H)$, which is always invertible.
Note that the engineering literature proposes matrix versions of the Kalman filter for matrix state-space models like the one in (ref)-(ref) (e.g. choukroun2006kalman). However, these approaches heavily rely on the $\textrm{vec}\left(\cdot\right)$ operator and offer only minor computational advantages, primarily due to algebraic simplifications. Similarly to our approach, also yu2024dynamic utilize the vectorized Kalman filter.
Under joint Gaussianity of $\mathsf F_T$ and $\mathsf Y_T$, it is known that $\mathsf{f}^{(n)}_{t|T}$, $\mathbf{\Pi}^{(n)}_{t|T}$, and $\mathbf{\Delta}^{(n)}_{t|T}$ are estimators of the first and second conditional moments of $\mathsf{f}_t$ given $\mathsf Y_T$, obtained when computing expectations using the estimated parameters $\widehat{\boldsymbol{\theta}}^{(n)}$. As mentioned above, here we do not make any Gaussianity assumption. Nevertheless, we show that in the present high-dimensional setting the Kalman smoother delivers consistent estimates of the factors, thus providing a good approximation (see the results in Section (ref)). \paragraph{E-step.} In the E-step, we use the output of the Kalman smoother to compute the expected quasi log-likelihood of the approximate DMFM. Spcifically, given $\mathsf{Y}_T$ and $\widehat{\boldsymbol{\theta}}^{(n)}$, by Bayes' rule, we have\footnote{Notice that conditioning on $\mathsf{Y}_T$, which is a vector, is equivalent to conditioning on the sequence of matrices $\{\mathbf Y_1,\ldots,\mathbf Y_T\}$, hence, we can write the argument of the log-likelihood in both ways. The same applies for $\mathsf{F}_T$ and the sequence of matrices $\{\mathbf F_1,\ldots,\mathbf F_T\}$. Therefore, in order to avoid introducing further notation, hereafter, we use only the vector notation $\mathsf{Y}_T$ and $\mathsf{F}_T$ to indicate the conditioning random variables and the arguments of the log-likelihoods, even when the latter are expressed explicitly as function of matrix valued time series.} \[ \ell\left(\mathsf{Y}_T; \boldsymbol{\theta}\right) = \underbrace{\mathbb{E}_{\widehat{\boldsymbol{\theta}}^{(n)}}\left[\ell\left(\mathsf{Y}_T |\mathsf{F}_T; \boldsymbol{\theta}\right) |\mathsf{Y}_T \right] + \mathbb{E}_{\widehat{\boldsymbol{\theta}}^{(n)}}\left[\ell\left(\mathsf{F}_T; \boldsymbol{\theta}\right) |\mathsf{Y}_T \right]}_{\mathcal{Q}(\boldsymbol{\theta}, \widehat{\boldsymbol{\theta}}^{(n)})} - \mathbb{E}_{\widehat{\boldsymbol{\theta}}^{(n)}}\left[\ell\left(\mathsf{F}_T |\mathsf{Y}_T; \boldsymbol{\theta}\right)|\mathsf{Y}_T\right]. \]
As proved in dempster1977maximum, maximizing $\ell\left(\mathsf{Y}_T; \boldsymbol{\theta}\right)$ is equivalent to maximizing $\mathcal{Q}(\boldsymbol{\theta}, \widehat{\boldsymbol{\theta}}^{(n)})$, and we thus need to compute only the latter. Specifically, we have (see Appendix (ref) for the derivation)
Notice that these log-likelihoods depend directly on the data in its matrix form.
\paragraph{M-step.} In the M-step, we maximize (ref) and (ref) to obtain a new estimate of the parameters $\widehat{\boldsymbol{\theta}}^{(n+1)}$. In particular, at any $n\geq 0$ iteration, the row and column loadings estimators are given by (see Appendix (ref) for the derivation):
Clearly, the estimators of $\mathbf R$ and $\mathbf C$ depend on each other. Here, we choose to first estimate $\mathbf R$ conditional on the previous iteration estimator of $\mathbf C$. Since, as shown in the next section, all these estimators are consistent at any iteration $n\ge 0$, provided we correctly initialize the EM algorithm, we ensure, in this way, that the bilinear structure of the DMFM is preserved. We can of course equivalently choose to first estimate $\mathbf C$ conditional on the previous iteration estimator of $\mathbf R$.
By using $\widehat{\mathbf{R}}^{(n+1)}$ and $\widehat{\mathbf{C}}^{(n+1)}$ we can compute the estimators of ${\mathbf{H}}$ and ${\mathbf{K}}$ which are enforced to be diagonal matrices in agreement with the considered misspecified log-likelihood given in (ref). Thus, for $i=1,\ldots, p_1$, we have:
and $[\widehat{\mathbf H}^{(n+1)}]_{ij}=0$ if $i\ne j$. Likewise, for $i=1,\ldots, p_2$, we have:
and $[\widehat{\mathbf K}^{(n+1)}]_{ij}=0$ if $i\ne j$. As for the loadings, the estimators of $\mathbf H$ and $\mathbf K$ depend on each other, and, in order to preserve the bilinear structure of the model, we first estimate $\mathbf H$ conditional on the previous iteration estimator of $\mathbf K$.
Finally, while the individual MAR matrices are part of the model's structure, we opt for a more streamlined approach by estimating their Kronecker product directly. This choice simplifies implementation, particularly since the Kalman smoother in our algorithm is applied to vectorized data. Thus, at each iteration we compute the estimators $\widehat{\mathbf{B}\otimes\mathbf{A}}^{(n+1)}$ and $\widehat{\mathbf{Q}\otimes \mathbf{P}}^{(n+1)}$ (see Appendix (ref) for their expressions and the expressions of the alternative estimators $\widehat{\mathbf{A}}^{(n+1)}$, $\widehat{\mathbf{B}}^{(n+1)}$, $\widehat{\mathbf{P}}^{(n+1)}$, and $\widehat{\mathbf{Q}}^{(n+1)}$). Clearly, such estimators do not satisfy the constraints imposed by the bilinear structure of the MAR. Nevertheless, the asymptotic properties of the estimated loadings and factor matrices are unaffected by this choice. This is also confirmed by our simulations in Appendix (ref).
\paragraph{Initialization.} We use the projected estimator (PE) of yu2021projected to obtain initial estimates $\widehat{\mathbf{R}}^{(0)} $, $\widehat{\mathbf{C}}^{(0)}$, and $\widetilde{\mathbf {F}}_t$ of the row and column loadings and the factor matrices. Initial estimates of the idiosyncratic variances, i.e., the diagonals of $\widehat{\mathbf{K}}^{(0)}$, $\widehat{\mathbf{H}}^{(0)}$ can be obtained by computing the sample variances of the PE residual idiosyncratic components. Last, in agreement with the M-step, pre-estimators of the MAR parameters can be computed without imposing the bilinear structure, i.e., by fitting a VAR on $\widetilde{\mathsf{f}}_t\equiv \text{vec}(\widetilde{\mathbf {F}}_t)$, thus giving $\widehat{\mathbf{B\otimes\mathbf{A}}}^{(0)} $ and $\widehat{\mathbf{Q}\otimes\mathbf{P}}^{(0)}$. Expressions of all these pre-estimators are in Appendix (ref).
Finally, we initialize the Kalman filter by setting $\mathsf{f}^{(0)}_{0|0} = \mathbf{0}_{k_1k_2}$ at $n=0$ and $\mathsf{f}^{(n)}_{0|0} = \mathsf{f}^{(n-1)}_{0|T}$ at $n\geq1$, and $\mathbf{\Pi}^{(n)}_{0|0} = \mathbb{I}_{k_1k_2}$ at $n\geq0$.
\paragraph{Convergence.} As in doz2012quasi, we run the EM algorithm for a finite pre-specified number of iterations $n_{\textrm{\tiny max}}$. For a given tolerance level $\epsilon$, the algorithm is stopped at the first iteration $n^*<n_{\textrm{\tiny max}}$ such that $ \Delta\mathcal{L}_{n^*}={|\mathcal{L}(\mathsf{Y}_t,\widehat{\boldsymbol{\theta}}^{(n^*+1)})-\mathcal{L}(\mathsf{Y}_t;\widehat{\boldsymbol{\theta}}^{(n^*)})|}/{\frac{1}{2}|\mathcal{L}(\mathsf{Y}_T;\widehat{\boldsymbol{\theta}}^{(n^*+1)}) +\mathcal{L}(\mathsf{Y}_t;\widehat{\boldsymbol{\theta}}^{(n^*)})|} < \epsilon, $ where $ \mathcal{L}(\mathsf{Y}_T;\boldsymbol{\theta})$ is the one-step-ahead prediction error log-likelihood, computed via the Kalman filter.
\paragraph{Final estimators.} Once the EM algorithm reaches convergence, we define the EM estimator of the model parameters as $\widehat{\boldsymbol{\theta}}\equiv \widehat{\boldsymbol{\theta}}^{(n^*+1)}$. In particular, the estimated factor loadings are given by $\widehat{\mathbf{R}}\equiv\widehat{\mathbf{R}}^{(n^*+1)}$ and $\widehat{\mathbf{C}}\equiv\widehat{\mathbf{C}}^{(n^*+1)}$. Finally, we obtain a final estimate of the factor matrices by running the Kalman smoother one last time, that is $\widehat{\mathbf{F}}_t \equiv \text{unvec}({\mathsf{f}^{(n^*+1)}_{t|T}})$ for any $t=1,\dots,T$.
We make the following assumptions on the loadings and the factors.
Assumptions (ref)-(ref) matches Assumptions B and C in yu2021projected. Assumption (ref) guarantees the stationarity of the MAR(1) model. Assumption (ref) requires the factor innovations $\{\mathbf{U}_t\}$ to have positive definite covariance matrix.
We characterize the idiosyncratic component through the following assumption.
Assumptions (ref) closely matches Assumptions A and D in yu2021projected. In particular, Assumption (ref) controls serial idiosyncratic dependence by requiring them to be $\alpha$-mixing (see also chen2021statistical). Assumption (ref) imposes finite absolute fourth moments and requires $\mathbf{H}$ and $\mathbf{K}$ to be positive definite matrices, and, since $\mathbf{H}$ and $\mathbf{K}$ are only determined up to a positive constant, and only their Kronecker product, $\mathbf{K} \otimes \mathbf{H}$, is uniquely defined, we impose the constraint $\textbf{tr}(\mathbf{H}) = p_1$ (see, e.g., viroli2012matrix). Assumption (ref) controls cross-sectional idiosyncratic dependence across both rows and columns (see also chen2021statistical), while Assumption (ref) bounds fourth order cumulants to allow for consistent estimation of the second order moments.
Finally, the dependence between common and idiosyncratic components is controlled through the following assumption, which matches Assumption E in yu2021projected.
Under the above assumptions we can then derive theoretical results on the convergence rates of the EM estimators for the loading and factor matrices $\widehat{\mathbf{R}}$, $\widehat{\mathbf{C}}$, and $\widehat{\mathbf{F}}_t$, defined in Section (ref).
These rates are comparable with those of the PE, which we use to initialize the EM algorithm. In particular, from yu2021projected, the PE are such that, as $\min\left\{p_1,p_2,T\right\} \rightarrow \infty$,
Consider, for example, the error we have for the initial estimator ${\widehat{\mathbf{R}}^{(0)}}$. While the first term $\sqrt {T p_1}$ is the same as the one we find in Proposition (ref), we also have a slower comparable term $\sqrt{Tp_2}$ coming from the initial estimator ${\widehat{\mathbf{C}}^{(0)}}$ of the columns. This is due to the fact that at each iteration we need also to estimate the idiosyncratic variances, which require estimating both the row and the column loadings first. The price to be paid is however negligible or even null when $p_1$ and $p_2$ are of the same order of magnitude.
The rate in yu2021projected for the PE is instead $\min(\sqrt T p_1,\sqrt T p_2, \sqrt{p_1p_2})$. The first term $\sqrt{p_1p_2}$ is the same and corresponds to the case of known parameters, while the other rates, due to the estimation of the parameters, are slower because they inherit the slower rates of the EM estimator of the loadings loadings. Nevertheless, if ${p_2}/{T}\to 0$ and ${p_1}/{T}\to 0$ as $\min(p_1,p_2,T)\to\infty$, then the Kalman smoother and the PE have the same $\sqrt{p_1p_2}$ rate, which would be obtained by vectorizing the data and estimating the factors via projection onto the true loadings.
We conclude by noting that many applications in economics and finance may involve settings in which one dimension of the matrix $\mathbf Y_t$, say $p_1$, together with the sample size $T$, diverges ($p_1, T \to \infty$), while the other dimension, say $p_2$, remains fixed ($p_2 < \infty$). Alternatively, both dimensions may diverge ($p_1, p_2 \to \infty$) while the sample size is fixed. An example of the former case arises in asset pricing applications, where the number of traded assets can be large, whereas the number of liquidity or volatility proxies is typically finite. An example of the latter case is firm-level balance sheet data, where observations are available for many firms and variables but only at low frequency (e.g., annually). Our estimator remains consistent under both scenarios.
If the data contains missing values, the estimation of the factors and their second moments with the Kalman smoother is still possible (see, e.g., durbin2012time for details). Now, since $\mathbb{E}_{\widehat{\boldsymbol{\theta}}_{(n)}} \left[\ell(\mathsf{F}_T;\boldsymbol{\theta}) | \mathsf{Y}_T \right]$ depends only on the factors but not on the data, its expression remains unchanged. Thus, the estimators of $\mathbf A\otimes \mathbf B$ and $\mathbf P\otimes \mathbf Q$ are also unchanged. However, $\mathbb{E}_{\widehat{\boldsymbol{\theta}}_{(n)}} \left[\ell(\mathsf{Y}_T|\mathsf{F}_T;\boldsymbol{\theta}) | \mathsf{Y}_T \right]$ depends on the data and therefore its expression is affected by missing values. It follows that we need to adjust the estimators of $\mathbf{R}$, $\mathbf{C}$, $\mathbf{H}$, and $\mathbf{K}$ in the M-step accordingly. To this end, we extend the procedure of banbura2014maximum to the matrix setting.
Let $\mathbf{W}_t$ be a $p_1 \times p_2$ matrix with $(i,j)$ entry equal to zero if $y_{tij}$ is missing and equal to one otherwise. For any iteration $n\geq 0$, the estimators of the row and column loadings are then modified to (see Appendix (ref) for the derivation of these expressions): \[ {0mu} {0mu} {0mu}
\]
Since $\mathbf{Y}_t$ contains missing observations, the initialization procedure described in Section (ref) cannot be applied directly. To address this, we introduce an additional preliminary step in which the missing entries are imputed using a suitable imputation method. To this end, an obvious choice consists in applying the imputation procedure for tensor factor models proposed by cen2024tensor.
Let $w_{t,i,j}$ be the entry $(i,j)$ of $\mathbf W_t$ and let $\eta$ be the fraction of missing data, i.e., such that $(1-\eta)\le \min_{i=1,\ldots, p_1}\min_{i=j,\ldots, p_2}T^{-1}{\sum_{t=1}^T w_{t,i,j}}$. Then, from cen2024tensor, we see that, under our assumptions, we have initial estimators such that, as $\min\left\{p_1,p_2,T\right\} \rightarrow \infty$,
These results extend the findings of xiong2023large from the vector to the matrix setting. They are comparable to the results for the initial PC-based estimator studied by yu2021projected and chen2021statistical. However, because no projection step is involved in the method of cen2024tensor, the rates in (ref)-(ref) are not directly comparable to those for the PE analyzed in yu2021projected.
From the discussion after Proposition (ref) it is clear that, when dealing with missing data and initializing the EM algorithm with the estimator by cen2024tensor, the consistency rates for our estimated loadings will be the minimum of the rates in (ref) and (ref). Specifically, the loadings estimated via the EM algorithm are such that, as $\min\left\{p_1,p_2,T\right\} \rightarrow \infty$,
The same rates hold for the single row or column estimators $\widehat{\mathbf r}_{i}$, $i=1,\ldots, p_1$, and $\widehat {\mathbf c}_{j}$, $j=1,\ldots, p_2$.
By the same arguments, we expect the Kalman smoother computed using the EM estimator of the parameters to be a consistent estimator of the factors with a rate given by the minimum of the rates in (ref) and (ref) and $\sqrt{p_1p_2}$ corresponding to the rate for known parameters. Hence, as $\min\left\{p_1,p_2,T\right\} \rightarrow \infty$,
which matches the rate in cen2024tensor.
A formal proof of the statements (ref), (ref), and (ref) would follow verbatim the same steps of the proofs of Propositions (ref) and (ref), respectively, but when using as initial estimators those satisfying (ref)-(ref) instead of the PE which satisfy (ref)-(ref). Hence, such proof is omitted.
In this section, we consider the case in which $\{\text{vec}(\mathbf Y_t)\}$ is no more covariance stationary but is instead an $I(1)$ process, following a matrix factor model as in (ref), but under the assumption that $\{\text{vec}(\mathbf F_t)\}$ is $I(1)$ and $\{\text{vec}(\mathbf E_t)\}$ is $I(0)$.
In a macroeconomic context, it is reasonable to assume that the elements of $\mathbf F_t$ are driven both by common trends, which are $I(1)$, and stationary components, which we could consider as common cycles (see, e.g., the applications in Section (ref) and in barigozzi2023measuring). In such a case, there exist a $k_1\times k_1$ invertible matrix $\bm{\mathcal R}$ and a $k_2\times k_2$ invertible matrix $\bm{\mathcal C}$ such that our model can be rewritten as (see Appendix (ref) for the explicit expressions)
where $\{\text{vec}(\mathbf G_{1t})\}$ is an $I(1)$ process with $\mathbf G_{1t}$ being $r_1\times r_2$, $\mathbf R_1$ being $p_1\times r_1$ and $\mathbf C_1$ being $p_2\times r_2$, while $\{\text{vec}(\mathbf G_{0t})\}$ and $\{\text{vec}(\mathbf E_t)\}$ are $I(0)$, with $\mathbf G_{0t}$ being $q_1\times q_2$, $\mathbf R_0$ being $p_1\times q_1$ and $\mathbf C_0$ being $p_2\times q_2$, so that $k_1=r_1+q_1$ and $k_2=r_2+q_2$.
The model on the rightmost side of (ref) is introduced by chen2025inference who propose PC-type estimators of both $\mathbf G_{1t}$ and $\mathbf G_{0t}$. Here, instead we focus on the estimation of the common component, i.e., $\mathbf S_t=\mathbf R\mathbf F_t\mathbf C'$, which does not require identification of the common trends. If, according to (ref), we make the assumption that $\{\text{vec}(\mathbf F_t)\}$ is a cointegrated process driven by $r_1r_2$ common trends, then the correct specification for its dynamics is either via a VECM or a VAR in levels. Hence, the DMFM must be estimated by applying the EM algorithm and the Kalman smoother on the levels of the data, i.e., without differencing them in order to achieve stationarity.
Under the assumption that $\{\text{vec}(\mathbf E_t)\}$ is stationary, we can still adopt the same initialization as in the stationary case, i.e., we can still use the PE as described in Section (ref). From chen2025inference we see that, under our assumptions plus the assumption of cointegrated factors, we have initial estimators such that, as $\min\left\{p_1,p_2,T\right\} \rightarrow \infty$,
These results generalize to the matrix case the results by bai2004estimating for the vector case.
Once again, our results can then be directly adapted to this setting. From the discussion after Proposition (ref) it is clear that the consistency rates for our estimated loadings will be the minimum of the rates in (ref) and (ref), i.e., the loadings estimated via the EM algorithm are such that, as $\min\left\{p_1,p_2,T\right\} \rightarrow \infty$,
The same rates hold for the single row or column estimators $\widehat{\mathbf r}_{i}$, $i=1,\ldots, p_1$, and $\widehat {\mathbf c}_{j}$, $j=1,\ldots, p_2$. A formal proof of the statements (ref) and (ref) would follow verbatim the same steps of the proofs of Propositions (ref) and (ref), respectively, but when using as initial estimators those satisfying (ref)-(ref) instead of the PE which satisfy (ref)-(ref). Hence, such proof is omitted.
By the same arguments, we expect the Kalman smoother computed using the EM estimator of the parameters to be a consistent estimator of the factors with rate the minimum between the rates in (ref) and (ref), divided by $\sqrt T$ due to non-stationarity, and $\sqrt{p_1p_2}$ corresponding to the rate for known parameters. Hence, as $\min\left\{p_1,p_2,T\right\} \rightarrow \infty$,
which matches the rate in chen2025inference. A formal proof of this statement is, however, less straightforward and, thus, it should be regarded just as an informed conjecture.
We conclude with three remarks. First, in the presence of missing observations, the EM algorithm can still be applied using the update modifications discussed in Section (ref). Since no iputation method exists for the non-stationary case, we propose to initialize the algorithm by running the EM procedure on a fully observed subset of the original matrix $\mathbf{Y}_t$. Simulation results in Appendix (ref) confirm the effectiveness of this approach.
Second, the number of common trends can be determined by following the same approach proposed in chen2025inference and based on eigenvalue ratios of suitable second moment matrices.
Third, if all or some of the idiosyncratic components were non-stationary due to the presence of stochastic trends, then, the above approach would not be consistent. Indeed, in that case the loadings should be estimated from the differenced data as explained in bai2004panic and barigozzi2021large in the vector case. In this case, the estimated loadings would retain the same rates as the PE for stationary data given in (ref) and (ref). Moreover, the Kalman smoother should be run by adding as additional latent states all those idiosyncratic components which are $I(1)$, in a way similar to the approach proposed by banbura2014maximum for serially correlated, but stationary, idiosyncratic components. This case is left for further research.
We perform Monte Carlo simulations in order to assess the finite sample properties of the proposed EM estimator and the Kalman smoother. For $t=1,\dots,T$, we generate observations according to the following DMFM: \[
\] where $\mathfrak{D}_{k_1,k_2}(\mathbf{0}_{k_1,k_2}, \mathbf{P}, \mathbf{Q}) $ and $\mathfrak{D}_{p_1,p_2}(\mathbf{0}_{p_1,p_2}, \mathbf{H}, \mathbf{K})$ denote general matrix distributions of dimensions $k_1\times k_2$ and $p_1\times p_2$, centered on zero, and with covariance matrices $\mathbf{P}, \mathbf{Q}$ and $\mathbf{H}, \mathbf{K}$, respectively. We consider $\mathfrak{D}$ either to be a matrix normal (N) or a matrix skew-t (St) distribution with 4 degrees of freedom. The loading matrices are such that $\left[\mathbf{R}\right]_{ij},\left[\mathbf{C}\right]_{ij} \sim \mathcal{U}(-1,1)$. The matrix of latent factors follows a MAR(1) process with $\mathbf{B}=\mu \frac{\mathbf{B}^*}{\lvert \nu^{(1)}\left(\mathbf{B}^{*} \otimes \mathbf{A}\right) \rvert}$ where $[\textbf{B}^*]_{ii}, [\textbf{A}]_{ii} \sim\mathcal{U}(0.7,0.9)$ and $[\textbf{B}^*]_{ij}, [\textbf{A}]_{ij} \sim \mathcal{U}(0,0.5)$ for $i\ne j$. Note that $\mu$ defines the maximum eigenvalue of the matrix $\mathbf{B}\otimes\mathbf{A}$ allowing us to control whether the matrix factor process is stationary or not. In particular when $\mu=1$, the simulated factors are driven by one common $I(1)$ trend. Throughout, we set $k_1=2$ and $k_2=2$.
The idiosyncratic components follow a MAR(1) process with \[ \left[D\right]_{ij},\left[G\right]_{ij} =
\qquad \left[H\right]_{ij}, \left[K\right]_{ij} =
\] with $\tau$ and $\delta$ controlling the degree of cross-sectional and serial correlation, respectively.
For each performance measure considered, we report its average and standard deviation over 100 replications. We use the column space distance $\mathcal{D}(\mathbf{R}, \widehat{\mathbf{R}})$ and $\mathcal{D}(\mathbf{C}, \widehat{\mathbf{C}})$ to evaluate the loadings matrices estimators, which, for any $m\times n$ matrix $\mathbf{A}$, is defined as $ \mathcal{D}(\mathbf{A}, \widehat{\mathbf{A}}) = \big\lVert \widehat{\mathbf{A}} \left(\widehat{\mathbf{A}}'\widehat{\mathbf{A}}\right)^{-1}\widehat{\mathbf{A}}'-\mathbf{A}\left(\mathbf{A}'\mathbf{A}\right)^{-1}\mathbf{A}' \big\rVert. $ We also consider the mean squared error in recovering the signal $\mathbf{S}_t=\mathbf{R}\mathbf{F}_t\mathbf{C}'_t$, defined as $ \textrm{MSE}_{\textbf{S}} = (Tp_1p_2)^{-1} \sum_{t=1}^{T} \lVert \widehat{\mathbf{S}}_t - \mathbf{S}_t \rVert^2_{\textrm{F}}, $ where $\widehat{\mathbf{S}}_t= \widehat{\mathbf{R}}\widehat{\mathbf{F}}_t\widehat{\mathbf{C}}$ denotes the estimated signal, as described in Sections (ref) or (ref).
In Table (ref) we compare the performance of the EM estimators for the loading and factor matrices with those of the PE both in the stationary and the $I(1)$ cases. The EM algorithm improves upon PE across all the different settings. Furthermore, in Figure (ref) we show for one replication the log-likelihood as function of the number of iterations of the EM algorithm. As expected the log-likelihood increases monotonically and the first few iterations seem to be the most important ones.
We then introduce missing observations in the data generating process. After simulating the matrix $\mathbf Y_t$ with no missing values as described above we introduce two patterns of missing observations widely seen in empirical application: (i) randomly missing, i.e., removing at each point in time observations of $\mathbf Y_t$ at random with a constant probability $\pi=\{25\%,50\%\}$; (ii) block missing, i.e., when a fixed portion $\pi=\{25\%,50\%\}$ of $\mathbf{Y}_t$ is removed for a given period of time. For the case of block missing we remove the bottom-right quarter and the right-half of $\mathbf{Y}_t$ for the first half of the time series when $\pi=25\%$ and $\pi=50\%$, respectively.
In this case, besides $\mathcal{D}(\mathbf{R}, \widehat{\mathbf{R}})$, $\mathcal{D}(\mathbf{C}, \widehat{\mathbf{C}})$, and $\textrm{MSE}_{\textbf{S}}$, we also investigate the goods of our imputation method by computing $ \textrm{MSE}_{\textbf{Y}^{(0)}} =(Tp_1p_2)^{-1} \sum_{t=1}^{T} \lVert (\widehat{\mathbf{S}}_t - \mathbf{Y}_t)\circ (\mathbf{1}_{p_1,p_2}-\mathbf{W}_t) \rVert^2_{\textrm{F}}$, where $\widehat{\mathbf{S}}_t= \widehat{\mathbf{R}}\widehat{\mathbf{F}}_t\widehat{\mathbf{C}}$ denotes the estimated signal, as described in Section (ref), and $\mathbf{W}_t$ is the binary matrix indicating observed entries.
Following the discussion in Section (ref), we adopt the imputation method proposed by cen2024tensor to fill in missing values prior to initializing the EM algorithm. Because this method requires stationarity, we restrict the analysis to stationary settings. Table (ref) reports summary statistics for the relative performance of the EM estimator compared to the PE estimator applied to the imputed data. The results indicate that the EM algorithm yields systematically improved estimates over PE. Additional simulation results based on initialization using a balanced subpanel are in Appendix (ref).
\paragraph{Forecasting volatilities.}
Despite the abundant use of high-frequency data in the financial econometrics literature, their availability is often limited to major equity indices or large U.S. stocks bollerslev2018risk, limiting the chance of building high-frequency-based estimates of volatility for a large number of traded companies. Given that volatility measures tend to covary across assets barigozzi2016generalized, a natural question is whether high-frequency-based volatility measures on a set of assets can be used to improve volatility estimates for a set of assets for which only daily observations are available.
We collect daily returns and realized measures for $30$ assets listed in the S&P500 under the Financial GICS sector. The data covers the period that goes from the beginning of 2006 to the end of 2010, covering the Great Financial Crisis. We consider $10$ realized measures of the daily integrated volatility. In particular, we have 7 high-frequency measures based on intra-daily data\footnote{These are: 5-min and 15-min realized variance, autocorrelation-corrected 5-min realized variance hansen2006realized, realized range christensen2007realized, realized kernel barndorff2008designing, pre-averaged realized variance jacod2009microstructure, maximum likelihood realized variance xiu2010quasi.} and three low-frequency proxies based on the opening (O), highest (H), lowest (L), and closing (C) daily prices (OHLC hereafter).\footnote{These are: the daily range $(H - L)^2/(4\log 2)$, the O/C adjusted daily range $0.5(H - L)^2 - (2\log 2 - 1)(C - O)^2$, and the O/C adjusted daily range $(H -C)(H -O)+(L-C)(L-O)$.} These measures are available only for half of the stocks in the sample, as we have access to high-frequency data solely for those assets. For the remaining stocks we only have daily data, and can therefore compute just the three OHLC variance proxies. We thus obtain a matrix time series of $p_2=10$ daily variance proxies on $p_1=30$ assets for $T=1259$ days, with a block of missing observations corresponding to 35% of the total number of possible observations which is $p_1p_2T$.
Our data can be modeled as a 2-layers hierarchical factor model which in turn is equivalent to a matrix factor model. First, let $\sigma^2_{i,t}$ be the $t$th day latent variance of the $i$th asset and define $\boldsymbol{\widetilde{\sigma}}^2_{t}=(\widetilde{\sigma}^2_{1,t},\dots,\widetilde{\sigma}^2_{p_1,t})'$, with $\widetilde{\sigma}^2_{i,t}=\sigma^2_{i,t}/\bar{\sigma}^2_{i,t}$ and $\overline{\sigma}^2_{i}=(\prod^T_{t=1}\sigma^2_{i,t})^{{1}/{T}}$, for all $i=1,\dots,p_1$. We assume that the vector of centered latent log-variances for all assets, $\log(\boldsymbol{\widetilde{\sigma}}^2_{t})$, follows a factor model with $\mathbf{f}_{t}$ being a vector of $k_1$ common factors, e.g., representing the stock market, that is
where $\mathbf{R}$ is a $p_1 \times k_1$ loading matrix and $ \boldsymbol{\varepsilon}_{t}$ contains the idiosyncratic component for each asset.
Second, let $\mathbf{s}_{i,t}$ be the vector of $p_2$ variance proxies for the $i$th asset on the $t$th day and define $\widetilde{\mathbf{s}}_{i,t}=(\widetilde{s}_{i,1,t},\dots,\widetilde{s}_{1,p_2,t})'$, with $\widetilde{s}_{i,j,t}=s_{i,j,t}/\bar{s}_i$ and $\overline{s}_i=(\prod^T_{t=1}\prod_{j=1}^{p_2} s_{i,j,t})^{{1}/{(Tp_2)}}$. It is reasonable to assume that the centered vector of log-variance proxies of asset $i$ follows a one factor model, where the common factor is the latent volatility $\log (\widetilde{\sigma}^2_{i,t})$ of asset $i$, that is
where $\textbf{c}$ is a $p_2$-dimensional loading vector and $\boldsymbol{\epsilon}_{i,t}$ is a $p_2$-dimensional vector contaning the measurement errors of all variance proxies of asset $i$.
It follows that the $p_2$ vector of observed centered log-transformed variance proxies for the $i$th asset follows a 2-layer factor model. Indeed, by substituting the transposed of (ref) into (ref), we have
where $\mathbf r_i'$ is the $i$th row of $\mathbf R$ and ${\varepsilon}_{i,t}$ is the $i$th element of $\boldsymbol{\varepsilon}_{t}$. By letting $\mathbf{Y}_t = \left(\log \left(\widetilde{\mathbf{s}}_{1,t}\right)', \dots, \log \left(\widetilde{\mathbf{s}}_{p_1,t}\right)' \right)'$, we see that (ref) is equivalent to the matrix factor model in (ref) with $\mathbf{E}_t= \boldsymbol {\varepsilon}_{t}\textbf{c}' + (\boldsymbol{\epsilon}_{1,t}'\cdots\boldsymbol{\epsilon}_{p_1,t}')'$. For economic reasons we fix the number of columns factors to $k_2=1$, indeed, this corresponds to the number of latent variance factor underlying all proxies. As for the number of row factors the eigenvalue-ratio criterion by cen2024tensor suggests to set $k_1=1$.
We then conduct a forecasting exercise. We define an in-sample window of $750$ observations for the models estimation and leave $509$ observations for the out-of-sample forecast evaluation. We estimate a DMFM on the in-sample window using our proposed EM algorithm modeling $f_t$, which, since $k_1=k_2=1$ is now a scalar, as an AR(1), and obtain one-step-ahead forecasts of $\widetilde{\sigma}^2_{i,t}$ as $\widehat{\widetilde{\sigma}}^2_{i,t|t-1} = \exp(\widehat{r}_{i} \, \widehat{A} \, \widehat{f}_{t-1|t-1})$, where $\widehat{r}_i$ is the estimated row loading for the $i$th asset and $\widehat{A}$ is the estimated autoregressive coefficient. For comparison, we also estimate an analogous DMFM on the in-sample window using the proposed EM algorithm, but restricted to the $15 \times 3$ sub-matrix of assets for which only low-frequency volatility measures are available, i.e., to a balanced subpanel of the considered dataset.
Table (ref) reports the out-of-sample MSE ratios comparing the model estimated on the reduced matrix to that estimated on the full dataset, along with the p-values from the diebold1995comparing test for each financial asset at the daily frequency. The out-of-sample MSE for the model estimated on the reduced matrix is higher for twelve out of fifteen assets, reaching up to 7% in some cases. According to the Diebold-Mariano test of equal predictive accuracy, these differences are statistically significant for nine assets. This finding underscores the advantage of incorporating high-frequency volatility proxies from assets that covary with those for which we only have access to low-frequency measures, and thus shows the importance of having a method which allows us to deal with panels with missing observations.
\paragraph{Macroeconomic trends in the Euro Area.} We analyze a collection of macroeconomic indicators from EA countries.\footnote{The data is available at \href{https://zenodo.org/doi/10.5281/zenodo.10514667}{https://zenodo.org/doi/10.5281/zenodo.10514667}.} Specifically, we consider 39 real macroeconomic indicators across three categories: National Accounts, Labor Market Indicators, and Industrial Production and Turnover. These indicators are collected at either monthly or quarterly frequency for eight countries: Austria, Belgium, Germany, Spain, France, Italy, the Netherlands, and Portugal, resulting in a matrix-valued time series of dimensions $(p_1, p_2) = (8, 39)$. The dataset spans the period from January 2000 to November 2024 ($T = 299$).
By applying the eigenvalue ratio criterion by yu2021projected on the differenced data we find evidence of one row factor and three column factors, but we cannot say whether any of these is $I(1)$ or stationary. To this end we can instead apply the eigenvalue ratio approach proposed by chen2025inference on the non-differenced data, showing evidence of just one $I(1)$ common factor, i.e., a common trend. Since the factor matrix is actually a 3-dimensional vector this implies that the process of latent factors is indeed cointegrated with two cointegrating relations. As explained before, and differently from chen2025inference, here we are not interested in identifying the trend or the other factors separately, but we are interested in recovering the whole common component of the data, i.e., $\widehat{\mathbf{S}}_t= \widehat{\mathbf{R}}\widehat{\mathbf{F}}_t\widehat{\mathbf{C}}$. Hence, we can apply the methodology described in Section (ref). Moreover, since the considered dataset contains both monthly and quarterly varaibles, we apply our method when also imputing missing values as described in Section (ref). Figure (ref) reports the GDP of Germany, France, Spain, and Italy (in black), which are quarterly, together with their estimated common components $\widehat{\mathbf{S}}_t$ (in red), which are monthly time series. While the GDPs of Germany and France are strongly related to the common EA factors, Spain and Italy display more idiosyncratic behavior, hinting at a different level of commonality among EA countries.
This paper introduces a methodology for estimating a large approximate DMFM using the EM algorithm combined with Kalman filtering. We establish the consistency of the spaces spanned by the estimated loadings and factors as $\min\left(p_1, p_2, T\right) \to \infty$. Our estimation framework accommodates missing observations and unit root data.
Our approach can be readily adapted to include additional constraints on the model parameters chen2020constrained or to explicitly model the dynamics of the idiosyncratic components which can be modeled as additional latent states banbura2014maximum. Moreover, the proposed approach can be further and straightforwardly extended to tensor data of higher order, enhancing its applicability to more complex data structures.