Factors in Fashion: Factor Analysis towards the Mode

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.

91,346 characters · 16 sections · 105 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

abstractThe modal factor model represents a new factor model for dimension reduction in high dimensional panel data. Unlike the approximate factor model that targets for the mean factors, it captures factors that influence the conditional mode of the distribution of the observables. Statistical inference is developed with the aid of mode estimation, where the modal factors and the loadings are estimated through maximizing a kernel-type objective function. An easy-to-implement alternating maximization algorithm is designed to obtain the estimators numerically. Two model selection criteria are further proposed to determine the number of factors. The asymptotic properties of the proposed estimators are established under some regularity conditions. Simulations demonstrate the nice finite sample performance of our proposed estimators, even in the presence of heavy-tailed and asymmetric idiosyncratic error distributions. Finally, the application to inflation forecasting illustrates the practical merits of modal factors. {\sl JEL classification:} C38; C52; C55. {\sl Keywords}: Alternating maximization; Factor model; Mode estimation; Information criteria; Rank estimation.

\pagenumbering{Alph} \thispagestyle{empty}

\pagenumbering{arabic} \setcounter{page}{1}

Introduction

Factor model has become one of the most important tools in analyzing high dimensional time series, due to its capability of dimension reduction and feature extraction through a small number of common factors, especially in the era of big data. Theoretical advancements in factor analysis have been made using principal component analysis BN2002,Bai2003,FHLR2000,FHLR2005,AH2013,HKYZ2022 and maximum likelihood approach BL2016,Wang2022, and so on. In the meanwhile, high dimensional factor models have found practical applications in a wide range of financial and economic studies, such as modeling monetary policy BB2003, break and threshold detection MT2023a,MT2023b, group structure identification AB2017,AGP2020,ZPYZ2023, forecasting excess stock returns LN2007, bond returns LN2009 and macroeconomic time series SW2002a,SW2002b,BN2006,CH2015,TL2019,GKP2016, and among many others. See FLL2022 for a recent review and references in the above studies for more related literature.

The majority of the above theoretical contributions to factor analysis have been confined to extracting common features that explain the (conditional) mean of the observed high dimensional time series, the factors obtained from which may be referred to as mean factors. While there is no dispute that mean is one of the most commonly used location parameters for a random variable, other location measures, such mode, median, quantiles, and so on, are also frequently seen in empirical studies as they could contain alternative unique distributional information as well. To enable factor analysis across the whole (conditional) distribution, CDG2021 recently put forward quantile factors that are derived under a quantile factor model. The quantile factors are allowed to vary across the quantile level, and can completely characterize features that shift any part of the conditional distribution. The empirical evidence provided by CDG2021 demonstrates that quantile factors are very informative for density forecasting of the inflation rate and real GDP growth.

This paper advocates a modal factor model (MFM) in order to capture common features that explain the mode of the (conditional) densities of the observed high dimensional time series. This leads to modal factors that are the most likely to appear in the (conditional) densities, which are thus referred to as “fashionable factors”. It is worth noting that the modal factors are in nature quantile factors corresponding to the quantile level at which the densities of the observed time series reach their peak, just like that mode is a specific quantile at which the density achieves its maximum. In this sense, the modal factor model is nested in the quantile factor model as a special case. However, the latter does not automatically produce the former, because the quantile factors are defined only for a given quantile level $\tau$, whose value at the mode is unknown unless in certain (impractical) scenarios such as that the density is symmetric ($\tau=0.5$) or is fully known. Consequently, the quantile factor model could fail to reveal how the conditional mode of the high dimensional time series depends on the modal factors directly to detect the “most likely” effect and may produce low density point predictions, as similarly argued by UWY2023 in the regression setting. As a result, the modal factor model is potentially a very useful tool that can be of interest in itself, or used to complement the PCA and quantile factors in dimension reduction for high dimensional data.

This paper contributes to the literature in several aspects. First, a modal factor analysis (MFA) procedure, called “alternating modal expectation-maximization”, is proposed to provide a basis on which statistical inference on modal factor model can be conducted. The loss function we use to derive the modal factor and loading estimators involves a kernel function with vanishing bandwidth, which adapts that designed to obtain modal regression estimators YL2014,KS2012. The estimation procedure marries the alternating maximization algorithm used in factor estimation CDG2021 and the modal expectation-maximization algorithm adopted in estimating modal regressions YL2014. The resulting algorithm is computationally efficient and easy to implement with the choice of a normal kernel function, which largely alleviates the practical challenge that there is no analytical closed-form solution for the MFA estimators.

Second, asymptotic properties of the proposed MFA estimators are established. We derive the average convergence rate of the MFA estimators, establish their asymptotic normality, and obtain consistent estimators for the associated asymptotic variances. The results are obtained under the condition that, given the factors, the errors are independent cross-sectionally but follow an $\alpha$-mixing process time serially, without any restriction imposed on the existence of the error moments. The time serial dependence allowed largely relaxes the independence requirement made by CDG2021. We show that the MFA estimators converge at the fastest possible rate $L_{NT}^{2/7}$, where $L_{NT}=\min\{N,T\}$, with $N$ and $T$ being the cross-sectional dimension and the time length, respectively. This rate is slower than the convergence rate $L_{NT}^{1/2}$ for PCA factor estimators Bai2003 and quantile factor estimators CDG2021. The slower convergence stems from the nature of nonparametric inference due to the use of a vanishing bandwidth, and is the cost we have to pay for estimating the conditional mode without the knowledge of the density functional form Parzen1962. The optimal order of bandwidth choice is also discussed.

Third, two data-driven model selection methods, based on the rank of a certain matrix and the information criterion, respectively, are proposed to determine the number of modal factors. We characterize the conditions for the tuning parameters under which the selection for the factor number can achieve consistency. Although these selection criteria bear similarity to those used by CDG2021, there are important distinctions in the theoretical development in current modal factor models, which are outlined in the remarks. Examples of tuning parameter choices that meet the consistency requirement are also provided.

Fourth, numerical evidences are provided to demonstrate the nice finite sample performance of the proposed estimators. In particular, simulated examples show that, in spite of the reduced convergence rate, the MFA estimators can effectively capture the true factor space in a variety of parameter settings, and tend to outperform the PCA and quantile factor estimators (at $\tau=0.5$), especially when the errors are heavy-tailed. Moreover, the two factor number selection criteria can select the correct number of factors with high probability. The above simulation findings are robust to the presence of heavy-tailed or skewed errors. Finally, empirical applications illustrate that MFA factors contain valuable information in enhancing the predictive accuracy of U.S. inflation rate.

This paper is also related to the growing literature on modal regressions. KS2012 and YL2014 consider the linear modal regression through maximizing a kernel-based objective function with a vanishing bandwidth, and largely extend the pioneer work of Lee1989, Lee1993 by developing asymptotic results under skewed error distributions. For nonparametric modal regressions, YLL2012 estimate the global conditional mode by local polynomial smoothing, while CGTW2016 estimate the collection of all conditional modes based on a kernel density estimate. Recently, UWY2021 study the fixed effects modal regression for panel data, UWY2023 consider a semiparametric partially linear varying coefficient modal regression, and Wang2024 investigates the nonlinear modal regression for dependent data. It is worth emphasizing that the above studies only involve observed regressors, while both the factors and loadings are unknown and need to be estimated in the current setup. There has been no study that considers the inference on the conditional mode in factor models so far. We note that SH2020 propose a modal PCA that uses the probability density value of the mode as a measure of concentration, the direction maximizes which is regarded as the minor component (factor) direction. In that way, their formulation is notably different from ours, and they do not consider the asymptotic properties of the estimators. These differences highlight the new contribution of current paper.

The outline of the rest of this paper is as follows. Section (ref) introduces MFM, provides a list of illustrative examples of MFM, presents the MFA estimators and the computational algorithm, and proposes two methods for selecting the number of factors. Section (ref) establishes the asymptotic properties of the proposed estimators. Section (ref) evaluates the finite sample performance of the estimators using Monte Carlo simulations. Section (ref) assesses the predictive power of the MFA factors in forecasting U.S. inflation rate. The proofs of Theorems (ref), (ref) and (ref) are contained in the Appendix, while the proofs of Theorems (ref) and (ref), together with some additional simulation results, are relegated to the Supplementary Material.

Notations. For any real number $a$, $\operatorname{sgn}(a)=1$ if $a\geq 0$ and $-1$ if $a< 0$. For any matrix $\mathbf A$, let $\operatorname{rank}(\mathbf A)$, $\operatorname{tr}(\mathbf A)$, $\mathbf A^{\prime}$, $\|A\|=[\operatorname{tr}(A^{\prime}A)]^{1/2}$ and $\operatorname{vech}(\mathbf A)$ denote its rank, trace, transpose, Frobenius norm and the vectorization of $A$, respectively. For any square matrix $\mathbf B$ with real eigenvalues, denote $\rho_{\mathrm{min}}(\mathbf B)$ (resp. $\rho_{\mathrm{max}}(\mathbf B)$) as its minimum (resp. maximum) eigenvalue, $\mathbf B_{jj}$ as its $j$-th diagonal element, and $\operatorname{sgn}(\mathbf B)$ as a diagonal matrix whose $j$-th diagonal element equals $\operatorname{sgn}({\mathbf{B}}_{jj})$. We use $\mathbf B>0$ (resp. $\mathbf B<0$) to signify that $\mathbf B$ is positive (resp. negative) definite.

The model and estimation

Section (ref) presents the modal factor models. Section (ref) defines the estimators for the factors and the loadings, and introduces an algorithm to obtain the estimators. Section (ref) proposes two criteria for selecting the number of factors.

The modal factor model

Suppose that the observed variables $\{X_{it};i=1,\ldots,N,t=1,\ldots,T\}$ satisfy the following modal factor model (MFM):

equation[equation omitted — 173 chars of source]

