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.
74,703 characters · 20 sections · 49 citation commands
Matrix Quantile Factor Model
\def\spacingset#1{ {#1}} \spacingset{1}
\if11 \fi
\if01 {
} \fi
{\it Keywords:} Two-way factor model; Row factor space; Column factor space;Check loss function
\spacingset{1.9}
The present paper studies the matrix sequence data with a latent low-rank structure. Instead of modeling the mean functionals conditional on the latent factors as in recent works (wang2019factor; yu2022projection; Jing2021Community), we model the conditional quantiles by an interactive effect of the row and column sections. The parameters, row and column factor loadings and factor matrices, are learnt by minimizing the empirical check loss under an identifiability constraint. For the first time, we derive the statistical accuracy of the estimated factor loadings and the factors. In the modeling side, this is essentially a work where the matrix factor structure meets the quantile feature representation.
On the matrix factor structure, we assume in the present paper that each matrix observation is driven by a much lower dimensional factor matrix, and that the two cross sections along the row and column dimensions interact with each other and thus generates the entries of the $t$-th data matrix ${\mbox{\boldmath $X$}}_t=(X_{ijt})_{p_1\times p_2}$. For example, in recommending system, ${\mbox{\boldmath $X$}}_t$ is a rating matrix of $p_1$ customers and $p_2$ commodities, and the scores in ${\mbox{\boldmath $X$}}_t$ are high when the latent common consumption preferences of the customers match the latent common features of the $p_2$ goods. We focus on matrix sequence rather than large vectors appeared in standard statistics for two-fold reasons. First, many recent data sets in financial market, medical research, social networks and electronic business platform, are themselves well organized intrinsically in matrices. The matrix factor structure is empirically found in these data sets and works well in applications, c.f, Liu2023Simultaneous for function magnetic resonance imaging data and Jing2022Community for political blog network. Second, modeling the matrix-value data with a low-rank structure, e.g. model ((ref)) below, makes the model parsimonious and statistical inference efficient once the structure is interpretably reasonable. A naive approach to analyze the data matrix ${\mbox{\boldmath $X$}}_t$ is to “flatten" it into a long vector $vec({\mbox{\boldmath $X$}}_t)$ by piling down column by column or row by row. After that, existing vector factor models in stock2002forecasting, bai2002determining, Trapani2018A, Barigozzi2020Sequential, fan2013large, kong2017number, kong2018systematic, kong2019factor and Chen2021quantile, can be applied. However, the flattened vector factor modeling easily misses the interplay between the row and column sections, and has parameter complexity of order $O(p_1p_2+T)$ while the row-column interaction model (see ((ref)) below) of order $O(p_1+p_2+T)$. This is also where the efficiency gain of the present paper comes from compared with the vector quantile factor modeling. For more motivation to study matrix or tensor sequence data, we refer to recent interesting works: wang2019factor, chen2023testing, Chang2023Modelling, He2023Iterative,Yuan2023Two-way, zhang2024tucker and ZhangLiuGuoYuenWelsh2024JASA-modeling.
On the quantile feature representation, mathematically, a $\tau$-quantile for a random variable $Y$ is $Q_{\tau}(Y)=\inf\{y; P(Y\leq y)\geq \tau\}$. With increasing complexity of data sets, how to understand the co-movement of the quantiles of large-dimensional random vectors evolving in time is of vital importance in theory and applications. To the best of our knowledge, Chen2021quantile is the first paper that models the $\tau$ quantile of a large vector by a {\sc vector} factor structure. Ando2021Quantile extended Chen2021quantile to allow for observed covariates in modeling the panel quantiles. But, so far, no works are done to investigate the co-movement of the quanitles of a matrix sequence or even more generally tensor sequence. The more parsimonious interactive quantile factor representation, compared to the vector quantile factor model, is still not well understood in achieving higher statistical estimation precision. There is yet a challenge in establishing the second-order asymptotic theory for the estimated loadings, for example, the technique for the vector quantile factor model in the existing works can not be trivially extended to the matrix sequences. There is no guarantee that the Hessian matrix of the empirical check loss function penalized by the commonly used identifiability constraints is locally positive definite around the true parameters, making a second-order expansion of the penalized loss function difficult and hence the difficulty of the central limit theorems of the estimated row and column factor loadings.
In this paper, we estimate the row and column factor loadings and factors by minimizing the empirical check loss function under constraint. Our theory demonstrates that our estimates converge at rate $O_p((\min\{p_1T,p_2T,p_1p_2\})^{-1/2})$ in the sense of averaged Frobenius norm, if the quantile interactive mechanism is effective. Our theoretical rate is faster than $O_p((\min\{p_1p_2,T\})^{-1/2})$, the rate expected from the vector quantile factor analysis by vectorizing ${\mbox{\boldmath $X$}}_t$, which is more pronounced when the sequence length $T$ is short. To the best of our knowledge, this is the first result on the estimation of the matrix quantile factor model and reveal of the interactive effect in reducing the estimation error. Our theory also shows that the convergence rates are reached without any moment constraints on the idiosyncratic errors, hence robust to the heavy tails of the heterogeneous idiosyncratic errors. To derive the central limit theorems of the smoothed versions of the loading estimates, we introduce an augmented Lagrangian function that not only takes the identification rotation constraints but also a cleverly constructed extra term into consideration. We have proved the reversibility of Hessian matrix for the smoothed quantile loss function with penalty locally around the true parameters, which is crucial to obtain the stochastic expansion of smoothed estimates and hence the central limit theorems. To solve the minimization of the check loss function under constraint, we present an iterative algorithm to find an approximate solution. Extensive simulation studies show that the numerical solutions are close enough to the true parameters, and demonstrate the robustness to the heavy tails. To determine the pair of the row and column factor numbers, we present three criteria, which are proved to be consistent and verified by simulations.
The present paper is organized as follows. Section (ref) gives the matrix quantile factor model and the estimation method. Main results on estimating the cross-sectional factor spaces and set-up Assumptions are provided in Section (ref). Section (ref) presents three model selection criteria to determine the numbers of row and column factors. Section (ref) presents a smoothed version of the loading estimates and the central limit theorems. Section (ref) conducts simulations and Section (ref) presents an empirical data analysis. Section (ref) concludes. The technical proofs are relegated to the supplementary material.
We model the co-movement of the quantiles of all entries in each matrix by the following matrix quantile factor model.
where ${\mbox{\boldmath $R$}}_{\tau}$, ${\mbox{\boldmath $C$}}_{\tau}$ and ${\mbox{\boldmath $F$}}_{t,\tau}$ are the $p_1\times k_{1,\tau}$ row factor loading matrix, $p_2\times k_{2,\tau}$ column factor loading matrix and $k_{1,\tau}\times k_{2,\tau}$ common factor matrix, respectively, and ${\mbox{\boldmath $E$}}_{t,\tau}$ is an error matrix. Obviously, $Q_{\tau}({\mbox{\boldmath $E$}}_{t,\tau}|{\mbox{\boldmath $F$}}_{t,\tau})=0$. The subscript $\tau$ emphasizes the dependence on $\tau$. That being said, the low-rank quantile structure is heterogeneous across different quantile levels, as seen in our real data analysis. Model ((ref)) demonstrates that the entries of ${\mbox{\boldmath $X$}}_t$ depends on how close the rows of ${\mbox{\boldmath $R$}}_{\tau}$ are to the rows of ${\mbox{\boldmath $C$}}_{\tau}$, i.e, an interactive effect between the row and column sections of variables. We refer to ${\mbox{\boldmath $R$}}_{\tau}{\mbox{\boldmath $F$}}_{t,\tau}{\mbox{\boldmath $C$}}_{\tau}^{\prime}$ and ${\mbox{\boldmath $E$}}_{t,\tau}$ as the common and idiosyncratic components, respectively. Model ((ref)) includes the two-way quantile fixed effect model as a special example. In particular, setting ${\mbox{\boldmath $R$}}_{\tau}=({\mbox{\boldmath $\alpha$}}_{p_1\times 1}(\tau), (\tilde{{\mbox{\boldmath $R$}}}_{\tau})_{p_1\times (k_{1,\tau}-2)}(\tau), {\bf{1}}_{p_1\times 1})$, ${\mbox{\boldmath $F$}}_{t,\tau}=\mbox{diag}\{1, (\tilde{{\mbox{\boldmath $F$}}}_{t,\tau})_{(k_{1,\tau}-2)\times (k_{2,\tau}-2)}, 1\}$ and ${\mbox{\boldmath $C$}}_{\tau}=({\bf{1}}_{p_2\times 1}, (\tilde{{\mbox{\boldmath $C$}}}_{\tau})_{p_2\times (k_{2,\tau}-2)}, {\mbox{\boldmath $\beta$}}_{p_2\times 1}(\tau))$, $$ {\mbox{\boldmath $X$}}_t={\mbox{\boldmath $\alpha$}}(\tau) {\bf{1}}_{1\times p_2}+{\bf{1}}_{p_1\times 1}{\mbox{\boldmath $\beta$}}'(\tau)+\tilde{{\mbox{\boldmath $R$}}}_{\tau}\tilde{{\mbox{\boldmath $F$}}}_{t,\tau}\tilde{{\mbox{\boldmath $C$}}}_{\tau}^{\prime}+{\mbox{\boldmath $E$}}_{t,\tau}, $$ where ${\mbox{\boldmath $\alpha$}}(\tau)$ and ${\mbox{\boldmath $\beta$}}(\tau)$ represent the time-invariant quantile fixed effects along the row and column dimensions, respectively. They can be heterogeneous across the rows and/or columns.
While the vector factor model is conceptually a generative mechanism for a single cross-section of variables that are closely related in nature, the matrix factor model in ((ref)) is a two-way joint generative modeling in two totally different cross-section of variables. Though different in interpretations, model ((ref)) can be mathematically rewritten in the form of a vector factor model
where $\mbox{vec}(\cdot)$ is the vectorization operator that stacks the columns of a matrix into a long vector and $\otimes$ stands for the Kronecker product operator. A general vector factor model for an observed vector ${\mbox{\boldmath $x$}}_t$ is typically expressed as
where ${\mbox{\boldmath $L$}}$, ${\mbox{\boldmath $f$}}_t$ and ${\mbox{\boldmath $\epsilon$}}_t$ are the loading matrix, factor vector and idiosyncratic error vector, respectively. That is, ((ref)) can be mathematically regarded as a vector factor model with parameter restrictions ${\mbox{\boldmath $L$}}={\mbox{\boldmath $C$}}_{\tau}\otimes {\mbox{\boldmath $R$}}_{\tau}$, $p=p_1p_2$ and $k=k_1k_2$. When the Kronecker structure ${\mbox{\boldmath $C$}}_{\tau}\otimes {\mbox{\boldmath $R$}}_{\tau}$ is latent in the matrix sequence, a simple vectorization and vector principal component analysis would yield consistent estimate of the factor loading matrix ${\mbox{\boldmath $L$}}$ (and hence ${\mbox{\boldmath $C$}}_{\tau}\otimes {\mbox{\boldmath $R$}}_{\tau}$) up to orthogonal transformation in the sense of averaged Frobenius norm. Expected from the vector quantile factor analysis in Chen2021quantile, the consistent rate for estimating ${\mbox{\boldmath $L$}}$ is $(\min\{p_1p_2,T\})^{-1/2}$. To recover the row and column factor spaces spanned by ${\mbox{\boldmath $R$}}_{\tau}$ and ${\mbox{\boldmath $C$}}_{\tau}$, a further nearest Kronecker decomposition has to be done, c.f., van2000kronecker, but the resulting estimates of ${\mbox{\boldmath $R$}}_{\tau}$ and ${\mbox{\boldmath $C$}}_{\tau}$ depend on the estimation error for ${\mbox{\boldmath $L$}}$. The other way around with vector quantile factor analysis is to minimize the empirical check loss function by restricting ${\mbox{\boldmath $L$}}={\mbox{\boldmath $C$}}_{\tau}\otimes {\mbox{\boldmath $R$}}_{\tau}$, but the number of restrictions is diverging which leads to complex computation. The matrix form ((ref)) gives a neat joint modeling of a two-way structure to start from.
Coming back to the general model ((ref)), the row factor loading matrix ${\mbox{\boldmath $R$}}_{\tau}$, the column factor loading matrix ${\mbox{\boldmath $C$}}_{\tau}$ and the factor matrix ${\mbox{\boldmath $F$}}_{t, \tau}$ are not separately identifiable, though the common component itself is under some signal conditions. Indeed, there exists orthonormal square matrices ${\mbox{\boldmath $O$}}_R$ and ${\mbox{\boldmath $O$}}_C$, such that ${\mbox{\boldmath $R$}}_{\tau}{\mbox{\boldmath $F$}}_{t, \tau}{\mbox{\boldmath $C$}}_{\tau}^{\prime}={\mbox{\boldmath $R$}}^*_{\tau}{\mbox{\boldmath $F$}}^*_{t, \tau}{\mbox{\boldmath $C$}}^{*\prime}_{\tau}$ where ${\mbox{\boldmath $R$}}^*_{\tau}={\mbox{\boldmath $R$}}_{\tau}{\mbox{\boldmath $O$}}_R$, ${\mbox{\boldmath $C$}}^*_{\tau}={\mbox{\boldmath $C$}}_{\tau}{\mbox{\boldmath $O$}}_C$ and ${\mbox{\boldmath $F$}}_{t,\tau}^*={\mbox{\boldmath $O$}}_R^{\prime}{\mbox{\boldmath $F$}}_{t,\tau}{\mbox{\boldmath $O$}}_C^{\prime}$. Without loss of generality, we assume throughout the paper that
To estimate the parameters, we propose to minimize the empirical check loss function subject to ((ref))
with respect to $\theta=\{\mathbf{r}_1,..., \mathbf{r}_{p_1};\mathbf{c}_1,...,\mathbf{c}_{p_2};{\mbox{\boldmath $F$}}_1,...,{\mbox{\boldmath $F$}}_T\}$, where $\rho_{\tau}(u)=(\tau-I\{u\leq 0\})u$, and $\mathbf{r}_i^{\prime}$ and $\mathbf{c}_j^{\prime}$ are the $i$-th row of ${\mbox{\boldmath $R$}}_{\tau}$ and $j$-th row of ${\mbox{\boldmath $C$}}_{\tau}$, respectively. Our estimates, denoted by $\hat{{\mbox{\boldmath $R$}}}_{\tau}$, $\hat{{\mbox{\boldmath $F$}}}_{t,\tau}$ and $\hat{{\mbox{\boldmath $C$}}}_{\tau}$, are simply the minimizers of the above empirical check loss function assuming that $k_{1, \tau}$ and $k_{2, \tau}$ are known numbers of factors a priori. Later, we will give consistent estimates of $k_{1, \tau}$ and $k_{2, \tau}$ by three methods. Notice that the empirical check loss function is not a convex function jointly in ${\mbox{\boldmath $R$}}_{\tau}$, ${\mbox{\boldmath $F$}}_{t, \tau}$ and ${\mbox{\boldmath $C$}}_{\tau}$, but it is a marginally convex function when the other two are fixed. Hence, we propose to optimize it via an iterative algorithm; see Algorithm (ref) below.
Although $\mathbb{M}_{p_1p_2T}(\theta)$ is not a joint convex function, it is convex in each iteration in one component of $(\hat{{\mbox{\boldmath $R$}}}(h), \hat{{\mbox{\boldmath $C$}}}(h), \hat{{\mbox{\boldmath $F$}}}_t(h))$ with the other two given. Motivated by Ge2017No, we set the initial values in Algorithm (ref) by random initialization. Our simulation shows that the algorithm converges fast and leads to accurate estimation.
In this section, we present a main result on the estimation accuracy of the estimated row and column factor loading matrices. Before stating the theorem, we give some technical assumptions. Without confusion, we suppress the dependence on $\tau$ of the notation $k_{1,\tau}$ and $k_{2,\tau}$, and write them simply as $k_1$ and $k_2$.
Assumption (ref)-1 is standard in the literature, e.g., the compactness of the parameters were assumed in Chen2021quantile, and the existence of the limits in ((ref)) and ((ref)) is guaranteed by the law of large numbers under various weak-correlation conditions. Assumption (ref)-2 assumed the existence of density functions which are uniformly bounded from below in compact sets, see also similar conditions in Chen2021quantile. Assumption (ref)-3 restricts that the idiosyncratic errors are conditionally independent but maybe dependent unconditionally, see the same condition in Chen2021quantile. Even if ((ref)) is satisfied, the columns of loading matrices ${\mbox{\boldmath $R$}}_{\tau}$ and ${\mbox{\boldmath $C$}}_{\tau}$ are identifiable only up to a positive or negative sign. We henceforth make a convention that the first nonzero entry of each column of ${\mbox{\boldmath $R$}}_{\tau}$ and ${\mbox{\boldmath $C$}}_{\tau}$ is positive.
Theorem (ref) demonstrates that the plug-in estimate $\hat{{\mbox{\boldmath $C$}}}_{\tau}\otimes\hat{{\mbox{\boldmath $R$}}}_{\tau}$ has convergence rate:
Expected from Chen2021quantile, the convergence rate of estimating the loading space spanned by ${\mbox{\boldmath $C$}}_{0,\tau}\otimes{\mbox{\boldmath $R$}}_{0,\tau}$ under the framework of the vector quantile factor model ((ref)) is $O_p((\min\{p_1p_2,T\})^{-1/2})$ by piling down the columns of each observed matrix into a long vector. A simple comparison shows that the latter rate is no faster than ours, and in particular, when $p_1p_2$ dominates $T$, ours is strictly faster than the rate by vectorizing the matrix. This is intuitively interpretable because the structure restriction of ${\mbox{\boldmath $L$}}$ in ((ref)) is not observed in the vector quantile analysis.
We propose three different methods to select the numbers of factors. The first selects the numbers of factors by rank minimization (RM), the second uses the information criterion (IC), while the third implements the eigenvalue ratio thresholding approach (ER). As before, the dependence on $\tau$ in all mathematical notations are suppressed for simplicity.
Let $K_{1}$ and $K_{2}$ be two positive integers larger than $k_{1}$ and $k_{2}$, respectively. Let $\mathcal{A}^{K_{1}}$ be compact subset of $\mathbb{R}^{K_{1}}$, $\mathcal{F}^{K_{1}\times K_{2}}$ be compact subset of $\mathbb{R}^{K_{1}\times K_{2}}$ and $\mathcal{B}^{K_{2}}$ be compact subset of $\mathbb{R}^{K_{2}}$. Assume that
for all $t$, $(\mathbf{r}_{0i}, \bf{0})\in \mathcal{A}^{K_1}$ and $(\mathbf{c}_{0j}, \bf{0})\in\mathcal{B}^{K_2}$. Let $\mathbf{r}_{i}^{K_{1}} \in \mathbb{R}^{K_{1}}$, ${\mbox{\boldmath $F$}}_{t}^{K_{1}\times K_{2}} \in \mathbb{R}^{K_{1}\times K_{2}}$, $\mathbf{c}_{j}^{K_{2}} \in \mathbb{R}^{K_{2}}$ for all $i, t, j$ and write
Consider the following normalization,
Define
and
Moreover, write $\widehat{{\mbox{\boldmath $F$}}}^{K_{1}\times K_{2}}=(\widehat{{\mbox{\boldmath $F$}}}^{K_{1}\times K_{2}}_{1},\cdots,\widehat{{\mbox{\boldmath $F$}}}^{K_{1}\times K_{2}}_{T})$ and
The rank minimization estimator of the numbers of factors, $k_{1}$ and $k_{2}$, are defined as
where $C_{p_{1}p_{2}T}$ is a sequence that goes to 0 as $p_{1},p_{2},T\rightarrow \infty$. That being said, $\widehat{k}_{1}^r$ and $\widehat{k}_{2}^r$ are, respectively, the numbers of the diagonal elements of $$ \sum_{t=1}^{T}\widehat{{\mbox{\boldmath $F$}}}^{K_{1}\times K_{2}}_{t}(\widehat{{\mbox{\boldmath $F$}}}^{K_{1}\times K_{2}}_{t})'/T,\quad \sum_{t=1}^{T}(\widehat{{\mbox{\boldmath $F$}}}^{K_{1}\times K_{2}}_{t})'\widehat{{\mbox{\boldmath $F$}}}^{K_{1}\times K_{2}}_{t}/T $$ that are larger than the threshold $C_{p_{1}p_{2}T}$.
The second estimator of $(k_{1}, k_{2})$ is similar to the IC-based estimator of bai2002determining, but is adaptive to the matrix observation and the check loss function. For $(l_1, l_2)\in \mathcal{P}=\{0,..., K_1\}\times \{0,..., K_2\}$, we search the minimizer of a penalized empirical check loss function.
The IC-based estimator of $(k_{1}, k_{2})$ is defined as
where $\hat{\mathbf{\theta}}^{l_1l_2}$ is similarly defined as $\hat{\mathbf{\theta}}^{K_1K_2}$ except for replacing $(K_1, K_2)$ by $(l_1, l_2)$, pretending that there are $l_1$ row factor and $l_2$ column factors.
Due to the assumption of ${\mbox{\boldmath $F$}}_{t}^{K_{1}\times K_{2}}$ in section 4.1, we expect $(\hat\sigma_{T,k_1+1}^{K_1},\ldots,\hat\sigma_{T,K_1}^{K_1})$ and $(\hat\sigma_{T,k_2+1}^{K_2},\ldots,\hat\sigma_{T,K_2}^{K_2})$ to be redundant and negligible. Therefore, motivated by the eigenvalue ratio approach in Ahn2013eigenvalue, a direct estimator for $(k_1, k_2)$ is given by $$ \hat k_1^{ER}=\arg\max_{1\le k\le K_1-1}\frac{\hat\sigma_{T,k}^{K_1}}{\hat\sigma_{T,k+1}^{K_1}+c_0L_{p_1p_2T}^{-2}}, \ \hat k_2^{ER}=\arg\max_{1\le k\le K_2-1}\frac{\hat\sigma_{T,k}^{K_2}}{\hat\sigma_{T,k+1}^{K_2}+c_0L_{p_1p_2T}^{-2}}, $$ where $c_0$ is a small positive constant so that the denominator is always larger than 0. In our simulation studies and real data analysis, we set $c_0=10^{-4}$.
The non-smoothness of the check loss function and the incidental-parameter problem make it difficult to derive the asymptotic distribution of the estimators $\widehat{\theta}$. As in the asymptotic analysis of quantile regression, one way to overcome these difficulties is to expand the expected score function and obtain a stochastic expansion for $\widehat{\mathbf{r}}_{i}-\mathbf{r}_{0i}$.
We proceed by defining a new estimator of $\mathbf{\theta}_{0}$, denoted as $\widetilde{\mathbf{\theta}}$ which relies on the following smoothed quantile optimization:
where $$ \mathbb{S}_{p_{1}p_{2}T}(\theta)=\frac{1}{p_{1}p_{2}T}\sum_{i=1}^{p_{1}}\sum_{j=1}^{p_{2}}\sum_{t=1}^{T} \Big[\tau-K\Big(\frac{X_{ijt}-\mathbf{r}_{i}'{\mbox{\boldmath $F$}}_{t}\mathbf{c}_{j}}{h}\Big)\Big](X_{ijt}-\mathbf{r}_{i}'\mathbf{F}_{t}\mathbf{c}_{j}), $$ such that $K(z)=1-\int_{-1}^{z}k(z)dz$, $k(z)$ is a continuous kernel function with support $[-1,1]$ and $h$ is a bandwidth parameter that goes to 0 as $p_{1},p_{2}$ and $T$ grow.
To derive the central limit theorem, we instead introduce an augmented Lagrangian function which is equivalent to ((ref)), that is
where
where $b_1$, $b_2$ and $b_3$ are positive Lagrangian multipliers, $\mathbf{F}_{tp.}$ is the $p$-th row of $\mathbf{F}_{t}$, $\mathbf{F}_{t.q}$ is the $q$-th column of ${\mbox{\boldmath $F$}}_{t}$ and $f_{t,pq}$ is the element of the $p$-th row and $q$-th column of ${\mbox{\boldmath $F$}}_{t}$.
The equivalence between ((ref)) and ((ref)) stems from the nonnegative property of $\mathbb{P}_{p_{1}p_{2}T}(\mathbf{\theta})$. Theoretically the minima of ((ref)) are achieved if and only if $\mathbb{P}_{p_{1}p_{2}T}(\mathbf{\theta})=0$. The first two terms $\mathbb{P}_{1}(\mathbf{\theta})$ and $\mathbb{P}_{2}(\mathbf{\theta})$ are associated with the typical four rotation constraints in ((ref)) under which the row and column factor loadings and the factor matrices are uniquely identified up to signs. The additional augmented term $\mathbb{P}_{3}(\mathbf{\theta})$, however, is carefully designed to ensure positive definiteness of the Hessian matrix of the penalized loss function, making the augmented Lagrangian function ((ref)) convex locally in $R$, $C$ and $F_t$'s around their true values. This renders a feasible second-order expansion of the estimation errors which makes the central limit theorems of the estimates easily derived. Since in the matrix factor model ((ref)), there are two cross-sections, row and column, only making use of the rotation identification constraint, i.e. $\mathbb{P}_{1}(\mathbf{\theta})+\mathbb{P}_{2}(\mathbf{\theta})$, as traditionally done in the quantile {\sc vector} factor model, is difficult to prove the local convexity of the penalized loss function. The proposed augmented Lagrangian loss function ((ref)) skillfully solved this problem, see the Appendix for the technical details.
Before stating the central limit theorem, we define, for all $i, j, t$.
The above conditions are standard in smoothed quantile optimization, with the exception of Assumption (ref)-5. Note that, as in Galvao and Kato (2016), we require $k(z)$ to be a higher-order kernel function to control the higher-order terms in the stochastic expansions of the estimators. However, Galvao and Kato (2016) assumed that $m^{-1}<c<1/3$, while we need $m^{-1}<c<1/6$ . This arises from the fact that the incidental parameters, $\mathbf{r}_{0i}$, ${\mbox{\boldmath $F$}}_{0t}$ and $\mathbf{c}_{0j}$, in quantile factor models enter the model interactively, but no interactive fixed effects appear in the panel quantile models considered by these authors.
We generate data from the following matrix series,
where ${\mbox{\boldmath $R$}}$ and ${\mbox{\boldmath $C$}}$ are $p_1\times k_1$ and $p_2\times k_2$ matrices, respectively. We set $k_1=2$ and $k_2=3$. The factor process follows an autoregressive model such that ${\mbox{\boldmath $F$}}_t=0.2{\mbox{\boldmath $F$}}_{t-1}+\Xi_t$. $g_t$ is a scalar random variable satisfying $g_t=0.2g_{t-1}+\epsilon_t$. The entries in ${\mbox{\boldmath $R$}}$, ${\mbox{\boldmath $C$}}$, $\{\Xi_t\}$ and $\{\epsilon_t\}$ are all generated from i.i.d. $\mathcal{N}(0,1)$. The entries of $\{{\mbox{\boldmath $E$}}_t\}$ are i.i.d. from $\mathcal{N}(0,1)$, or $t$ distributions with degree of freedom being 3 or 1, covering both light-tailed and heavy-tailed distributions. $\theta^*$ is a parameter controlling the signal-to-noise ratio.
To ensure the identification condition ((ref)), a normalization step should be applied to the loading and factor score matrices. For instance, when $\tau=0.5$, do singular-value decomposition to ${\mbox{\boldmath $R$}}$ and ${\mbox{\boldmath $C$}}$ as \[ {\mbox{\boldmath $R$}}={\mbox{\boldmath $U$}}_R{\mbox{\boldmath $D$}}_R{\mbox{\boldmath $V$}}_R={\mbox{\boldmath $U$}}_R{\mbox{\boldmath $Q$}}_R,\quad {\mbox{\boldmath $C$}}={\mbox{\boldmath $U$}}_C{\mbox{\boldmath $D$}}_C{\mbox{\boldmath $V$}}_C={\mbox{\boldmath $U$}}_C{\mbox{\boldmath $Q$}}_C. \] Further define \[ \tilde{\mbox{\boldmath $\Sigma$}}_1=\frac{1}{Tp_1p_2}\sum_{t=1}^T{\mbox{\boldmath $Q$}}_R{\mbox{\boldmath $F$}}_t{\mbox{\boldmath $C$}}^\prime{\mbox{\boldmath $C$}}{\mbox{\boldmath $F$}}_t^\prime{\mbox{\boldmath $Q$}}_R^\prime,\quad \tilde{\mbox{\boldmath $\Sigma$}}_2=\frac{1}{Tp_1p_2}\sum_{t=1}^T{\mbox{\boldmath $Q$}}_C{\mbox{\boldmath $F$}}_t^\prime{\mbox{\boldmath $R$}}^\prime{\mbox{\boldmath $R$}}{\mbox{\boldmath $F$}}_t{\mbox{\boldmath $Q$}}_C^\prime, \] and the eigenvalue decompositions \[ \tilde{\mbox{\boldmath $\Sigma$}}_1=\tilde{\mbox{\boldmath $\Gamma$}}_1\tilde{\mbox{\boldmath $\Lambda$}}_1\tilde{\mbox{\boldmath $\Gamma$}}_1^\prime, \quad \tilde{\mbox{\boldmath $\Sigma$}}_2=\tilde{\mbox{\boldmath $\Gamma$}}_2\tilde{\mbox{\boldmath $\Lambda$}}_2\tilde{\mbox{\boldmath $\Gamma$}}_2^\prime. \] Then, the normalized loading and factor score matrices are \[ \tilde{\mbox{\boldmath $R$}}=p_1^{1/2}{\mbox{\boldmath $U$}}_R\tilde{\mbox{\boldmath $\Gamma$}}_1,\quad \tilde{\mbox{\boldmath $C$}}=p_2^{1/2}{\mbox{\boldmath $U$}}_C\tilde{\mbox{\boldmath $\Gamma$}}_2,\quad \tilde{\mbox{\boldmath $F$}}_t=\tilde{\mbox{\boldmath $\Gamma$}}_1^\top{\mbox{\boldmath $Q$}}_R{\mbox{\boldmath $F$}}_t{\mbox{\boldmath $Q$}}_C^\prime\tilde{\mbox{\boldmath $\Gamma$}}_2. \] We are actually estimating $\tilde{\mbox{\boldmath $R$}}$, $\tilde{\mbox{\boldmath $C$}}$ and $\tilde{\mbox{\boldmath $F$}}_t$. Moreover, in the iterative algorithm, we will normalize the estimators similarly in each step, so that condition ((ref)) is always satisfied.
The simulation results for $\tau\ne 0.5$ and the results for dependent idiosyncratic errors are postponed to the Supplementary Material to save space and comply with the page limit.
This section aims to verify the effectiveness of the proposed methods for estimating the numbers of row and column factors, when $\tau=0.5$. Table (ref) reports the frequencies of exact estimation when $(T,p_1,p_2)$ grows gradually and the noises are sampled from different distributions. The approaches proposed in Chen2023Statistical and yu2022projection are taken as competitors, which are also designed for matrix factor models. Another natural idea is to first vectorize the data matrices ${\mbox{\boldmath $X$}}_t$ and then use the approach in Chen2021quantile, which expects to lead to an estimation of the total number of $k=k_1k_2$ factors in theory. $K_1, K_2$ are set as 6 for matrix factor models while $k_{\max}=12$ for Chen2021quantile's method. The $\theta^*$ is set to be 3.
Following Chen2021quantile, for rank-minimization we set $C_{p_1p_2T}=\delta L_{p_1p_2T}^{2/3}$, where $\delta=(\hat\sigma_{T,1}^{K_1}+\hat\sigma_{T,1}^{K_2})/2$. For the information criterion, we actually use an accelerated algorithm in the simulation rather than direct grid search in $\{1,...,K_1\}\times\{1,...,K_2\}$. In detail, we first fix $l_2=K_2$ and estimate $k_1$ by grid search in $\{1,...,K_1\}$. Next, we fix $l_1=\hat k_1$ and estimate $k_2$ by grid search in $\{1,...,K_2\}$. The thresholding parameter for the information criterion is set as $C_{p_1p_2T}=\delta L_{p_1p_2T}$, which is slightly smaller than that for rank-minimization.
By Table (ref), when the noises are from the standard normal distribution, the proposed three approaches with matrix quantile factor model perform comparably with the $\alpha$-PCA ($\alpha=0$) by Chen2023Statistical and the projected estimation (PE) by yu2022projection. On the other hand, when the noises are from heavy-tailed distributions $t_3$ or $t_1$, the $\alpha$-PCA and PE methods gradually lose accuracy, while the proposed three methods remain reliable, due to the robustness of check loss functions. The vectorized method doesn't work in this example mainly because the dimensions are much smaller compared with the settings in Chen2021quantile and we are considering weak signals with large $\theta^*$. Moreover, the data matrix after vectorization is severely unbalanced ($p_1p_2\gg T$), making the idiosyncratic errors matter too much.
Next, we investigate the accuracy of the estimated loadings and factor scores by different approaches. We use the similar settings as in Table (ref) and let $\tau=0.5$. Note that the minimizers to the check loss function are not unique, so the estimated loading matrices converge only after a rotation. Due to such an identification issue, we will mainly focus on the estimation accuracy of the loading spaces. Let ${\mbox{\boldmath $R$}}_0$ and $\hat {\mbox{\boldmath $R$}}$ be the true and estimated loading matrices respectively, both satisfying the identification condition in ((ref)). We define the distance between the two loading spaces by \[ \mathcal{D}({\mbox{\boldmath $R$}}_0,\hat{\mbox{\boldmath $R$}})=\bigg(1-\frac{1}{k_1p_1^2}\text{tr}(\hat{\mbox{\boldmath $R$}}^\prime{\mbox{\boldmath $R$}}_0{\mbox{\boldmath $R$}}_0^\prime\hat{\mbox{\boldmath $R$}})\bigg)^{1/2}. \] It's easy to see that $\mathcal{D}({\mbox{\boldmath $R$}}_0,\hat{\mbox{\boldmath $R$}})$ always takes value in the interval $[0,1]$. A smaller value of $\mathcal{D}({\mbox{\boldmath $R$}}_0,\hat{\mbox{\boldmath $R$}})$ indicates more accurate estimation of ${\mbox{\boldmath $R$}}_0$. When $\mathcal{D}({\mbox{\boldmath $R$}}_0,\hat{\mbox{\boldmath $R$}})=0$, the two loading spaces are exactly the same. Similar distance can be defined between ${\mbox{\boldmath $C$}}_0$ and $\hat {\mbox{\boldmath $C$}}$. Let ${\mbox{\boldmath $W$}}_0={\mbox{\boldmath $C$}}_0\otimes{\mbox{\boldmath $R$}}_0$, $p=p_1p_2$ and $\hat{\mbox{\boldmath $W$}}$ be an estimate of ${\mbox{\boldmath $W$}}_0$. Similarly, we define $$ \mathcal{D}({\mbox{\boldmath $W$}}_0,\hat{\mbox{\boldmath $W$}})=\bigg(1-\frac{1}{kp^2}\text{tr}(\hat{\mbox{\boldmath $W$}}^\prime{\mbox{\boldmath $W$}}_0{\mbox{\boldmath $W$}}_0^\prime\hat{\mbox{\boldmath $W$}})\bigg)^{1/2}. $$ The existing vector quantile factor analysis estimates ${\mbox{\boldmath $W$}}_0$ by $\hat{\mbox{\boldmath $W$}}=\hat{\mbox{\boldmath $L$}}$ given in Chen2021quantile. The matrix quantile factor analysis estimates ${\mbox{\boldmath $W$}}_0$ by the plug-in estimator $\hat{\mbox{\boldmath $W$}}=\hat{{\mbox{\boldmath $C$}}}\otimes\hat{{\mbox{\boldmath $R$}}}$.
Table (ref) reports the estimation accuracy of the loading spaces by different methods over 500 replications. The conclusions almost follow those in (ref). The estimation based on matrix quantile factor models (“mqf”) is accurate and stable under all settings, while $\alpha$-PCA and “PE” only work for light-tailed cases. Even under the normal cases, “mqf” can outperform “PE”, mainly because the latter only contains one-step iteration thus relying on a good initial projection direction. There are some enormous errors for $\alpha$-PCA, “PE” and the vectorized method in the table.
We verify the asymptotic distributions of the smoothly estimated factor loadings in Theorem (ref). We set $\theta^*=g_t=1$ in ((ref)) and generate ${\mbox{\boldmath $F$}}_t$ from i.i.d. $\mathcal{N}(0,1)$. Figure (ref) plots the empirical density of $\hat {\mbox{\boldmath $R$}}_{11}$ after standardization according to Theorem (ref), over 1000 replications with $\tau=0.5$ when the entries of ${\mbox{\boldmath $E$}}_t$ are i.i.d. from $\mathcal{N}(0,1)$, $t_3$ or $t_1$. Figure (ref) clearly shows the asymptotic normality of the estimators with well fitted variances.
In this section, we apply the proposed matrix quantile factor model and associated estimators to the analysis of a real data set, the Fama-French 100 portfolios data set. This an open resource provided by Kenneth R. French, which can be downloaded from the website \url{http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/data_library.html}. It contains monthly return series of 100 portfolios, structured into a $10\times 10$ matrix according to 10 levels of market capital size (S1-S10) and 10 levels of book-to-equity ratio (BE1-BE10). Considering the missing rate, we use the data from 1964-01 to 2021-12 in this study, covering 696 months. Similar data set has ever been studied in wang2019factor and yu2022projection. The data library also provides information on Fama-French three factors and excess market returns. Following the preprocessing steps in wang2019factor and yu2022projection, we first subtract the excess market returns and standardize each of the portfolio return series. In the first step, we provide some descriptive information of the data set in the Supplementary Material to save space.
In the second step, we fit the matrix quantile factor model. The numbers of row and column factors should be determined first. Table (ref) provides the estimated $(k_1,k_2)$ using the proposed three approaches at different quantiles $\tau$. The results by the vectorized method with Chen2021quantile are also reported in the table, which leads to the estimation of total number of factors. By Table (ref), the proposed eigenvalue ratio method and information criterion always lead to an estimate of $\hat k_1=\hat k_2=1$, while the rank minimization approach gives more row and/or column factors when $\tau\in [0.15,0.9]$. The vectorized method leads to an estimate of 2 factors in total at most quantiles. Based on the results, there should be at least one powerful row factor and column factor in the system, and potentially one weak row factor and/or column factor. When $\tau$ is at the edge, the leading factor becomes more influential. On the other hand, the approaches in wang2019factor and yu2022projection will both lead to $\hat k_1=\hat k_2=1$. In this example, it might be a good choice to use $\hat k_1=\hat k_2=2$ when $\tau$ is around $0.5$, and $\hat k_1=\hat k_2=1$ when $\tau$ is at the edge.
The next step is to estimate the loading matrices and factor scores with $\hat k_1=\hat k_2=2$. It's worth noting that the quantile factor models can handle missing values naturally by optimization only with non-missing entries, e.g., in the completely random missing case by defining \[ \mathbb{M}_{p_1p_2T}(\theta)=\frac{1}{p_1p_2T}\sum_{(i,j,t)\in \mathcal{M}}\rho_{\tau}(X_{ijt}-\mathbf{r}_i^\prime{\mbox{\boldmath $F$}}_t \mathbf{c}_j), \] where $\mathcal{M}$ indicates the index set of all non-missing entries. However, the $\alpha$-PCA and “PE” methods require to impute the missing entries first. Considering that the missing rate is small in this example ($0.23\%$), we use simple linear interpolation method to impute missing data. To measure the similarity of two estimated loading spaces, we define the following indicator: \[ S(\hat {\mbox{\boldmath $R$}}_1,\hat {\mbox{\boldmath $R$}}_2)=\frac{1}{k_1}\text{tr}\bigg(\frac{1}{p_1^2}\hat{\mbox{\boldmath $R$}}_1^\prime\hat{\mbox{\boldmath $R$}}_2\hat{\mbox{\boldmath $R$}}_2^\prime\hat{\mbox{\boldmath $R$}}_1\bigg), \] where $\hat {\mbox{\boldmath $R$}}_1$ and $\hat{\mbox{\boldmath $R$}}_2$ are the estimated $p_1\times k_1$ row loading matrices. Note that the columns of $\hat {\mbox{\boldmath $R$}}_1$ and $\hat{\mbox{\boldmath $R$}}_2$ are orthogonal after scaling. Therefore, the value $p_1^{-1}\hat{\mbox{\boldmath $R$}}_2\hat{\mbox{\boldmath $R$}}_2^\prime\hat{\mbox{\boldmath $R$}}_1$ is actually the projection matrix of $\hat{\mbox{\boldmath $R$}}_1$ to the space of $\hat{\mbox{\boldmath $R$}}_2$. The value of $S(\hat {\mbox{\boldmath $R$}}_1,\hat {\mbox{\boldmath $R$}}_2)$ will always be in the interval $[0,1]$. When the two loading spaces are closer to each other, the value of $S(\hat {\mbox{\boldmath $R$}}_1,\hat {\mbox{\boldmath $R$}}_2)$ will be larger.
The last four columns of Table (ref) report the similarity of estimated loading spaces by matrix-quantile-factor-model and two competitors, “PE” and the vectorization approach. For the vectorization, we calculate similarity by considering the Kronecker product $\hat{\mbox{\boldmath $C$}}\otimes \hat{\mbox{\boldmath $R$}}$. It's seen that the similarity indicators for the matrix-quantile-factor based approach and “PE” approach are very close to 1, implying that the estimated loading spaces are almost the same, especially when $\tau$ is near 0.5. However, when $\tau$ is at the edge, the difference of the estimated loading spaces becomes more significantly. For the vectorization approach, the estimated loading space is always not similar to that from the matrix models, consistent with our findings from the simulation study.
Now we aim to interpret the matrix quantile factors in this example. Table (ref) presents the estimated $\hat {\mbox{\boldmath $R$}}$ and $\hat{\mbox{\boldmath $C$}}$ by matrix quantile factor model at $\tau=0.5$, $0.05$, $0.95$ as well as those by “PE`”. By Table (ref), the effects of the row factors and column factors are closely related to market capital sizes and book-to-equity ratios. From the perspective of size ($\hat{\mbox{\boldmath $C$}}$), under the matrix quantile factor model with $\tau=0.5$, the small-size portfolios load more heavily on the first factor than the large-size portfolios. Moreover, the second factor has opposite effects on small-size portfolios and large-size ones. Similar results are found for the “PE” method, although the values of loadings are not exactly the same. Taking $\tau$ at edge will lead to different finding, where the first factor has more significant effect on the large-size portfolios. In other words, the edge quantile factors show disparate information of the data.
From the perspective of book-to-equity ratio ($\hat{\mbox{\boldmath $R$}}$), with $\tau=0.5$, the large-BE portfolios load more heavily on the first factor than small-BE ones, while the second factor has opposite effects on the two classes. The “PE” factors show similar trend after orthogonal transformation (changing sign). When $\tau=0.05$ and $\tau=0.95$, the first factor depends on the average of all portfolios. It's worth noting that the reported row and column factors are highly suggestive, because they coincide with financial theories. The capital size and book-to-equity ratio are known to be two important factors affecting portfolio returns in negative collaboration. The row and column factors in this example might be closely related to the SMB and HML factors in portfolio theory.
We also verify the robustness of the check losss function by the matrix quantile factor model with this real example. The results are postponed in the Supplementary Material.
By Table (ref), when $\tau$ is at the edge, the similarity indicator decreases, suggesting that considering the edge quantiles might be helpful for extracting extra information from the data. However, by Figure 3 in the supplementary material, the low similarity can also potentially results from the reduced stability. Therefore, to justify the usefulness of the proposed model, we construct a rolling prediction procedure as follows. Let $y_t$ be any of the Fama-French three factors at month $t$. We consider a forecasting model for $y_t$: \[ y_{t+1}=\alpha+\beta y_{t}+\gamma^\prime {\mbox{\boldmath $F$}}_{t+1}+e_{t+1}, \] where ${\mbox{\boldmath $F$}}_t$ is a vector of estimated factors from the Fama-French 100 portfolio data set. We estimate $\alpha,\beta,\gamma$ using ordinary least squares. For ${\mbox{\boldmath $F$}}_t$, we consider eight specifications: (i) ${\mbox{\boldmath $F$}}_t=0$, which is the benchmark AR(1) model, (ii) ${\mbox{\boldmath $F$}}_t$ from “PE”, (iii) ${\mbox{\boldmath $F$}}_t$ from “PE” and matrix quantile factor model at $\tau=0.05$, (iv) ${\mbox{\boldmath $F$}}_t$ from “PE” and matrix quantile factor model at $\tau=0.95$, (v) ${\mbox{\boldmath $F$}}_t$ from “PE” and matrix quantile factor model at $\tau=0.05$ and $\tau=0.95$, (vi) to (viii) generate ${\mbox{\boldmath $F$}}_t$ similarly to (iii) to (v) but replacing matrix quantile factors with vectorized quantile factors. To control the dimension of the design matrix, we use $k_1=k_2=1$ in this part, and ignore the case $\tau=0.5$ because the estimated loading space is very close to that from “PE” by Table (ref). To predict $y_{t+1}$, we first estimate all the factors using historical data before (inclusive) $(t+1)$ with a rolling window of 60 months, and then fit the predicting model using data only before $(t+1)$. The predictor $\hat y_{t+1}$ then follows the fitted model. Table (ref) reports the root of mean squared error (RMSE) and mean absolute error (MAE) for the prediction over all the periods, from different predicting models and for different Fama-French factors. As shown in the table, adding estimated factors into the model helps reduce the error for the SMB factor and the HML factor, while considering the edge quantiles further improves the prediction performance. This is consistent with our interpretation in Section (ref). The estimated row and column factors from the matrix quantile factor model are closely related to the Fama and French SMB and HML factors. In this example, the matrix quantile factor model leads to smaller MAE while the vectorized model leads to smaller RMSE. But for the RF factor, the Benchmark AR(1) model works already the best. Adding more factors into the predictors only results in more errors, mainly because the market excess return has already been removed from the data.
In our last experiment, we investigate the performance of matrix quantile factor model in imputing missing entries. We deliberately kick out a proportion of entries from the data, and treat them as missing values. Then, we fit a factor model, estimate the loading and factor score matrices, and impute the missing entries by the estimated common components. We calculate the imputing error in terms of RMSE, denoted by $a_1$. As a benchmark, we also calculate the imputing error when simply imputing the missing values with 0 (the data are standardized), denoted as $a_0$. For robustness check, we repeat the procedure 50 times and report the averaging $a_1/a_0$ under different factor models in Figure (ref), as the kicking-out proportion increases. It's seen that the matrix-quantile factor model with $\tau=0.5$ leads to the lowest imputing error in all scenarios.
In this study, we proposed a matrix quantile factor model that is a hybrid of quantile feature representation and a low rank structure for matrix data. By minimizing the check loss function under rotation constraints, we obtain estimates of the row and column factor spaces that are proved to be consistent in the sense of Frobenious norm. Three model selection criteria were given to consistently determine simultaneously the numbers of row and column factors. Central limit theorems are derived for the smoothed loading estimates by novelly introducing an equivalent augmented Lagrangian function. There are at least three problems that are worthy of being studied in the future. First, a statistical test for the presence of the low-rank matrix structure in the matrix quantile factor model is of potential usefulness as a model checking tool. Second, the latent factor structure here can be extended to the case where both observable explanatory variables and latent factors are incorporated into modeling the quantiles of matrix sequences. Third, the computation error with the algorithm, that parallels to the statistical error given in our theorem, is still unknown. We leave all these to our future research work.