where ${\bf{f}}_{0t}$ is an $r\times1$ vector of random common factors, $\bm{{\bm{\lambda}}}_{0i}$ is an $r\times1$ vector of non-random factor loadings, with the conditional mode function of $X_{it}$ given ${\bf{f}}_{0t}$ denoted as

equation*[equation* omitted — 132 chars of source]

and $f_{X_{it}}(\cdot|{\bf{f}}_{0t})$ being the conditional distribution of $X_{it}$ given ${\bf{f}}_{0t}$. Let $e_{it}^0=X_{it}-\bm{{\bm{\lambda}}}_{0i}^{\prime} \mathbf{f}_{0t}$ denote the idiosyncratic error, then (ref) can be equivalently represented as

equation*[equation* omitted — 288 chars of source]

where $g_{e_{it}^0}(\cdot|{\bf{f}}_{0t})$ is the conditional density of $e_{it}^0$ given ${\bf{f}}_{0t}$.

If the $g_{e_{it}^0}(\cdot|{\bf{f}}_{0t})$ is symmetric about 0 across both $i$ and $t$, then $\operatorname{Mode}(X_{it}|{\bf{f}}_{0t})=E(X_{it}|{\bf{f}}_{0t})$. In this case, the factors and loadings from the above model are the same as those of the approximate factor models (AFMs), and thus can be estimated by principal component analysis (PCA) as studied by BN2002 and Bai2003, provided that certain moment conditions hold for the error terms. However, if $g_{e_{it}^0}(\cdot|{\bf{f}}_{0t})$ is skewed for some $i$ and $t$, then these components from the two models become different. The following example provides more illustrations on the relationship among our MFM, the AFM BN2002 and the quantile factor model CDG2021.

exampleLet $\{e_{it}\}$ be independent and identically distributed (i.i.d.) errors with density function $g_e(\cdot)$ and cumulative distribution function $G_e(\cdot)$. Let $Q_e(\tau)=G_e^{-1}(\tau)=\inf\{c: G_e(c)\geq \tau)$ be the quantile function of $e_{it}$. Additionally, let ${\bf{f}}_{1t}\in \mathbb{R}^{r_1}, {\bf{f}}_{1t}\in \mathbb{R}^{r_2}$, where $r_1$ and $r_2$ are positive constants, and $\{e_{it}\}$ is independent of $\{{\bf{f}}_{1t}\}$ and $\{{\bf{f}}_{2t}\}$. Without loss of generality, suppose that $E(e_{it})=0$, $Q_{e}(\tau_0)=0$ for some $\tau_0\in [0,1]$, and $\operatorname{Mode}(e_{it})=e_m$, which corresponds to the $\tau_m$-th quantile, i.e., $Q_{e}(\tau_m)=e_m$ for some $\tau_m\in [0,1]$. Consider the following location-scale-shift factor model: $$X_{it}={\bm{\lambda}}_{i}^{\prime}{\bf{f}}_{1t}+{\bm{\alpha}}_i^{\prime}{\bf{f}}_{2t}e_{it},$$ where ${\bm{\lambda}}_i\in\mathbb{R}^{r_1}, {\bm{\alpha}}_i\in\mathbb{R}^{r_2}, {\bm{\alpha}}^{\prime}_{i}{\bf{f}}_{2t}>0$. Let ${\bf{f}}_t=[{\bf{f}}^{\prime}_{1t},{\bf{f}}^{\prime}_{2t}]^{\prime}.$ \begin{enumerate} • If ${\bf{f}}_{2t}$ and ${\bf{f}}_{1t}$ do not share any common element, then \begin{equation*} E(X_{it}|{\bf{f}}_{t})={\bm{\lambda}}_{i}^{\prime}{\bf{f}}_{1t},\qquad \operatorname{Mode}(X_{it}|{\bf{f}}_{t})={{\bm{\lambda}}_i^M}^{\prime}{\bf{f}}_{t},\qquad Q_{X_{it}}(\tau|{\bf{f}}_{t})={{\bm{\lambda}}_i(\tau)}^{\prime}{\bf{f}}_{t}, \end{equation*} where ${\bm{\lambda}}^M_i=[{\bm{\lambda}}_i^{\prime}, e_m{\bm{\alpha}}_{i}^{\prime}]^{\prime}$, ${\bm{\lambda}}_i(\tau)=[{\bm{\lambda}}_{i}^{\prime},Q_e(\tau){\bm{\alpha}}_i^{\prime}]^{\prime}$. • If ${\bf{f}}_{2t}={\bf{f}}_{1t}$, then \begin{equation*} E(X_{it}|{\bf{f}}_{t})={\bm{\lambda}}_{i}^{\prime}{\bf{f}}_{1t},\qquad \operatorname{Mode}(X_{it}|{\bf{f}}_{t})={{\bm{\lambda}}_i^M}^{\prime}{\bf{f}}_{1t},\qquad Q_{X_{it}}(\tau|{\bf{f}}_{t})={{\bm{\lambda}}_i(\tau)}^{\prime}{\bf{f}}_{1t}, \end{equation*} where ${\bm{\lambda}}^M_i={\bm{\lambda}}_i+e_m{\bm{\alpha}}_i, {\bm{\lambda}}_i(\tau)={\bm{\lambda}}_i+Q_e(\tau){\bm{\alpha}}_i$. • If $e_m=0$ and $\tau=\tau_0$, then \begin{equation*} E(X_{it}|{\bf{f}}_{t})= \operatorname{Mode}(X_{it}|{\bf{f}}_{t})= Q_{X_{it}}(\tau|{\bf{f}}_{t})={\bm{\lambda}}_{i}^{\prime}{\bf{f}}_{1t}. \end{equation*} \end{enumerate}

The above example illustrates that only in case (c) where $e_m=0$ and $\tau=\tau_0$, AFM, MFM and QFM have the same representation, while in general the three models are different. In particular, the above example demonstrates that different characteristics of $X_{it}$ can be driven by different common factors. For instance, in case (a), the conditional mean of $X_{it}$ is only affected by ${\bf{f}}_{1t}$, while the conditional mode and quantiles of $X_{it}$ are affected by both ${\bf{f}}_{1t}$ and ${\bf{f}}_{2t}$ if the location and scale factors are different. From this perspective, MFM can capture (scale) factors that are missed by AFM.

Estimating the factors and loadings

For the observed sample $\{X_{it}\}$, we take the fixed-effect approach that treats $\{\mathbf{f}_{0t}\}$ and $\{\bm{\lambda}_{0i}\}$ as unknown parameters to be estimated, while the asymptotic analysis is conditional on $\{\mathbf{f}_{0t}\}$. We first assume that the factor number $r_0$ is known in this subsection. The data-driven selection of $r_0$ will be presented in Section (ref).

Let $M=(N+T) r_0, \bm{{\bm{\theta}}}=\left({\bm{{\bm{\lambda}}}}_{1}^{\prime}, \ldots, {\bm{{\bm{\lambda}}}}_{N}^{\prime}, {\bf{f}}_{1}^{\prime}, \ldots, {\bf{f}}_{T}^{\prime}\right)^{\prime}$. Denote $\bm{{\bm{\theta}}}_{0}=\left(\bm{{\bm{\lambda}}}_{01}^{\prime}, \ldots, \bm{{\bm{\lambda}}}_{0 N}^{\prime}, \mathbf{f}_{01}^{\prime}, \ldots, \mathbf{f}_{0 T}^{\prime}\right)^{\prime} $ as the vector of true parameters, and write ${\bm{{\bm{\Lambda}}}}_0=(\bm{{\bm{\lambda}}}_{01},\cdots,{\bm{{\bm{\lambda}}}}_{0N})^{\prime}$, ${\bf{F}}_0=({\bf{f}}_{01},\cdots,{\bf{f}}_{0T})^{\prime}.$ The dependence of $\bm{{\bm{\theta}}} $ and $\bm{{\bm{\theta}}}_0 $ on $M$ is suppressed for notational convenience.

A well-known fact in factor models BN2002 is that $\left\{\bm{{\bm{\lambda}}}_{0i}\right\} $ and $ \left\{{\bf{f}}_{0t}\right\} $ cannot be separately identified without imposing certain normalizations. Without loss of generality, we adopt the following normalizations:

align[align omitted — 284 chars of source]

Let $\mathcal{A}, \mathcal{F} \subset \mathbb{R}^{r_0}$ and define the parameter space as

equation*[equation* omitted — 317 chars of source]

To estimate the modal factor model parameter $\bm{{\bm{\theta}}}_0$, we follow the strategy of KS2012 and YL2014 to use a kernel-based objective function in the regression setting. Specifically, define

equation[equation omitted — 319 chars of source]

where $K_h(u)=\frac{1}{h}K(\frac{u}{h})$, $K(\cdot)$ is a smooth kernel function, and $h$ is a bandwidth diminishing towards $0$ as $N, T\to \infty$. The estimator of $\bm{{\bm{\theta}}}_0$ is then defined as

equation[equation omitted — 338 chars of source]

The reason why we estimate ${\bm{{\bm{\theta}}}}_0$ in this way is related to the fact that, for any fixed $\bm{{\bm{\theta}}}$, $\mathbb{M}_{NT}(\bm{{\bm{\theta}}})$ can be seen as a kernel density estimator for the residuals $e_{it}=X_{it}-{\bm{{\bm{\lambda}}}}_i^{\prime}{\bf{f}}_t$ at $0$ (i.e. $g_{e_{it}}(0)$). Note that $g_{e_{it}}(0)=E\big[g_{e_{it}}(0|{\bf{f}}_{0t})\big]=E\big[f_{X_{it}}({\bm{{\bm{\lambda}}}}_i^{\prime}{\bf{f}}_t|{\bf{f}}_{0t})\big]\leq E\big[f_{X_{it}}({\bm{{\bm{\lambda}}}}_{0i}^{\prime}{\bf{f}}_{0t}|{\bf{f}}_{0t})\big]=E\big[g_{e_{it}^0}(0|{\bf{f}}_{0t})\big]$, provided that $f_{X_{it}}({\bm{{\bm{\lambda}}}}_i^{\prime}{\bf{f}}_t|{\bf{f}}_{0t})\leq f_{X_{it}}({\bm{{\bm{\lambda}}}}_{0i}^{\prime}{\bf{f}}_{0t}|{\bf{f}}_{0t})$ for all $\bm{\theta}\in{\bm{\Theta}}^{r_0}$, with a strict inequality when ${\bm{{\bm{\lambda}}}}_i^{\prime}{\bf{f}}_t\neq {\bm{{\bm{\lambda}}}}_{0i}^{\prime}{\bf{f}}_{0t}$. This suggests that $\bm{{\bm{\theta}}}_0$ can be well estimated by the vector at which the kernel density function $\mathbb{M}_{NT}(\bm{{\bm{\theta}}})$ reaches the peak. See KS2012 and YL2014 for illustrations of such criterion function choices in the regression setting.

As typical in modal regressions, the maximization in (ref) does not yield an analytical closed form for $\hat{\bm{{\bm{\theta}}}}$. For practical implementation, we propose an algorithm, which we call alternating modal expectation-maximization (AMEM), to obtain the numerical estimate for an observed sample $\{X_{it}\}$. As both the factor and loading vectors are unknown and need to be estimated, the algorithm estimates one by maximizing the objective function in (ref) given the other and then alternates, until some stopping rule is satisfied. Each maximization involved in the iterated process is essentially a linear modal regression estimation, which can be conveniently solved by the modal expectation-maximization algorithm proposed by YL2014.

To precisely describe the algorithm, let $\bm{{\bm{\Lambda}}}=({\bm{{\bm{\lambda}}}}_1,\cdots, {\bm{{\bm{\lambda}}}}_N)^{\prime}, {\bf{F}}=({\bf{f}}_1,\cdots, {\bf{f}}_T)^{\prime}$, and write

equation*[equation* omitted — 304 chars of source]

Then $\mathbb{M}_{NT}(\bm{{\bm{\theta}}})=N^{-1}\sum_{i=1}^N\mathbb{M}_{i,T}(\bm{{\bm{\lambda}}}_i,\mathbf{F})=T^{-1}\sum_{t=1}^T\mathbb{M}_{t,N}(\bm{{\bm{\Lambda}}},{\bf{f}}_t)$.

\noindentThe AMEM algorithm:

enumerate[{\bf Step 1}] • Choose random starting parameters: $\mathbf{F}^{(0)}, \bm{\Lambda}^{(0)}$. • For $l\geq1$, with given $\mathbf{F}^{(l-1)}$, solve $\bm{{\bm{\lambda}}}_i^{(l)}=\arg\underset{\bm{{\bm{\lambda}}}}{\max}\,\,\mathbb{M}_{i,T}(\bm{{\bm{\lambda}}},\mathbf{F}^{(l-1)})$ for $i=1,\cdots,N$. With initial value $\bm{{\bm{\lambda}}}_i(0)={\bm{\lambda}}_i^{(l-1)}$, the solution $\bm{{\bm{\lambda}}}_i^{(l)}$ can be found by repeating the following two steps until convergence: E-Step: Calculate weights $\pi(t|\bm{{\bm{\lambda}}}_i(k),\mathbf{F}^{(l-1)}),$ for $t=1,\cdots,T$ as \begin{equation} \pi(t|\bm{{\bm{\lambda}}}_i (k),\mathbf{F}^{(l-1)})=\frac{K_h(X_{it}-\bm{{\bm{\lambda}}}_i (k)^{\prime} \mathbf{f}_t^{(l-1)})}{\sum_{t=1}^{T}K_h(X_{it}-\bm{{\bm{\lambda}}}_i (k)^{\prime} \mathbf{f}_t^{(l-1)})}\propto K_h(X_{it}-\bm{{\bm{\lambda}}}_i (k)^{\prime} \mathbf{f}_t^{(l-1)}). \end{equation} M-Step: Update $\bm{{\bm{\lambda}}}_i (k+1)$ as \begin{align} \bm{{\bm{\lambda}}}_i (k+1)&=\arg\underset{\bm{{\bm{\lambda}}}}{\max}\,\,\sum\limits_{t=1}^{T}\big\{\pi(t|\bm{{\bm{\lambda}}}_i (k),\mathbf{F}^{(l-1)})\log K_h(X_{it}-\bm{{\bm{\lambda}}}^{\prime} \mathbf{f}_t^{(l-1)})\big\}\nonumber\\ &=(\mathbf{F}^{(l-1)^{T}}\mathbf{W}_k\mathbf{F}^{(l-1)})^{-1}\mathbf{F}^{(l-1)^{T}}\mathbf{W}_k\mathbf{X}_i, \end{align} where $\mathbf{X}_i=(X_{i1},\cdots,X_{iT})^{\prime}$ and $\mathbf{W}_k$ is a $T\times T$ diagonal matrix with $t$-th diagonal element $\pi(t|\bm{{\bm{\lambda}}}_i (k),\mathbf{F}^{(l-1)})$. • Given $\bm{{\bm{\Lambda}}}^{(l)}$, solve $\mathbf{f}_t^{(l)}=\arg\underset{\mathbf{f}}{\max}\,\,\mathbb{M}_{t,N}(\bm{{\bm{\Lambda}}}^{(l)}, \mathbf{f})$ for $t=1,\cdots,T$ following a similar procedure as outlined in Step 2, with initial value $\bm{{\bf{f}}}_t(0)=\bm{{\bf{f}}}_t^{(l-1)}$. • For $l=1,2,\cdots$, iterate Steps 2-3 until $\mathbb{M}_{NT}(\bm{{\bm{\theta}}}^{(L)})$ is close to $\mathbb{M}_{NT}(\bm{{\bm{\theta}}}^{(L-1)})$ for some $L$, i.e., $|\mathbb{M}_{NT}(\bm{{\bm{\theta}}}^{(L)})-\mathbb{M}_{NT}(\bm{{\bm{\theta}}}^{(L-1)})|<\epsilon$, for some small positive $\epsilon$, where $\bm{{\bm{\theta}}}^{(l)}=\operatorname{vech}((\bm{{\bm{\Lambda}}}^{(l)})^{\prime},(\mathbf{F}^{(l)})^{\prime})$. • Normalize $\bm{{\bm{\Lambda}}}^{(L)}$ and $\mathbf{F}^{(L)}$ to satisfy the normalizations in (ref).

Note that the closed-form solution as shown in (ref) above during the M-Step is only obtained for the standard normal kernel choice $K(u)=\phi(u)$, which largely enhances the computational efficiency. The asymptotic results obtained in this article remain valid for other kernel choices as well, even though in general they do not produce explicit solution as in (ref), in which cases numerical optimization becomes necessary.

remark\begin{enumerate}[(i)] • The AMEM algorithm is a hybrid procedure that combines the alternating maximization (AM) algorithm and the modal expectation-maximization (MEM) algorithm. The former refers to the process of alternately maximizing $\mathbb{M}_{NT}(\bm{{\bm{\theta}}})$ with respect to $\bm{{\bm{\Lambda}}}$ or $\bf{F}$ given the other, while the latter algorithm is borrowed from YL2014 and solves each maximization as a linear modal regression. The idea of AM algorithm has been utilized in a variety of factor models in which no explicit solutions are available for the factor and loading estimators, including the quantile factor model of CDG2021 and the generalized factor model of Wang2022, among others. • Following the arguments for Theorem 2.1 of YL2014, we can establish the monotonically ascending property of the objective function in Steps 2 and 3, i.e., $\mathbb{M}_{i,T}(\bm{{\bm{\lambda}}}_i(k+1),{\bf{F}}^{(l-1)})\geq \mathbb{M}_{i,T}({\bm{{\bm{\lambda}}}_i}(k),{\bf{F}}^{(l-1)})$, $ \mathbb{M}_{t,N}({\bm{{\bm{\Lambda}}}}^{(l)},{\bf{f}}_t(k+1))\geq \mathbb{M}_{t,N}({\bm{{\bm{\Lambda}}}}^{(l)},{\bf{f}}_t(k))$ for $i=1,\cdots,N; t=1,\cdots,T$, and for each $l$ and $k$. Consequently, we have that $\mathbb{M}_{NT}(\bm{{\bm{\theta}}}|{\bm{{\bm{\Lambda}}}}^{(l+1)},{\bf{F}}^{(l+1)})\geq \mathbb{M}_{NT}(\bm{{\bm{\theta}}}|{\bm{{\bm{\Lambda}}}}^{(l+1)}, {\bf{F}}^{(l)})$ $\geq \mathbb{M}_{NT}(\bm{{\bm{\theta}}}|{\bm{{\bm{\Lambda}}}}^{(l)},{\bf{F}}^{(l)})$, for $l=0,1,\cdots, L$, which guarantees the convergence of the AMEM algorithm. • For both the AM algorithm and the MEM algorithm, the solutions obtained upon convergence are necessarily local maxima. As a result, the converged value obtained by the AMEM algorithm depends on the starting points. To find the global optimum, it is necessary to run the algorithm multiple times with different initial values of ${\bf{F}}^{(0)}$, $\bm{{\bm{\Lambda}}}^{(0)}$, and then choose the best local maximum. \end{enumerate}

Selecting the number of factors

The previous subsection assumes that the number of factors $r_0$ is known, which is needed to outline the estimation algorithm for the factors and loadings. In practice, data-driven procedures are desired to select the correct $r_0$. In the following, we propose two methods that can consistently determine $r_0$ based on the observed sample.

For the sake of exposition, we first introduce some notations. Let $\bar{r}$ be a large positive integer such that $r_0<\bar{r}<\infty$. For any $r=1,\ldots, \bar{r}$, let $\mathcal{A}^r$ and $\mathcal{F}^r$ be compact subsets of $\mathbb{R}^r$. Let ${\bm{{\bm{\lambda}}}}_i^r, {\bf{f}}_t^r\in\mathbb{R}^r$ for $i=1,\cdots, N, t=1,\cdots,T$, and write ${\bm{{\bm{\theta}}}}^r=({{\bm{{\bm{\lambda}}}}_1^r}^{\prime},\cdots,{{\bm{{\bm{\lambda}}}}_N^r}^{\prime}, {{\bf{f}}_1^r}^{\prime}, \cdots, {{\bf{f}}_T^r}^{\prime})^{\prime}, {\bm{{\bm{\Lambda}}}}^r=({{\bm{{\bm{\lambda}}}}_1^r},\cdots,{{\bm{{\bm{\lambda}}}}_N^r})^{\prime}, {\bf{F}}^r=({{\bf{f}}_1^r}, \cdots, {{\bf{f}}_T^r})^{\prime}$. Similar to (ref), we adopt the following normalizations:

align[align omitted — 298 chars of source]

Define ${\bm{\Theta}}^{r}=\left\{\bm{{\bm{\theta}}}^r : \bm{{\bm{\lambda}}}^r_i\in \mathcal{A}^r, \mathbf{f}^r_{t} \in \mathcal{F}^r \text{ for all } i,t,\left\{\bm{{\bm{\lambda}}}^r_{i}\right\} \text{and} \left\{\mathbf{f}^r_{t}\right\} \text{satisfy } (\ref{81}) \right\}$, and denote

equation[equation omitted — 392 chars of source]

Further, let $\hat{\mathbf{{\bm{\Lambda}}}}^{\bar{r}}=(\hat{\bm{{\bm{\lambda}}}}^{\bar{r}}_1,\cdots,\hat{\bm{{\bm{\lambda}}}}^{\bar{r}}_N)^{\prime}$ and write

equation*[equation* omitted — 230 chars of source]

Our first estimator of $r_0$ is related to the rank of the above matrix. In particular, the rank estimator $\hat{r}_{\mathrm{rank}}$ is defined as

equation*[equation* omitted — 132 chars of source]

where $P_{1,N T}$ is a sequence that approaches 0 as $N, T \rightarrow \infty$. It is easily seen that $\hat{r}_{\mathrm{rank }}$ equals the number of diagonal elements in $\left(\hat{\bm{{\bm{\Lambda}}}}^{\bar{r}}\right)^{\prime} \hat{\bm{{\bm{\Lambda}}}}^{\bar{r}}/N$ that exceed the threshold $P_{1,N T}$. Consequently, $\hat{r}_{\mathrm{rank}}$ can be interpreted as the rank estimator of $\left(\hat{\bm{{\bm{\Lambda}}}}^{\bar{r}}\right)^{\prime} \hat{\bm{{\bm{\Lambda}}}}^{\bar{r}} / N$. It will be shown that $\hat{r}_{\text {rank }}$ approaches $r_0$ as $N, T\to \infty$, provided that the tuning parameter $P_{1,N T}$ diminishes at certain rates.

Our second estimator of $r_0$ is derived from the information criterion (IC), inspired from BN2002. To be specific, we consider the following IC:

equation[equation omitted — 148 chars of source]

where $P_{2,NT}$ is a sequence that approaches 0 as $N, T \rightarrow \infty$. The IC-based estimator $\hat{r}_{\mathrm{IC}}$ of $r_0$ is defined as

equation*[equation* omitted — 108 chars of source]

It is worth noting that the computation for $\hat{r}_{\mathrm{rank}}$ is significantly less demanding than that for $\hat{r}_{\mathrm{IC}}$. This is because for $\hat{r}_{\mathrm{rank}}$ the MFM only needs to be estimated once with the number of factors setting as $\bar{r}$, while for $\hat{r}_{\mathrm{IC}}$ the MFM needs to be estimated $\bar{r}$ times, with the number of factors specified as $r=1,\ldots,\bar{r}$, respectively. This indicates that $\hat{r}_{\mathrm{rank}}$ is practical preferred in terms of computational time, particularly with large sample sizes.

Asymptotic properties

This section establishes the consistency and asymptotic normality of the proposed estimators for the factors and loadings, and the selection consistency of the two criteria for selecting the factor number.

Consistency

The following assumptions are needed to facilitate the theoretical development.

assumption(Factors and factor loadings) Suppose that the parameter spaces $\mathcal{A}$ and $\mathcal{F}$ are compact, and ${\bm{{\bm{\theta}}}}_{0}$ is an interior point of ${\bm{\Theta}}^{r_0}$. Further, it holds that \begin{enumerate}[(i)] • $T^{-1}\sum_{t=1}^T{\bf{f}}_{0t}^{\prime}{\bf{f}}_{0t}\overset{p}{\to}\mathbb{I}_{r_0}$, and there exists a finite positive constant $M_1$, such that $\|\mathbf{f}_{0t}\|\leq M_1$ for all $t=1,\ldots,T$. • $N^{-1}\sum_{i=1}^{N}\bm{{\bm{\lambda}}}_{0i}^{\prime}{\bm{\lambda}}_{0i}=\operatorname{diag}(\sigma_{N1},\cdots,\sigma_{Nr_0})$, $\sigma_{N1}\geq\sigma_{N2}\cdots\geq \sigma_{Nr_0},$ and $\sigma_{Nj} \to \sigma_{j}$ as $N \to \infty$ for $j=1,\cdots, r_0$, with $\infty>\sigma_{1}>\sigma_{2}\cdots >\sigma_{r_0}>0$. In addition, there exists a finite positive constant $M_2$ such that $\|\bm{{\bm{\lambda}}}_{0i} \|\leq M_2 $ for all $i=1,\ldots,N$. \end{enumerate}

assumption(Cross-section dependence and heteroskedasticity) \begin{enumerate}[(i)] • Given $\{\mathbf{f}_{0t}, 1\leq t\leq T\}, \{e^0_{it}, i=1,\cdots,N; t=1, \cdots, T\} $ are independent across $i$. • The sequence $\{{\bf{f}}_{0t},{\bf{e}}^0_t\}$ is $\alpha$-mixing, with mixing coefficients $\alpha(k)\leq B\rho^k$, for $\rho\in (0,1)$ and $B>0$, where ${\bf{e}}^0_t=(e^0_{1t},\cdots,e_{Nt}^0)^{\prime}.$ \end{enumerate}

assumption(Conditional density) For notational convenience, let $g_{it}(\cdot)=g_{e^0_{it}}(\cdot|{\bf{f}}_{0t})$, \begin{enumerate}[(i)] • $g_{it}(\cdot)$ is continuous and uniformly bounded for all $i,t$. • For any positive real number $C$, there exists $\underline{g}>0$ (depending on $C$) such that $g_{it}(0)-\sup_{|u|\geq C}g_{it}(u)\geq \underline{g}$ for all $i,t$. • $g_{it}(\cdot)$ is three times continuously differentiable, with $g_{it}^{(v)}(u)={\partial}^vg_{it}(u)/\partial u^v$ being uniformly bounded for $v=1,2,3$ and for all $i,t$. • $g_{it}^{(1)}(0)=0$, $-g_{it}^{(2)}(0)\geq \underline{g_1}$ for some $\underline{g_1}>0$, and for all $i,t$. \end{enumerate}

assumption(Joint density) Let $g_{it,js}(\cdot, \cdot)$ denote the joint conditional density function of $e^0_{it}$ and $e^0_{js}, (i,t)\neq (j,s)$, given ${\bf{f}}_{0t}$ and ${\bf{f}}_{0s}$. $g_{it,js}(\cdot, \cdot)$ is uniformly bounded for all $i,t,j,s, (i,t)\neq (j,s)$.

assumption(Kernel function) $K(\cdot): \mathbb{R}\to \mathbb{R}$ is a twice continuously differentiable kernel function such that (i) $\int_{-\infty}^{\infty}K(u)du=1$; (ii) $K(\cdot)$ is symmetric about $0$. (iii) $\lim_{u\to \pm \infty} K(u)=0$; (iv) $\sup_{u\in \mathbb{R}}|K(u)|=c_0<\infty$; (v) $\sup_{u\in \mathbb{R}}|K^{(1)}(u)|=c_1<\infty$, where $K^{(1)}(u)=dK(u)/du$; (vi) $\sup_{u\in \mathbb{R}}|K^{(2)}(u)|=c_2<\infty$, where $K^{(2)}(u)=d^2K(u)/du^2$; (vii) $\int_{-\infty}^{\infty}|K(u)|u^2du=L_0<\infty$; (viii) $\int_{-\infty}^{\infty}|K^{(1)}(u)|du=L_1<\infty$.

assumption(Bandwidth) Let $L_{NT}=\min\{N,T\},U_{NT}=\max\{N,T\}.$ As $N, T \to \infty$, we have $ \log L_{NT}/(L_{NT}h^3)\to 0$ and $U_{NT}h^{\gamma}\to 0$ for some $\gamma>3$.\footnote{Assumption (ref) requires $U_{NT}\leq L_{NT}^{\gamma/3}$ for some $\gamma>3$, which is quite weak when $\gamma$ is large.}

Assumption (ref) places conditions on the factors and loadings following Bai2003. Assumption (ref) (i) requires the factors to be non-degenerate, while Assumption (ref) (ii) ensures that each factor has a nontrivial contribution and can be ordered according to their contributions. The requirement that true parameter is an interior point of the compact parameter space has been similarly made in BL2016 and Wang2022, and is common in nonlinear models where the estimators have no explicit forms.

Assumption (ref) imposes restrictions on the dependence of the idiosyncratic errors. In particular, Assumption (ref) (i) requires the errors to be cross-sectionally independent, and Assumption (ref) (ii) assumes that the errors follow an $\alpha$-mixing process over time. The mixing property resembles that allowed in UWY2021 and largely relaxes the independence requirement imposed by KS2012 and YL2014 for modal regressions and by CDG2021 for quantile factor models. It is essential for establishing the bound for the sum of dependent variables. Additionally, both cross-sectional and time-series heteroscedasticity are allowed.

Assumption (ref) imposes smoothness conditions on the conditional density of the error term. Assumption (ref) (ii) requires that the idiosyncratic error has a well defined unique global mode at $0$ KS2012,UWY2021, which is necessary for parameter identification. Note that $g_{it}(\cdot)$ is allowed to be heterogeneous over $i$ and $t$, and it is not required to be unimodal. Assumption (ref) (iii) assumes the derivatives of the conditional density function up to three order are uniformly bounded KS2012, which is needed to control the remainder term in the Taylor expansion. Assumption (ref) (iv) supposes the conditional density function is concave over a neighbourhood of the mode KPS2020, which is required to determine the sign of a specific term in the Taylor expansion. It is worth emphasizing that the existence of moments of the errors is not needed, in contrast to Bai2003. Assumption (ref) stipulates that the joint conditional density of the error terms given the factors should be bounded, which is needed to control the sum of covariances when calculating the variance for the sum of dependent variables Masry1996.

The kernel function under Assumption (ref) is a bounded density function that is symmetric about the mode YL2014, with tail approaching $0$ and has bounded first two derivatives KS2012,KPS2020. Assumptions (ref) (vii) and (viii) are needed when calculating the moments of some random quantities KPS2020. Assumption (ref) specifies the convergence rate for the bandwidth required for asymptotic analysis, which is similar to those made by Romano1988 and KS2012.

Write $\hat{\bm{{\bm{\Lambda}}}}=(\hat{\bm{{\bm{\lambda}}}}_1,\cdots,\hat{\bm{{\bm{\lambda}}}}_N)^{\prime}$, $\hat{\bf{F}}=(\hat{\bf{f}}_{1},\cdots,\hat{\bf{f}}_{T})^{\prime}.$ For any $\bm{\theta}_{a},\bm{\theta}_{b} \in \bm{\Theta}^{r_0},$ let $\bm{\theta}_{a}=({\bm{\lambda}}^{\prime}_{a1}, \cdots, {\bm{\lambda}}^{\prime}_{aN},$ $ {\bf{f}}^{\prime}_{a1},\cdots,{\bf{f}}^{\prime}_{aT})^{\prime}$, $\bm{\theta}_{b}=({\bm{\lambda}}^{\prime}_{b1}, \cdots, {\bm{\lambda}}^{\prime}_{bN}, {\bf{f}}^{\prime}_{b1},\cdots,{\bf{f}}^{\prime}_{bT})^{\prime}$, and define $$d({\bm{\theta}}_a,\bm{\theta}_{b})=\frac{1}{NT}\sum_{i=1}^N\sum_{t=1}^T({\bm{\lambda}}^{\prime}_{ai}{\bf{f}}_{at}-{\bm{\lambda}}^{\prime}_{bi}{\bf{f}}_{bt})^2,$$ which measures the distance between the common components of ${\bm{\theta}}_a$ and ${\bm{\theta}}_b$. The following theorem establishes the average convergence rate of $\hat{\bf{F}}$ and $\hat{\bm{\Lambda}}$.

theoremUnder Assumptions (ref)-(ref), as $N,T \to \infty$, it holds that \begin{enumerate}[(a)] • $\|\hat{{\bm{\Lambda}}}-{\bm{\Lambda}}_0 \hat{\bf{S}}\|/\sqrt{N}=O_p\left(\frac{1}{\sqrt{L_{NT}h^3}}+h^2\right)$; • $\|\hat{\mathbf{F}}-\mathbf{F}_0 \hat{\bf{S}}\|/\sqrt{T}=O_p\left(\frac{1}{\sqrt{L_{NT}h^3}}+h^2\right)$; • $d(\hat{\bm{\theta}},{\bm{\theta}}_0)=O_p\left(\frac{1}{\sqrt{L_{NT}h^3}}+h^2\right)$, \end{enumerate} with $\hat{\bf{S}}=\operatorname{sgn}(\hat{\bf{F}}^{\prime}{\bf{F}}_0)$, which appears because the value of $\hat{\bm{{\bm{\lambda}}}}_i^{\prime}\hat{\bf{f}}_t$ remains unchanged if both $\hat{\bm{{\bm{\lambda}}}}_i$ and $\hat{\bf{f}}_t$ are multiplied by $-1$.

remark\begin{enumerate}[(i)] • There are two notable differences between our theoretical study for MFM and those in the modal regression setting. First, unlike in the modal regressions KS2012,YL2014 where regressors are observed, both factors and loadings are unobserved in our setup. This calls for simultaneous inference on both the factor and loading estimators. Second, the number of parameters of interest in the regression setting is often of finite dimension, while the number of parameters in MFM diverges along with both $N$ and $T$. Such distinctions prevent the use of proof strategies designed for modal regression in our theoretical development. • The proof of Theorem (ref) borrows asymptotic techniques developed by CDG2021 for analyzing QFM, given that both estimation procedures involve iterative estimation of diverging number of parameters. However, compared to CDG2021, there exist at least two major innovative aspects in the theoretical development. First, the crucial inequality to bound $d^2(\hat{\bm{\theta}},{\bm{\theta}}_0)$ in CDG2021 does not hold in our setup. Instead, Taylor expansion and properties of the error density functions are used to establish the bound. Second, we remove the restrictive time serial independence error assumption imposed by CDG2021 and replace it with a mixing condition in the asymptotic development, under which the exponential-type inequality Bosq2012 remains effective to control the tail probabilities. • Theorem (ref) reveals that the optimal bandwidth order for estimation is $h_{opt}=O_p(L_{NT}^{-1/7})$, with which the fastest average convergence rates of $\hat{\bm{\Lambda}}$ and $\hat{\bf{F}}$ are both $L_{NT}^{2/7}$. This rate is slower than the typical rate $L_{NT}^{1/2}$ obtained for factor and loading estimators for the AFM BN2002, the QFM CDG2021 or generalized factor model Wang2022. This finding aligns with the result in the regression framework, where the fastest convergence rates for the modal regression estimator KS2012,YL2014 is $n^{2/7}$ and that for the mean regression estimator is $n^{1/2}$, with $n$ representing the sample size. As explained by UWY2022, such reduced convergence rate is attributed to the use of a shrinking bandwidth, which makes the modal estimators rely only on the observations in a small neighbourhood of the mode. In spite of this, the simulations in Section (ref) indicate that the proposed estimators enjoy desirable estimation accuracy compared to other alternatives. \end{enumerate}

Asymptotic normality

We next study the distributional properties of the factor and loading estimators. The following additional assumptions are needed.

assumption(Kernel function) $K(\cdot)$ is three times continuously differentiable, such that (i) $\sup_{u\in \mathbb{R}}|K^{(3)}(u)|=c_3<\infty$, where $K^{(3)}(u)=d^3K(u)/du^3$; (ii) $\int_{-\infty}^{\infty}|K(u)| |u|^5du=L_2<\infty$; (iii) $\int_{-\infty}^{\infty}|K^{(1)}(u)|^2|u| du=L_3<\infty$; (iv)$\int_{-\infty}^{\infty}|K^{(2)}(u)|^2du=L_4<\infty$.

assumption(Conditional density) \begin{enumerate}[(i)] • $g_{it}(\cdot)$ is six times continuously differentiable, with $g_{it}^{(v)}(u)={\partial}^vg_{it}(u)/\partial u^v$ being uniformly bounded for $v=4,5,6$ and for all $i,t$. • $g_{it}^{(v)}(0)=0$ for $v=3,5$ and for all $i,t$. \end{enumerate}

assumption(Bandwidth) As $N,T\to \infty, N \propto T$, $Th^{11}\to\infty$, $Th^{13}\to 0$.

Assumptions (ref)-(ref) strengthen the conditions used in establishing the consistency results earlier. Assumptions (ref) and (ref) require that the kernel function and the conditional density of the error terms should exhibit higher order smoothness. Assumption (ref) is standard in modal regression for asymptotic normality KS2012,KPS2020. Assumption (ref) is stronger than Assumption B3 in KS2012, and is required to control the higher order terms in the stochastic expansions of the estimators. Assumption (ref) specifies a suitable rate for the bandwidth, which ensures the modal estimators to be asymptotically unbiased.

Define

equation*[equation* omitted — 297 chars of source]

and

equation*[equation* omitted — 291 chars of source]

where $L=\int_{-\infty}^{\infty}|K^{(1)}(u)|^2du<\infty$.

assumption(Matrix definiteness) \begin{enumerate}[(i)] • ${\bm{{\bm{\Phi}}}}_i<0$ and ${\bm{{\bm{\Psi}}}}_t<0$ for all $i,t$. • ${\bm{{\bm{\Sigma}}}}_i>0$ and ${\bm{{\bm{\Omega}}}}_t>0$ for all $i,t$. \end{enumerate}

Assumption (ref) requires the matrices to be positive or negative definite, which is essential for establishing the asymptotic variances of the estimators CDG2021.

theoremLet $\hat{\bf{S}}=\operatorname{sgn}(\hat{\bf{F}}^{\prime}{\bf{F}}_0/T)$. Under Assumptions (ref)-(ref), as $N,T\to \infty$, \begin{equation*} \sqrt{Th^3}(\hat{\bm{\lambda}}_i-\hat{\bf{S}}{\bm{\lambda}}_{0i})\overset{d}{\rightarrow}\mathcal{N}(0,{\bm{\Phi}}_i^{-1}{\bm{\Sigma}}_i{\bm{\Phi}}_i^{-1}), \qquad \sqrt{Nh^3}(\hat{\bf{f}}_t-\hat{\bf{S}}{\bf{f}}_{0t})\overset{d}{\rightarrow}\mathcal{N}(0,{\bm{\Psi}}_t^{-1}{\bm{\Omega}}_t{\bm{\Psi}}_t^{-1}). \end{equation*}

remark\begin{enumerate}[(i)] • Theorem (ref) establishes the limiting distributions of $\hat{\bm{\lambda}}_i$ and $\hat{\bf{f}}_t$. The asymptotic variance for $\hat{\bm{\lambda}}_i$ and $\hat{\bf{f}}_t$ are similar to those for regression coefficient estimator in KS2012 with observed regressors. • The proof of this theorem relies on expanding the first order conditions around the true values of the factors and loadings. Taking the first result for example. Let $K_h^{(j)}(u)=d^jK_h(u)/du^j,j=1,2$. The first order condition $\sum_{t=1}^TK_h^{(1)}(X_{it}-\hat{\bm{\lambda}}_i^{\prime}\hat{\bf{f}}_t)\hat{\bf{f}}_t/T=0$ expanded around $(\hat{\bf{S}}{\bm{\lambda}}_{0i}, {\bf{F}}_0\hat{\bf{S}})$ gives rise to, after some simple calculations, \begin{align} &\frac{1}{T}\sum_{t=1}^TK^{(2)}_h(e_{it}^0){\bf{f}}_{0t}{\bf{f}}_{0t}^{\prime}(\hat{\bm{\lambda}}_i-\hat{\bf{S}}{\bm{\lambda}}_{0i})\nonumber\\ =&\frac{1}{T}\sum_{t=1}^TK^{(1)}_h(e_{it}^0)\hat{\bf{S}}{\bf{f}}_{0t}+\frac{1}{T}\sum_{t=1}^TK_{h}^{(1)}(e_{it}^{1})(\hat{\bf{f}}_t-\hat{\bf{S}}{\bf{f}}_{0t})-\frac{1}{T}\sum_{t=1}^TK_{h}^{(2)}(e_{it}^{2}){\bf{f}}_{0t}{\bm{\lambda}}_{0i}^{\prime}(\hat{\bf{f}}_t-\hat{\bf{S}}{\bf{f}}_{0t})\nonumber\\ +&o_p(\|\hat{\bm{\lambda}}_{i}-\hat{\bf{S}}{\bm{\lambda}}_{0i}\|), \end{align} where $e_{it}^0=X_{it}-{\bm{\lambda}}_{0i}^{\prime}{\bf{f}}_{0t},\hat{e}_{it}=X_{it}-\hat{\bm{\lambda}}_i^{\prime}\hat{\bf{f}}_t$, and $e_{it}^1, e_{it}^2$ lie between $e_{it}^0$ and $\hat{e}_{it}$. With the above expansion, the proof then proceeds in two steps. The first step involves demonstrating that both the second and the third terms on the right hand side of (ref) are $o_p(1/\sqrt{Th^3})$. Since $\hat{\bf{f}}_t-\hat{\bf{S}}{\bf{f}}_{0t}$ does not yield an analytical form, we require a stochastic expansion for $ \hat{\bf{f}}_t-\hat{\bf{S}}{\bf{f}}_{0t}$. This technical challenge is solved by showing that the expected Hessian matrix is asymptotically block diagonal. In this way, (ref) can be simplified as \begin{equation*} \frac{1}{T}\sum_{t=1}^TK^{(2)}_h(e_{it}^0){\bf{f}}_{0t}{\bf{f}}_{0t}^{\prime}(\hat{\bm{\lambda}}_i-\hat{\bf{S}}{\bm{\lambda}}_{0i})=\frac{1}{T}\sum_{t=1}^T{K}_h^{(1)}(e^0_{it})\hat{\bf{S}}{\bf{f}}_{0t}+o_p(1/\sqrt{Th^3})+o_p(\|\hat{\bm{\lambda}}_i-\hat{\bf{S}}{\bm{\lambda}}_{0i}\|). \end{equation*} The second step is then to establish the asymptotic normality of $\sum_{t=1}^T{K}_h^{(1)}(e^0_{it})\hat{\bf{S}}{\bf{f}}_{0t}/T$ under the mixing condition, following the proof strategy used by Masry1996. Consequently, the asymptotic normality of the loading estimators can be established. \end{enumerate}

Define

equation*[equation* omitted — 346 chars of source]

and

equation*[equation* omitted — 376 chars of source]

The following theorem establishes the consistency of the variance estimators for those of the factor and loading estimators.

theoremUnder Assumptions (ref)-(ref), as $N,T\to \infty$, it holds that \begin{equation*} \hat{\bm{\Phi}}_i^{-1}\hat{\bm{\Sigma}}_i\hat{\bm{\Phi}}_i^{-1}\overset{p}{\rightarrow} {\bm{\Phi}}_i^{-1}{\bm{\Sigma}}_i{\bm{\Phi}}_i^{-1}, \qquad \hat{\bm{\Psi}}_t^{-1}\hat{\bm{\Omega}}_t\hat{\bm{\Psi}}^{-1}_t\overset{p}{\rightarrow} {\bm{\Psi}}_t^{-1}{\bm{\Omega}}_t{\bm{\Psi}}_t^{-1}. \end{equation*}

Selection consistency

We next consider the consistency for the factor number estimators $\hat{r}_{\mathrm{IC}}$ and $\hat{r}_{\mathrm{rank}}$.

theoremLet $\delta_{NT}=(\sqrt{L_{NT}h^3})^{-1}+h^2$. Under Assumptions (ref)-(ref), if $P_{1,NT} \to 0$, $P_{1,NT}\delta_{NT}^2\to \infty,$ then $ P[\hat{r}_{\mathrm{rank}}=r_0]\to 1 $, as $ N, T \to \infty$.

remark\begin{enumerate}[(i)] • Theorem (ref) establishes the consistency of the rank estimator $\hat{r}_{\mathrm{rank}}$. The idea of the proof is as follows. For any $\mathbf{{\bm{\Lambda}}}^{r}$, $r_0<r\leq \bar{r}$, decompose $\mathbf{{\bm{\Lambda}}}^{r}=[\mathbf{{\bm{\Lambda}}}^{r,r_0},\mathbf{{\bm{\Lambda}}}^{r,-r_0}] $, where the submatrix $\mathbf{{\bm{\Lambda}}}^{r,r_0}$ collects the first $r_0$ columns of $\mathbf{{\bm{\Lambda}}}^{r}$, and $\mathbf{{\bm{\Lambda}}}^{r,-r_0}$ collects the remaining $r-r_0$ columns. To prove Theorem (ref), we establish that up to sign, it holds that \begin{equation} \|\hat{\bm{{\bm{\Lambda}}}}^{\bar{r}}-\bm{{\bm{\Lambda}}}^*_0 \|/\sqrt{N}=O_p(\delta_{NT}^{-1}), \end{equation} where ${\bm{\Lambda}}_0^*=[{\bm{\Lambda}}_0,\mathbf{0}_{N\times(\bar{r}-r_0)}]$. Then $\hat{\sigma}^{\bar{r}}_{N,j}={\sigma}_{Nj}+o_p(1) \overset{p}{\to}\sigma_{j}> 0$, for $j=1,\cdots, r_0$ by Assumption (ref) (ii), and $\hat{\sigma}^{\bar{r}}_{N,j}\leq \sum_{s=r_0+1}^{\bar{r}}\hat{\sigma}^{\bar{r}}_{N,s}=\|\hat{\bm{{\bm{\Lambda}}}}^{\bar{r},-r_0}\|^2/N=O_p(\delta_{NT}^{-2})$ for $j=r_0+1,\cdots,\bar{r}$. Therefore, $\hat{{\bm{\Lambda}}}^{{\bar{r}}^{\prime}}\hat{{\bm{\Lambda}}}^{\bar{r}}/N $ converges in probability to a matrix with rank $r_0$ at rate $\delta_{NT}^2$, which leads to the consistency of $\hat{r}_{\mathrm{rank}}$ naturally. • The proof of (ref) parallels that of Lemma 2 in CDG2021. The difference lies in that their proof requires $\rho_{\mathrm{min}}[{\bm{\Lambda}}_0^{\prime}{\bm{\Lambda}}_0]=\sigma_{Nr_0} \overset{p}{\to} \sigma_{r_0}>0$ as $N\to \infty$, which is not satisfied here, as $\rho_{\mathrm{min}}[{{\bm{\Lambda}}^*_0}^{\prime}{\bm{\Lambda}}^*_0]=0$. We address this discrepancy by showing that, for $\|{\bm{\Lambda}}^{r}{\bf{F}}^{r^{\prime}}-{\bm{\Lambda}}_0{\bf{F}}_0^{\prime}\|\leq \delta$ with $\delta>0$ sufficiently small and $r>r_0$, $\rho_{\mathrm{min}}[{\bm{\Lambda}}^{{r,r_0}^{\prime}}{\bm{\Lambda}^{r,r_0}}]$ is positively bounded below, together with leveraging properties of positive definite matrices. \end{enumerate}

theoremUnder Assumptions (ref)-(ref), if $P_{2,NT}\to 0$, $P_{2,NT}\delta_{NT}^2\to \infty$, then $P[\hat{r}_{\mathrm{IC}}=r_0]\to 1$, as $ N, T \to \infty$.

remark\begin{enumerate}[(i)] • Theorem (ref) establishes the consistency of the IC-based estimator $\hat{r}_{\mathrm{IC}}$. The proof of this result follows closely from BN2002. In particular, we consider two cases (i) $0<r<r_0$ and (ii) $r_0< r\leq \bar{r}$, respectively. For case (i) $0<r<r_0$, we prove that there exists $C>0$, such that $\mathbb{M}_{NT}(\hat{\bm{{\bm{\theta}}}}^{r_0})-\mathbb{M}_{NT}(\hat{\bm{{\bm{\theta}}}}^{r})\geq C+o_p(1)$. For case (ii) $r_0<r\leq \bar{r}$, we demonstrate $\mathbb{M}_{NT}(\hat{\bm{{\bm{\theta}}}}^{r_0})-\mathbb{M}_{NT}(\hat{\bm{{\bm{\theta}}}}^{r})=O_p(\delta_{NT}^{-2})$. Therefore, when $P_{2,NT}$ vanishes at a rate slower than $\delta_{NT}^2$, for any $r$ such that $0<r\neq r_0\leq \bar{r}$, IC$(r)$-IC$(r_0)$ will always be dominated by a positive term. Hence, $P[\hat{r}_{\mathrm{IC}}=r_0]\to 1$. • The conditions placed on the penalties in Theorem (ref) and Theorem (ref) are identical, therefore, the choices of $P_{1,NT}$ and $P_{2,NT}$ can be made at the same order. With the optimal bandwidth $h=O(L_{NT}^{-1/7})$, we obtain $\delta_{NT}=O(L_{NT}^{2/7})$. For $\hat{r}_{\mathrm{rank}}$, the choice \begin{equation*} P_{1,NT}=\hat{\sigma}^{\bar{r}}_{N,1}\cdot (L_{NT}^{4/7})^{-0.3}, \end{equation*} meets the rate requirement, which leads to high probability of correctly selecting the factor number as long as $\min\{N,T\}\geq 100$ in finite sample simulations later on. \end{enumerate} For $\hat{r}_{\mathrm{IC}}$, let $\mathbb{M}_{NT}=\frac{1}{NT}\sum_{i=1}^N\sum_{t=1}^TK_h(X_{it})$, then the choice \begin{equation*} P_{2,NT}=3/7\cdot (\mathbb{M}_{NT}(\hat{\bm{\theta}}^{\bar{r}})-\mathbb{M}_{NT})\cdot (L_{NT}^{4/7})^{-0.4}, \end{equation*} is desirable. Simulations indicate that $\hat{r}_{\mathrm{IC}}$ outperforms $\hat{r}_{\mathrm{rank}}$, and $\hat{r}_{\mathrm{IC}}$ can accurately select the number of factors in most settings even with sample size as small as $(N,T)=(60,60)$.

Simulation studies

This section conducts a set of Monte Carlo experiments to evaluate the finite sample performance of our proposed estimators, with a comparison to that of the estimators derived under AFM BN2002 and QFM CDG2021.

Data generating processes

The data are drawn from the following three-factor model:

equation*[equation* omitted — 122 chars of source]

where $f_{it},{\lambda}_{it}$ are independent $\mathcal{N}(0,1)$ variates. For the error term, we entertain three different error specifications.

\noindentS1 (Heavy-tailed errors). The error term ${e_{it}}^{\prime}s$ are independently drawn from $ t_\nu$, the students $t$-distribution with $\nu$ degrees of freedom, for $\nu=1,2,3$.

\noindentS2 (Dependent errors). Following BN2002, $e_{it}$ is generated according to

equation*[equation* omitted — 90 chars of source]

where ${v_{it}}^{\prime}s$ are independently drawn from $t_3$. The autoregressive coefficient $\rho$ reflects the serial correlation of $e_{it}$, while the parameters $\beta$ and $J$ reflect the cross-sectional dependence of $e_{it}$. In particular, the following three sets of error dependence parameters are considered.

(D1) Serially correlated errors: $\rho=0.2,\beta=0$.

(D2) Cross-sectionally correlated errors: $\rho=0,\beta=0.2,J=3$.

(D3) Serially and cross-sectionally correlated errors: $\rho=0.2,\beta=0.2,J=3$.

\noindentS3 (Skewed errors). Following YL2014, we entertain a mixture normal distribution for $e_{it}$,

equation*[equation* omitted — 90 chars of source]

where $\sigma\in \{2.6,3,3.4\}$. This distribution is skewed left with E$(e_{it})=0$, Mode$(e_{it})\approx 0.8$, Median$(e_{it})\approx \{0.53,0.57,0.60\}$, corresponding to the three values of $\sigma$. Additionally, as discussed in Example (ref), to ensure that AFM, MFM and QFM have the same representations, the error terms are accordingly normalized when obtaining the estimates for MFM ($e_{it}-\operatorname{Mode}(e_{it})$) and QFM ($e_{it}-\operatorname{Median}(e_{it})$).

Included in the comparison are the MFA estimator ${\hat{\bf{F}}}_M$ proposed in this paper, the PCA estimator ${\hat{\bf{F}}}_{P}$ studied by BN2002, and the QFA estimator at $\tau=0.5$, i.e., ${\hat{\bf{F}}}_{Q}^{0.5}$ proposed by CDG2021. To obtain ${\hat{\bf{F}}}_M$, we run the AMEM algorithm with two different sets of random starting parameters, and adopt the estimate that maximizes the objective function. For each error specification, we set $N,T \in \{60,100,200\}$. For the tuning parameters, we follow the discussion in Remark (ref) to set $h=c\cdot L_{NT}^{-1/7}$, $P_{1,NT}=\hat{\sigma}^{\bar{r}}_{N,1}\cdot (L_{NT}^{4/7})^{-0.3}$ and $P_{2,NT}=3/7\cdot (\mathbb{M}_{NT}(\hat{\bm{\theta}}^{\bar{r}})-\mathbb{M}_{NT})\cdot(L_{NT}^{4/7})^{-0.4}$. To evaluate the robustness of our MFA estimator to the bandwidth choice, we consider $c\in \{3,5,7\}$. Finally, we set the maximum number of factors $\bar{r}=8$ following BN2002.

Several commonly used criteria are employed to evaluate the performance of the estimators. First, to assess the performance of the estimated factors in capturing the space of the true factors, we consider the trace-ratio statistic adopted by BN2006 and CKKK2018, which measures the distance between the estimated factor space and the true factor space. Specifically, let $\hat{\mathbf{F}}$ and ${\mathbf{F}}_{0}$ denote the estimated and the true factor matrices, respectively. The trace-ratio statistic $\operatorname{tr}(\mathbf{\hat{F}})$ is defined as

equation*[equation* omitted — 248 chars of source]

It is seen that $\operatorname{tr}(\mathbf{\hat{F}})\in [0,1]$, and a larger $\operatorname{tr}(\mathbf{\hat{F}})$ indicates a smaller distance between the space spanned by $\mathbf{\hat{F}}$ and ${\bf{F}}_{0}$. Second, to evaluate the precision of model selection, we compute the average estimated number of factors and the frequency of correctly selecting the number of factors. All the results are obtained over $S=500$ replications.

Results

table[table omitted — 3,510 chars of source]

Table (ref) presents the factor estimation accuracy results under heavy-tailed errors S1. The main findings are summarized as follows. First, both the modal factor $\hat{\bf{F}}_M$ and the median factor ${\hat{\bf{F}}}_{Q}^{0.5}$ can effectively capture the true factor space for all three error choices. Second, our estimator ${\bf{\hat{F}}}_M$ is hardly affected by the choice of bandwidth and consistently outperforms the other two estimators across all parameter configurations. Third, ${\hat{\bf{F}}}_{P}$ is notably inferior to the other two estimators under $t_1$ and $t_2$ errors, but becomes comparable under $t_3$ errors. Overall, the performances of all three estimators tend to improve as the error distribution shifts from $t_1$ to $t_3$, and as both $N$ and $T$ increase.

The factor estimation accuracy results under dependent errors S2 and skewed errors S3 are presented in the Supplementary Materials. Additional to those discovered from Table (ref), there are several new findings which are summarized as follows. First, as the dependence in the error terms increases in S2, the performances of all three estimators deteriorate, but our estimator ${\hat{\bf{F}}}_M$ continues to provide the most accurate estimates in almost all cases. Second, the accuracy of ${\hat{\bf{F}}}_M$ and that of ${\hat{\bf{F}}}_Q^{0.5}$ remain nearly unchanged as the error terms become increasingly skewed in S3, while that of ${\hat{\bf{F}}}_P$ decreases notably. This finding demonstrates that ${\hat{\bf{F}}}_M$ and ${\hat{\bf{F}}}_Q^{0.5}$ exhibit robustness to the presence of skewness as compared to ${\hat{\bf{F}}}_P$. Third, the performance of ${\hat{\bf{F}}}_M$ is only slightly influenced by the choice of bandwidth in S2 and S3, even though the effect seems a bit larger than that in S1.

table[table omitted — 4,929 chars of source]

Table (ref) presents the results for model selection using the rank estimator $\hat{r}_{\mathrm{rank}}$ and the IC-based estimator $\hat{r}_{\mathrm{IC}}$ under S1, while those under S2-S3 are relegated to the Supplementary Material. In these tables, $\bar{\hat{r}}_{\mathrm{rank}}$ and $\mathit{Freq}_1$ denote the average number of factors selected and the frequency of selecting the correct number of factors by $\hat{r}_{\mathrm{rank}}$, respectively. $\bar{\hat{r}}_{\mathrm{IC}}$ and $\mathit{Freq}_2$ denote the corresponding estimates for ${\hat{r}}_{\mathrm{IC}}$. The main findings from Table (ref) are summarized as follows. First, although both $\hat{r}_{\mathrm{IC}}$ and $\hat{r}_{\mathrm{rank}}$ tend to underestimate the factor number with small sample sizes, the probability of correct estimation increases along with $N$ and $T$. In particular, both estimators can correctly estimate the true factor number when both $N$ and $T$ are as large as 200. Second, $\hat{r}_{\mathrm{IC}}$ is more precise with a smaller bandwidth across all scenarios, while $\hat{r}_{\mathrm{rank}}$ seems insensitive to the bandwidth setting. Third, $\hat{r}_{\mathrm{IC}}$ tends to outperform $\hat{r}_{\mathrm{rank}}$ in almost all cases. Fourth, $\hat{r}_{\mathrm{IC}}$ and $\hat{r}_{\mathrm{rank}}$ can provide accurate estimates as long as $\min\{N,T\}\geq 100$, while $\hat{r}_{\mathrm{rank}}$ achieves high accuracy in most settings even for sample sizes as small as $(N,T)=(60,60)$.

We finally evaluate the finite sample distributional behavior of the MFA estimators and check how closely the asymptotic distributions derived in Theorem (ref) approximate the finite sample distributions. Following Bai2003, we consider the following factor specification:

equation*[equation* omitted — 44 chars of source]

where $\lambda_i, f_t, e_{it}$ are $i.i.d$ $\mathcal{N}(0,1)$ for all $i,t$. The number of factors is set to $1$, and $(N,T)=(60,60),(100,100),(200,200)$. According to Assumption (ref), we set the bandwidth as $h=c \cdot T^{-1/12}, c\in\{3,5,7\}$.

Figure S.1 in the Supplementary Material displays standard normal density curve and the histogram for the estimates

equation*[equation* omitted — 119 chars of source]

at $t=T/2$ obtained from $1000$ repetitions.\footnote{We choose the sign of $\hat{f}_t$ such that $\hat{\bf{S}}=\mathbb{I}_{r_0}$} It can be seen that the standard normal density curve provides good approximations to the histograms in all cases. Therefore, the finite sample distribution of the standardized factor estimator is well approximated by standard normal distribution.

We further construct $95\%$ confidence intervals for the true factor process $\{f_{0t}, t=1,\cdots, 20\}$ in Figure S.2. The confidence intervals for the remaining time points are not presented for the sake of clarity. The solid curve in the middle of each plot represents the true factor process, and the dashed curves signify the estimated confidence intervals. It can be observed that, the confidence intervals are hardly sensitive to the choice of bandwidth and contain the true factors in most cases. In addition, they become narrower as $N, T$ grow.

Empirical applications

This section examines the prediction of U.S. inflation rate as studied by CDG2021, and explores the usefulness of the MFA factors as predictors. Specifically, we extract a set of common factors by MFA, PCA and QFA from a large panel of macroeconomic data, and evaluate the predictive power of each set of factors in forecasting the inflation rate. The data set we use is the FRED-QD dataset, which includes quarterly data for 211 U.S. macroeconomic time series from 1960Q1 to 2019Q2 ($N=211, T=238$), and is available from the Federal Reserve Economic Data database at \url{https://research.stlouisfed.org/econ/mccracken/fred-databases/}. Each series is first transformed to be stationary using MATLAB codes available on the FRED-QD data website. The transformed series is then demeaned and standardized to have zero mean and unit variance.

Let $y_{t}$ denote the realized value of U.S. inflation rate at period $t$.\footnote{We use CPI to measure the inflation rate. Let $x_t$ denote CPI at period $t$, then $y_t=\ln{x_t}-\ln{x_{t-4}}$.} The $s$-step-ahead forecasting model with factor-augmented predictors writes:

equation[equation omitted — 136 chars of source]

where ${\bf{F}}_t$ denotes the latent factors derived from the large data set. Based on the vector of estimated factors $\hat{\bf{F}}_t$, the least squares prediction of $y_{t+s}$ is obtained as:

equation*[equation* omitted — 128 chars of source]

where $\hat{\alpha},\hat{\beta}_{j}, \hat{\bm{\gamma}}$ are least squares estimates of the coefficients, and $\hat{p}$ is the optimal lag length determined by BIC.

For each time point $t$, we extract the MFA factors $\hat{\bf{F}}_{M,t}$, PCA factors $\hat{\bf{F}}_{P,t}$, as well as the QFA factors $\hat{\bf{F}}_{Q,t}^{\tau}$ at the $\tau$-th quantile from the data set. Then consider the specifications of the forecasting model: (i) $\hat{\bf{F}}_t=0$, which is the benchmark autoregressive (AR) model; (ii) AR+$\hat{\bf{F}}_{M,t}$; (iii) AR+$\hat{\bf{F}}_{P,t}$; (iv) AR+$\hat{\bf{F}}^{0.5}_{Q,t}$; (v) AR+$\hat{\bf{F}}_{P,t}+\hat{\bf{F}}^{\tau_1}_{Q,t}, \tau_1=0.9$; (vi) AR+$\hat{\bf{F}}_{P,t}+\hat{\bf{F}}^{\tau_2}_{Q,t}, \tau_2=0.99$. The model specifications in (v)-(vi) are taken from CDG2021.

A rolling window of $120$ quarters is adopted to estimate the coefficients and generate the rolling window forecasts. Within each window, the number of MFA, PCA, and QFA factors is determined by our IC-based estimator $\hat{r}_{\mathrm{IC}}$, the $PC_{p1}$ criterion of BN2002, and the rank-minimization estimator proposed by CDG2021, respectively. The maximum number of each kind of factor is set to $8$, and the maximum lag length is limited to 4, i.e., $p_{\mathrm{max}}=3$. Further, for the “AR+$\hat{\bf{F}}_{M,t}$" model, we set the bandwidth $h=c\cdot 120^{-1/7}$, and consider $c\in\{3,5,7\}$, as done in the simulations. Following CDG2021, the initial estimation period spans from 1960Q1 to 1989Q4 (120 quarters), and the forecast evaluation period is divided into two subperiods: the great moderation pre-crisis (1990Q1 to 2007Q2) and the financial crisis/recovery (2007Q3 to 2019Q2). The mean squared errors (MSEs) of forecasts from these models are calculated, and their relative MSEs (RMSEs) to that of the AR model are reported in Table (ref) for the full evaluation period and two subperiods.

table[table omitted — 4,195 chars of source]

The main findings in Table (ref) are summarized as follows. First, incorporating factors as predictors can improve the forecast accuracy of the benchmark AR model, and the “AR+$\hat{\bf{F}}_{M,t}$" model enjoys the best precision. Second, the “AR+$\hat{\bf{F}}_{M,t}$" model produces more accurate predictions than the AR model for both subperiods, which demonstrates the robustness of the “AR+$\hat{\bf{F}}_{M,t}$" model even in crisis scenarios. Third, the performance of the “AR+$\hat{\bf{F}}_{M,t}$" model is only slightly affected by the choice of bandwidth, and achieves optimal when the bandwidth constant $c$ is set as $5$.

To overcome the limitation of point forecasts that little is known regarding its accuracy, we next consider density forecasts of inflation rate based on (ref). Following CDG2021, we obtain these density forecasts using quantile regression (QR). Specifically, we first predict the conditional quantiles of the target variable by

equation*[equation* omitted — 146 chars of source]

for $\tau\in\{0.05, 0.25, 0.75,0.95\}$, where $\hat{\alpha}_{\tau},\hat{\beta}_{\tau,j}, \hat{\bm{\gamma}}_{\tau}$ are estimated by running QR of $y_{t+s}$ on $[1,y_t,\cdots, y_{t-p}, \hat{\bf{F}}_{t}]$. Given the predicted qunatiles $[\hat{q}_{0.05,t+s},\hat{q}_{0.25,t+s},\hat{q}_{0.75,t+s},\hat{q}_{0.95,t+s}]$, the predicted density of $y_{t+s}$ is constructed as the density of a skewed $t$-distribution by matching the predicted quantiles. Finally, the predictive score, i.e., the predicted density evaluated at the observed value of $y_{t+s}$, which measures the accuracy of density forecasts is obtained.\footnote{We refer readers to ABG2019 for further details of doing density forecasts, and to AC2003 for the definition and properties of the skewed $t$-distribution.} Note that higher scores indicate more accurate predictions. In line with before, we consider factor specifications (i)-(vi), and set $h=c\cdot 120^{-1/7}, c\in \{3,5,7\}$, $p=3$. Additionally, a rolling window of the most recent 120 observations is adopted to generate the out-of-sample density forecasts, and the evaluation period is from 1990Q1 to 2019Q2.

table[table omitted — 3,572 chars of source]

Table (ref) presents the average predictive scores of different models for predicting inflation rate over the whole evaluation period and two subperiods. It can be observed that, the “AR+$\hat{\bf{F}}_{M.t}$" model is hardly affected by the bandwidth setting, and its average predictive scores surpass those of the other models in the majority of cases. This indicates that the MFA factors provide valuable information for density forecasting of inflation rate.

Conclusions

This paper proposes a modal factor model to extract factors influencing the conditional mode of the distribution of the observables. An AMEM algorithm is developed to obtain the factor and loading estimators and two model selection methods are introduced for selecting the number of factors. The asymptotic properties of the proposed estimators are established and numerical results demonstrate the nice finite sample performance of these estimators in both simulations and empirical applications to forecasting macroeconomic variables.

The current paper can be extended in several directions. First, further investigation is needed on how to select the bandwidth in MFM, and the cross-validation method could be considered following CGTW2016. Second, the conditional cross-sectional independence assumption on the error terms may be further relaxed, as supported by the simulation results. Third, our results can be extended to cover the case where the conditional mode has a non-linear factor representation Wang2024,MT2023b. Fourth, it would be interesting to further consider inference on the role of modal factors in factor-augmented models. Fifth, it is possible to extend the static MFM to dynamic MFM by allowing the factor loadings to be time-variant SW2017 or by including lagged factors FHLR2005. These issues involve new technical challenges and deserve separate future efforts.

Acknowledgments

The authors would like to thank seminar and conference participants at Beihang University, Nanjing Audit University and Tsinghua University for comments that help improve the paper. The authors thank the partial support from National Natural Science Foundation of China (Grant 72425009, 72073002), the Center for Statistical Science at Peking University, and Key Laboratory of Mathematical Economics and Quantitative Finance (Peking University), Ministry of Education.