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.
83,065 characters · 11 sections · 40 citation commands
Factor-Based Imputation of Missing Values and Covariances in Panel Data of Large Dimensions
Keywords: risk management, covariance structure, matrix completion, incomplete data
JEL Classification: C1, C2.
\setcounter{page}{0} \thispagestyle{empty}
\baselineskip=18.0pt
Missing data is a problem that empirical researchers frequently encounter. For example, data can be missing due to attrition in longitudinal surveys such as the PSID and SIPP, to IPOs and bankruptcies in the case of stock returns, and to staggered releases of data by different institutions. Data sampled at a lower frequency will have missing values in an analysis that involves higher frequency data. Whatever is the reason, the way we handle missing data is not innocuous. Dropping rows or columns necessarily entails loss of information, and while the EM algorithm allows imputation of missing values, it is designed for low dimensional data and may not be easily scalable.
In recent years, progress has been made for situations when the missing values occur in large panels of data with a low rank component. Assuming incoherence conditions to ensure that the low rank structure is non-degenerate as in candes-recht:09, the machine learning literature has developed regularization based algorithms to solve matrix completion problems of the Netflix challenge type. On the econometrics front, recent work by su-missing xiong-pelger:19, and baing-miss:jasa propose different implementations that permit the entire common (low rank) component to be consistently estimated if the factor structure is strong. While the machine learning literature provides worst-case error bounds, the econometrics literature provides the distribution theory that is needed for inference, and is also the approach that this paper takes.
We introduce a Tall-Project estimator (or TP for short) for imputing missing values in a panel of data with $T$ rows and $N$ columns, where $N$ and $T$ are both large. The factors are estimated from a tall\ block consisting of complete data for $N_o\le N$ units, while the loadings are obtained from $N$ time series projections with series-specific sample size. We obtain two main results. First, under certain assumptions on the factor structure in the different blocks, we show that the tp estimates are consistent and asymptotically normal for every entry of the low rank component though the convergence rate is series specific. And while iteration is not needed for consistent estimation, one re-estimation can improve the convergence rate. Because of missing data, stronger conditions are needed for the factor estimates to be treated as known in factor-augmented regressions. \footnote{An R package that implements tw and \textsc{tp} is available for download at \url{htpps://github.com/cykbennie/fbi.}}
Our second result concerns estimation of covariances from imputed data. Covariances play an important role in portfolio analysis and in structural equation (covariance structure) modeling. While the imputed data is unbiased for its mean, the variance of the imputed data is biased because the idiosyncratic noise associated with the missing observations are set to zero. To remedy this problem, we repeatedly overlay the first step estimate of missing values with resampled idiosyncractic residuals before computing sample covariances. An average of these imputed covariances is then used to estimate the objects of interest. In simulations calibrated to CRSP data from 1990-2018, resampling from the own residuals (ie. without pooling) yields risk measures that compare favorably and sometimes outperform the ones based on a factor-based covariance estimator.
In what follows, let $X$ be a $T\times N$ panel of data, $X_i=(X_{i1},\ldots X_{iT})^\prime$ be a $T\times 1$ vector of random variables and $X=(X_1,X_2,\ldots, X_N)$ be a $T\times N$ matrix. We use $i=1,\ldots N$ to index cross-section units and $t=1,\ldots T$ to index time series observations. In practice, $X_{i}$ is transformed to be stationary, demeaned, and is often standardized. It is assumed that the normalized data $Z=\frac{X}{\sqrt{NT}}$ has singular value decomposition (svd) \[ Z=\frac{X}{\sqrt{NT}} =U DV{^\prime }=U_r D_r V_r{^\prime }+ U_{n-r} D_{n-r} V_{n-r}{^\prime }\] where $ D_r$ is a diagonal matrix of $r$ singular values ordered such that $ d_1\ge d_2\ldots \ge d_r$, while $U_r, V_r^\prime$ are the corresponding left and right singular vectors respectively. Analogously, $U_{n-r}$ and $V_{n-r}^\prime$ are $n\times (n-r)$ matrices of left and right eigenvectors associated with $d_{r+1},\ldots, d_n$. The low rank component $U_rD_rV_r^\prime$ can be defined without probabilistic assumptions..
We consider a factor model defined as
where $F$ is a $T\times r$ matrix of common factors, $\Lambda$ is a $N \times r$ matrix of factor loadings, and $e$ is a $T\times N$ matrix of idiosyncratic errors. We will let $(F^0,\Lambda^0)$ be the true values of $(F,\Lambda)$. The common component $C^0=F^0\Lambda^{0{^\prime }}$ has reduced rank $r$ because $F^0$ and $\Lambda^0$ both have rank $r$. We will estimate the factors and loadings by the method of static asymptotic principal components (APC).\footnote{In this paper, we rely heavily on the results given in stock-watson-jasa:02,baing-ecta:02,bai-ecta-03,baing-ecta:06.} The normalization $\frac{F{^\prime } F}T=I_r$ gives the APC estimates that are unique up to a column sign: \[ (\tilde F,\tilde\Lambda)=(\sqrt{T} U_r, \sqrt{N} V_r D_r). \] Since $r$ can be consistently estimated, we proceed as though $r$ is known.
We consider a large dimensional approximate factor model in which $T$ and $N$ are large and which satisfies the following assumptions:
\paragraph{Assumption A:} There exists a constant $M<\infty$ not depending on $N,T$ such that
Assumption A is used in bai-ecta-03 to develop an inferential theory for large dimensional factor analysis when $X$ is completely observed. The moment conditions in (b) ensure that the factor structure is strong and can be separated from the idiosyncratic errors which can be weakly correlated, both in the time and the cross-section dimensions. Let $\mathbb D_{r}^2$ and $ \mathbb V_{r } $ be the eigenvalues and eigenvectors of the $r\times r$ matrix $\Sigma_{\Lambda}^{1/2}\Sigma_F\Sigma_{\Lambda}^{1/2}$, respectively. As shown in bai-ecta-03, $\plim_{N,T\rightarrow\infty} D_r^2 = \mathbb D_r^2$ and $\plim_{N,T\rightarrow\infty} \frac{\tilde F{^\prime } F^0}{T}=\mathbb Q_r$, where $\mathbb Q_r=\mathbb D_{r} \mathbb V_{r}\Sigma_{\Lambda}^{-1/2}$. The following properties of the APC estimator are given in bai-ecta-03:
In what follows, we will obtain results for the factor estimates when $X$ has missing values.
A panel is said to be complete if it does not have missing values. In practice most data panels have missing values. The data are said to be {\em missing completely at random} (MCAR) if missingness is independent of the data whether or not they are observed. As the observed and missing data are drawn from the same underlying distribution and hence have no systematic differences, MCAR does not create bias. This is not the case with {\em missing at random} (MAR), though the systematic differences can be explained by other observables. For example, missing values in financial statements of smaller firms may arise because smaller firms are less regulated, in which case the missing data propensity can be related to market capitalization which is observed for publicly traded firms. But the issue remains that the observed values under MAR do not form a random sample.
To overcome the bias when the data are missing at random, one approach is to re-weigh the data (such as 'inverse probability weighting' used in survey design). A second approach is to rebalance the panel by amputation which can take one of three forms: (a) listwise (complete case) deletion under which the entire row with missing values is eliminated; (b) variable deletion under which a series with one or more missing values will be eliminated; and (c) pairwise deletion where only the cases with missing data involved will be deleted. Though the solution is simple and does not entail modeling assumptions, amputation leads to significant information loss. Furthermore, the resulting balanced sample may not be representative of the population and may yield biased estimates.
Another way to rebalance the panel is imputation. The simplest procedure is to replace the missing values with the sample mean or median. A more sophisticated procedure is to impute from the unconditional sample first and second moments via the EM algorithm. schneider:01 proposes a regularized EM algorithm that is quite popular in analysis of climate data. Imputation on the basis of unconditional moments does not require a model. A fully specified approach is to construct the likelihood using incomplete observed data and iteratively solve for the parameters of the conditional mean function. The estimates would be more efficient if the parametric assumptions were correct. These methods are designed for imputing missing values in a small number of variables/predictors. See robins-rotnitzky-zhao:95 li-shen-li-robins:13 and raghunathan:04 for a review.
Factor models provide a simple framework for imputing missing values of many variables. banbura-modugno:14, jungbacker-koopman, jkv:11 consider state-space modeling of a strict factor (i.e., a diagonal idiosyncratic error variance matrix) with missing data and uses non-linear filters to compute the likelihood as $F_t$ and $\Lambda_i$ are both random. giannone-reichlin-small initializes the missing values with estimates from the balanced panel and uses the Kalman filter to perform updating. stock-watson:handbook-16 discusses the issues with state space estimation of factor models with missing data. For large dimensional approximate factor models, stock-watson-di-wp suggests to fill missing values in $X$ with the APC estimates of the common component. These are all EM algorithms that use the factor structure to evaluate the conditional mean in the E-step, and principal components estimation in the M step.
Though many factor-based imputation methods have been used for some time, the theoretical properties of the imputed estimates are studied only recently in independent works by su-missing, xiong-pelger:19, and baing-miss:jasa. All three exploit the strong factor structure in a large panel setting to estimate the factors and the loadings from incompletely observed data. su-missing focuses on the case of missing at random and analyzes the EM estimator considered in stock-watson-di-wp. After initializing the missing values to zero, the factors are estimated from data reweighted by the frequency of no missing values. Though these estimates are consistent, they are not asymptotically normal without further iteration. xiong-pelger:19 also re-weights the data in principal components estimation of the loadings, but uses cross-sectional regressions to estimate the factors at every $t$. Similar to su-missing, the estimates are consistent and asymptotically normal when iterated until convergence. but the `all purpose' estimator proposed in xiong-pelger:19 allows for mechanisms other than missing at random. Though robustness comes at the cost of larger variances, efficiency can be improved if stronger assumptions on missingness are made.
While su-missing and xiong-pelger:19 reweigh the data in APC estimation, baing-miss:jasa implicitly re-oganizes the data into blocks. This is partly motivated by the concern that the missing at random assumption may at times be inappropriate. Examples of microeconomic applications are discussed in abdik. Furthermore, for macroeconomic panels such as FRED-MD where the data are collected from different sources, assuming that the missing values are due to the same missing mechanism seems restrictive.
Consider the left panel of Figure (ref). This systematic (non-random) missing pattern can arise if, for example, weekly data for related variables (such as inventories and orders) are released in the same week of every month. The right panel shows a case of `semi-block missing'. This can arise if, for example in a survey, respondent $N$ dropped out early, respondent $N-1$ did not respond in period $T-1$, while $N_o+1$ did not respond in period $T$. The observations can be missing for different reasons likely unknown to the researcher.
While the missing data mechanism can be difficult to verify, we do observe the missing data pattern, and organizing the data into blocks offers a new way of thinking about imputation. The insight in baing-miss:jasa is that the quantities needed for imputation can be obtained from two completely observed blocks. Their tall-wide (or tw for short) procedure uses the tall block consisting of complete data for $N_o<N$ units to estimate the factors, and a wide block in which data for all units are available over $T_o<T$ periods to estimate the factor loadings. To align the space spanned by the estimated factors and loadings, one then estimates a new rotation matrix by regressing the $N_o\times r$ matrix of loadings estimated from the tall block on a sub-block of $N\times r$ matrix of loadings estimated from the wide block. The procedure produces consistent and asymptotically normal estimates without iteration. Re-estimation using the completed (imputed) data will improve efficiency of the estimates in the balance block where the \textsc{tall} and \textsc{wide} blocks intersect to the fastest rate possible, which is the rate obtained in the complete data case.
The tw algorithm takes as given that the tall and wide blocks are 'big enough' for consistent estimation. Because the size of the missing data block is defined by number of time periods and units with complete cases, the size of the missing block is the same for both examples in Figure (ref). The algorithm uses information available efficiently when the missing pattern is homogeneous. But for the example in the right panel, more information could have been used.
More concerning are situations when a few series significantly reduce $N_o$ and $T_o$. As an example, consider Figure (ref). The left panel shows a case of reverse monotone missing. In multi-country panels, this pattern will arise if data for less developed countries are only available at a later date than developed countries. The right panel shows a case where the sample size is much smaller for unit $N$ than for units $N-1$ and $N-2$. In both examples, the missing block can contain many observed values.
In this section, we will develop a new procedure to obtain consistent estimates without iteration as in tw, but can accommodate flexible missing data patterns while making more efficient use information available. The point of departure is to partition the $T\times N$ matrix $X$ in a different way. We assume that there are $N_o$ series with no missing values so there is a $T\times N_o$ block labeled tall, and a $T\times (N-N_o)$ block of partially observed data labeled incomplete. There is no need to reorganize the data in practice, but Figure (ref) shows the data with the $N_o$ variables ordered first to help visualize the idea.
It will be helpful to define two locator sets as follows:
Let $N_{o_t}$ be the number of units observed at time $t$ and $T_{o_i}$ be the number of rows observed for series $i$. Then $J^t$ is a $N_{o_t}\times 1$ vector that keeps track of the units observed at $t$, and $J_i$ is a $T_{o_i}\times 1$ vector that keeps track of the periods that unit $i$ is observed. For a set $\mathcal A$, let $|\mathcal A|$ be its cardinality. In this notation, $N_o=| \cap_t J^t|$, $T_o=| \cap_i J_i|$, and $T_{o_i}=T$ for every unit in tall.
\paragraph{Algorithm Tall-Project (TP):}
Step (i) of tp is the same as tw since both procedures estimate $F$ from the tall block. If such a block does not exist, then neither tw or tp will be appropriate. Fortunately, for macroeconomic panels, a tall\ block is often available. Algorithm \textsc{tp} accommodates staggered and irregular patterns of missingness by estimating $\Lambda_i$ using a customized sample for each series. The two algorithms should be numerically identical if $ T_{o_i}=T_o$ for all $i$, but since $N_{o_t}\ge N_o$ and $T_{o_i}\ge T_o$, \textsc{tp} will utilize more information in general. The difference can be significant if $T_o$ is dictated by a few series with many missing observations. Note, however, that \textsc{tp} estimates the $H\Lambda_i$ matrix directly, while \textsc{tw} estimates $H$ and $\Lambda_i$ separately. Furthermore, the number of factors in Algorithm \textsc{tp} is determined by the \textsc{tall}\; block so Algorithm \textsc{tw} is more flexible in this regard.
We will need additional assumptions to analyze the properties of tp.
\paragraph{Assumption B:} $\frac{\sqrt{N}}{\min\{N_o,T_o\}}\rightarrow 0$ and $ \frac{\sqrt{T}}{\min\{N_o,T_o\}}\rightarrow 0$ as $N\rightarrow\infty$ and $T\rightarrow\infty$.
\paragraph{Assumption C:} There exists $M<\infty$ such that for all $N$ and $t$, $ \frac {N_{o_t}(N - N_{o_t})} {N N_o} \le M$.
\paragraph{Assumption D:}
For each $i$, $\Big(\frac 1 {T-T_{{o_i}}} \sum\limits_{s\notin J_i} F_s^0 F_s^{0\prime}\Big) \Big( \frac1 {T_{{o_i}}} \sum\limits_{s\in J_i} F_s^0 F_s^{0\prime} \Big)^{-1} \smash{\mathop{\longrightarrow}\limits^p} I_r,$ $\frac 1 {\sqrt{T_{o_i}} } \sum\limits_{s\in J_i} F_s^0 e_{is} \smash{\mathop{\longrightarrow}\limits^d} N(0,\Phi_i).$ For each $t$, $\Big(\frac 1 {N-N_{o_t}} \sum_{k\notin J^t} \Lambda_k^0\Lambda_k^{0\prime}\Big) \bigg(\frac{1}{N_o} \sum_{k \in \cap_s J^s} \Lambda_k^{0\prime}\Lambda_k^0\bigg)^{-1}\smash{\mathop{\longrightarrow}\limits^p} I_r$, $\frac{1}{\sqrt{N_{o_t}}} \sum_{k \in J^t} \Lambda_k^0 e_{kt}\smash{\mathop{\longrightarrow}\limits^d} N(0,\Gamma_t)$.
Assumption B essentially puts a lower bound on the observed number of rows and columns. Assumption C puts an upper bound on the extent of missing values at any $t$ and can be understood as a noise-to-signal constraint. Assumption C allows $N_{o_t}/N\rightarrow 0$, and especially $N_o/N\rightarrow 0$. For example, if $N_{o_t}=N_o$, then the ratio is bounded by 1. For the central limit theorems in Assumption D to hold, $T_{o_i}$ should be independent of $(F^0_{s}, e_s)$, and $N_{o_t}$ independent of $(\Lambda_k^0,e_k)$. The assumption is satisfied if the missing data are unrelated to the intrinsic properties of the series. For example, missing data due to mergers are allowed if the events are unrelated to $F_t^0$. However, missing data due to bankruptcies would not be allowed if the bankruptcy probability depends on $F_t^0$. These conditions are stronger than stationarity and need to be justified on an application by application basis. All factor-based imputation procedure requires a similar assumption because without some commonality between the observed and missing units, imputation would not be possible.
Part $i$ and $ii(a) $ are implied by Lemma (ref) because $\tilde F_t$ and $\tilde \Lambda_i$ are estimated from the tall block, which is a complete matrix of $T \times N_o$ dimension. Result (b) of part ii arises because the factor loadings for unit $i$ in the incomplete block are estimated from a regression with $T_{o_i}$ observations. Part (iii) shows that each $\tilde C_{it}$ has its own convergence rate that depends on $T_{o_i}$, the number of observed entries in $X_i$. The result is based on the fact for those with $X_{it}$ that are missing,
The error due to estimating $F$ from the tall block appears as the first term on the right hand side and is denoted $u_{it}$, while the error due to estimating $\Lambda_i$ by projections is summarized by the second term and is denoted $v_{it}$. The term $r_{it}=O_p(\delta^{-2}_{T_o,N_o})$ uniformly in $i$ and $t$ represents high order errors when estimating the factor and factor loadings. These quantities depend on $N$ and $T$ but the notation is suppressed for simplicity. The representation (shown in the Appendix) is useful in understanding how $N$ and $T$ affect factor-based imputation.
As $\tilde F$ is estimated using the tall block alone, re-estimation of the factors and the loadings from the completed (imputed) matrix $\tilde X$ appears desirable. But imputation error in $\tilde X$ must be taken into account. Using ((ref)), we have
A quantity that will play a role in APC estimation from $\tilde X$ is
Under Assumption C, $B^i_{t}$ is bounded. For units in the tall block (i.e., $i \in \cap_s J^s)$, $B_t^i$ is roughly inflated from $\frac {N_{o_t}} N I_r$ by $N_{o_t}(N-N_{o_t})/(NN_o)$ which can be thought of as the noise to signal ratio..
The proposition establishes the convergence rate of the factors and loadings constructed from $\tilde X$. There are three changes due to re-estimation. First, there is only one (instead of many) rotation matrices for the factor loadings. This is a consequence of the fact that the factors are now estimated from $\tilde X$ in its entirety, instead of sub-blocks. Second, re-estimation generates efficiency gains. The convergence rate for $\tilde F_t^+$ is improved from $\sqrt{N_o}$ to $\sqrt{N_{o_t}}$ or to $\sqrt{N}$. Though the rate for $\tilde \Lambda^+$ is unchanged, the convergence rate of $\tilde C^+_{it}$ is now $\min(\sqrt{N_{o_t}},\sqrt{T_{o_i}})$ instead of $\min(\sqrt{N_o},\sqrt{T_{o_i}})$. If $N_{o_t}=N$ and $T_{o_i}=T$ for a given pair $(i,t)$, the convergence rate for $\tilde C_{it}^+$ is $\min(\sqrt{N},\sqrt{T})$, the best rate possible (the same as in complete data). Third, the asymptotic variance of $\tilde C_{it}^+$ depends on the $r\times r$ matrix $B_{t}^i$. From this result, a 95% confidence interval for $C^0_{it}$ of $\tilde C^+_{it} \pm 1.96 \tilde {\text{se}}(\tilde C^+_{it})$ where $\tilde{ \text{se}}(\tilde C^+_{it})=\sqrt{\frac{\tilde {\mathbb V^+_{it}}(N_{o_t},T_{o_i})}{\min(N_{o_t},T_{o_i})}}$. The prediction interval for a missing $X_{it}$ is $\tilde X^+_{it}\pm 1.96 \tilde {\text{se}}(\tilde X^+_{it})$ where $\tilde{ \text{se}}(\tilde X^+_{it})=(\tilde \sigma_{ei}^2+\tilde{ \text{se}}(\tilde C^+_{it})^2)^{1/2}$, $\tilde\sigma_{ei}^2=\frac{1}{T_{o_i}}\sum_{s \in J_i} \tilde e_{is}^2$ is an estimate of the variance of $e_{it}$.
Consider the infeasible regression with observed covariates $W_t$ and a latent predictor $F_t$:
where $z_t=(F_t^\prime, W_t^\prime)^\prime)$ and $\delta=(\alpha^\prime,\beta^\prime)^\prime$. Suppose that $F_t$ is replaced by an estimate based on a small number of predictors. It is known that in general, the sampling uncertainty in a generated regressor will inflate standard errors in subsequent regressions. A useful result in baing-ecta:06 is that if $F$ is estimated by APC from a completely observed panel $X$, then $\tilde F$ can be used in a second step regression without the need for standard error adjustments if certain conditions on the sample size are satisfied. Our foregoing analysis suggests that even if $X$ has missing values, $\tilde F^+$ estimated by tp is still consistent for the space spanned by $F$ albeit at a slower convergence rate.
Part (i) is based on Lemma 3 of baing-miss:jasa obtained for tw which can be used as a worse case rate for tp since $T_o$ in tw cannot exceed $\min_i T_{oi}$ in tp. A consequence of (i) is that
These results are used to obtain (ii). In particular, let $\mathbf S_{\hat z\hat z}=\frac{1}{T}\sum_{t=1}^T \hat z_t\hat z_t'$. Then
provided $\frac{\sqrt{T}}{N_o}\rightarrow 0$, $\frac{\sqrt{T}}{T_o}\rightarrow 0$. If, in addition, $\frac{1}{\sqrt{T}} \sum_{t=1}^T z_t\epsilon_{t+h}\smash{\mathop{\longrightarrow}\limits^d} N(0,\mathbf \Sigma_{zz,\epsilon})$, the asymptotic variance of $\hat\delta$ is of the sandwich form $\mathbf A^{-1} \mathbf B \mathbf A^{-1}$ with $\mathbf A=\mathbf J\mathbf \Sigma_{zz}\mathbf J$ and $\mathbf B = \mathbf J\mathbf \Sigma_{zz,\epsilon}\mathbf J'$, where $\mathbf J$ is the probability limit of $\mathbf J_{NT}$. Simplifying yields (ii). The result implies that as in baing-ecta:06, the estimated factors can be treated as though $F$ were observed albeit under more stringent conditions than the complete data case which only requires that $\sqrt{T}/N\rightarrow 0$. For macroeconomic analysis when the factors are estimated from panels with large $T_o$ and $N_o$, one can still expect $\hat\delta$ to be precisely estimated up to a rotation.
Consider portfolio analysis where $X$ is a matrix of returns and $\Sigma_X$ is its covariance. The population weights of a minimum variance portfolio are given by $ w_p = \frac{ \boldsymbol{ \Sigma_X^{-1} 1}}{ \boldsymbol{ 1\prime \Sigma_X^{-1} 1} }$ where $\boldsymbol{1}$ is the unit vector. If $X$ was completely observed and $x_i=X_i-\bar X_i$, the sample moments (SM) estimator of the $N\times N$ matrix is
It is well known that the sample covariance matrix is singular when $N>T$ and has large sampling uncertainty when $N$ and $T$ are large. Nonetheless, when the data admit a factor structure so that the decomposition $\mathbf \Sigma_X=\mathbf \Lambda\mathbf \Sigma_F \Lambda^\prime+\mathbf \Sigma_e$ holds, a factor-based covariance estimator can be defined as
where $\hat{\mathbf \Sigma}_{\tilde e}$ is the sample covariance of $\tilde {\mathbf e} $, which is never exactly diagonal even when $\mathbf \Sigma_e$ is diagonal. As argued in chamberlain-rothschild, an approximate factor model that allows for some correlation in the idiosyncratic errors is usually a better characterization of returns data. To distinguish a strict from an approximate factor structure, let $\mathbf \Psi_e$ be a diagonal matrix whose $i$-th entry is $\mathbf E[e_{it}^2]$ and which can be consistently estimated by
The strict-factor (SF) covariance estimator is defined as \[ \mathbf {\tilde \Sigma}_{ X}= \tilde \Lambda \tilde{\mathbf \Sigma}_F \tilde \Lambda^\prime+ \tilde { \mathbf \Psi}_{\tilde e} \tag{SF} \] Consistency of this estimator is studied by fan:11, among others, for large $N$ and large $T$, assuming that $X$ is completely observed.
We now turn to the incomplete data case. Missing data presents a challenge for risk management because the portfolio weights depend on $\mathbf \Sigma_X$. Variable deletion is not an option when the variances and covariances of the missing variables are the objects of interest. Practitioners often resort to listwise deletion of the sample data, deleting the entire record from the analysis if any single value is missing. The covariance estimates can be unstable when the reduction in the sample size is sufficiently large. Though it is possible to calculate the sample covariance matrix with pairwise-complete observations, there is no guarantee that the resulting matrix will be positive definite. This is problematic because a singular covariance matrix must have at least one eigenvalue that is zero, implying that it would be possible to construct an eigenportfolio that has zero volatility (risk) by using the corresponding eigenvector as the weighs, making it possible to have an infinite ex-ante Sharpe ratio. Not only is this unrealistic, the hedges implied by the eigenportfolio weights are spurious and would fail out of sample. For these reasons, singular covariance matrices are of little use in portfolio construction in practice. Most risk applications require an estimate of a non-singular covariance matrix. The EM algorithm is one possibility, but it is likelihood based and requires parametric assumptions.
Unfortunately, successful factor-based imputation of the level of the data is not enough for precise estimation of their covariances. Even if the strict factor structure is correctly imposed, both the sample-covariance estimator and factor-based covariance estimator based on the imputed data will not be consistent. The problem with both estimators is that $\tilde e_{it}$ is set to zero when $X_{it}$ is not observed, so the sample variance of series $i$ is biased whenever the series has missing values. We can replace $\tilde{\mathbf \Psi}_{\tilde e}$ in ((ref)) by \[ \hat {\mathbf \Psi}_{\tilde e} =\text{diag}\bigg(\frac{1}{T_{o_i}}\sum_{t\in J_i} \tilde e_{it}^2,i=1,2,...,N\bigg) \] to recognize that some $\tilde e_{it}$ are zero, and the estimate $\hat\Psi_{\hat e,ii}=\frac{1}{T_{o_i}}\sum_{t\in J_i} \tilde e_{it}^2$ is consistent for $\mathbf \Psi_{e,ii}$ if $T_{{o_i}}$ is sufficiently large. Define the strict-factor adjusted (SFA) covariance matrix estimator as \[ \mathbf {\hat \Sigma}_{\tilde X}= \tilde \Lambda \tilde{\mathbf \Sigma}_F \tilde \Lambda^\prime+ \hat { \mathbf \Psi}_{\tilde e} \tag{SFA} \]
The proof is given in the Appendix.
This subsection considers a sample covariance estimator as an alternative to the SFA covariance estimator just described. Motivated by the fact that the variance of an imputed variable is downward biased, we use a second imputation to rectify this problem. \paragraph{Algorithm Residual Overlay}
Let $\tilde {e} = \tilde X - \tilde C$ where $\tilde X$ is produced by Algorithm tp. For residual resampling scheme $j$ ($j=1,2,3,4)$ to be discussed below, repeat for $s =1,\ldots S$:
The term residual overlay is motivated by the fact that errors are added to $\tilde X$. The algorithm injects randomness to $\tilde X_{it}$ whenever $X_{it}$ is missing. It uses stochastic simulations to compensate for the lost variability, an idea briefly considered in enders:10 in a fixed $N$ or $T$ setting. Even though each $\hat \Sigma_{\hat X}(s,j)$ is large in dimension and has rank $\min(N,T)$, in our experience, the average $\bar{\bf \Sigma}_{\tilde X}$ is always full rank, making it a viable alternative estimator in the case when $N>T$.
The method of multiple imputation typically estimates the object of interset each time an imputed set of data is obtained, and produces as output the average over the multiple estimates. In contrast, we average the multiple imputed covariances to obtain $\bar{\mathbf \Sigma}_{\hat X}$ and compute the object of interest (for example, portfolio weights) from it. Averaging has the advantage of reducing sampling uncertainty especially when $N$ and $T$ are both large.
\paragraph{Step a: Four Sampling Schemes}
All four methods assume that the returns data are covariance stationary.\footnote{If the variances and correlations change over time, the covariance matrix can be calculated by estimating GARCH models, for example.} Methods sm1 and sm2 are non-parametric, while method sm3 and sm4 calibrate the first two moments to the estimated errors associated with the observed data. Methods sm1 and sm3 are better suited for homogeneous data and since $\hat e_{it}$ are sampled from the stacked up vector of estimated errors. Methods \textsc{sm2} and \textsc{sm4} are better suited for heterogeneous data as they sample from the errors of the corresponding series. Methods \textsc{sm1} and \textsc{sm2} do not make distributional assumptions about the errors. As many financial time series tend to have fat tails and skewed distributions, other parametric distributions can be used in place of the normal distribution in \textsc{sm3} and \textsc{sm4}. The four sampling schemes do not account for cross-correlation in the errors, which implicitly shrinks $\mathbf \Sigma_e$ towards zero. This can be desirable when estimating high dimensional covariance matrices.
Methods sm+1, sm+2, sm+3, sm+4 are likewise defined with $\tilde e_{it}$ replaced by $\tilde e_{it}^+=\tilde X_{it}-\tilde F_t^{+\prime} \tilde \Lambda_i^+$, where $\tilde F^+$ and $\tilde \Lambda^+$ are obtained from the principal components of the imputed data $\tilde X$. In principle, we can also compute a SF estimator $ \hat{\mathbf \Sigma}_{ \hat X}(j)=\tilde{\mathbf \Lambda}\tilde{\mathbf \Sigma}_{ F}\tilde{\mathbf \Lambda}'+\frac{1}{S}\sum_{s=1}^S \hat{\mathbf \Psi}_{\hat e}(s,j) $ with $ \hat{\mathbf \Psi}_{\hat e_{ii}}(s)=\frac{1}{T} \sum_{t=1}^T\hat e_{it}^2(s,j)$, but there seems no advantage over SFA which already gives a good estimate of the variance of $e_{it}$ when $T_{oi}$ is sufficiently large. This is indeed the case in simulations, and hence not considered.
Though injecting noise to imputed values is not new, overlaying the noise to $\tilde X$ imputed using SM appears to be new. The resampling involves bootstrap draws. Bootstrapping large dimensional matrices is not a trivial exercise. Our set up is different, as we only need to bootstrap the noise corresponding to the missing data, which is much smaller in dimension. However, we have to pre-estimate the factors and the loadings, and bootstrapping the factor estimates is also not a trivial problem as seen in goncalves-perron:20. The proof of consistency of this estimator is thus quite delicate and is beyond the scope of this analysis. However, we can use simulations to evaluate if the idea holds promise to be worthy of further investigation.
In this section, we consider four experiments: the first and second assess the accuracy of the asymptotic approximations for $\tilde C_{it}$ and $\tilde C_{it}^+$ given in Lemma 2 and Proposition 1. The third evaluates the residual overlay procedures assuming a strict factor structure, while the fourth uses observed SP 500 returns as complete data to mimic an approximate factor structure.
The first experiment compares the performance of tp with and without re-estimation. The design of the monte-carlo is identical to the one used in baing-miss:jasa. Data are generated from $F\sim N(0,D_r)$ and $\Lambda\sim N(0,D_r)$ with $r=2$, the diagonal entries in $D_r$ are equally spaced between 1 and $1/r$, and $e_{it}\sim N(0,1)$. For each replication, the error in estimating $C_{it}$ is computed for four locations of $(i,t)$. As benchmark, we consider the infeasible case of complete data labeled complete. Also reported are results for the EM algorithm in stock-watson:handbook-16 that uses the estimates from the balanced panel as initial values. The algorithm repeatedly regresses $X$ on $F$ and then $X$ on $\Lambda$ till convergence, but the converged estimates may not be mutually orthogonal. We consider three versions of each estimator: one applied to the raw data $X$, one to the demeanend data, one to the standardized data. These are labled TP2, TP1, TP0 in the tables reported. The means and standard deviations are computed using the observations available for each series.
Table (ref) reports the root-mean-square-error for four chosen $(i,t)$ pairs, one in each of the four blocks: tall is $T\times N_o$, wide is $T_o\times N$, bal is $T_o\times N_o$, and miss is $(T-T_o)\times (N-N_o)$, With $N=T=200$, results are reported for four configurations of missingness. Evidently, the error depends on observability of $X_{it}$. The estimation error is largest if $X_{it}$ is in the miss block and smallest when $X_{it}$ is in the bal block. For a given block, the estimator errors are smaller when the factors are re-estimated from $\tilde X$. This is consistent with the theory. Results for $(N,T)=(300,500)$ and $(500,300)$ are similar.
In the second experiment, we generate data from a model with two factors to assess the adequacy of the asymptotic approximations. Two configurations of $(T,N) $ are considered: $(300, 500)$ and $(500, 300)$ with $T_o=.4T$ and $N_o=.6N$. In both cases, about 15% of the observations are missing. The factors and the loadings are drawn from the standard normal distribution once and held fixed. In each of the 5000 replications, a new batch of idiosyncratic errors are drawn from the normal distribution, so $X$ varies across replications but the locations of the missing values do not change. We then evaluate estimates of $C_{it}$ at four different $(i,t)$ chosen in the neighborhood of $N_o,T_o$ so that they come from four blocks as defined above.
Row 1 of Table (ref) reports the mean estimate of $C_{it}$ for (i) when all data are observed, (ii) tw0 (no re-estimation), (ii) tw+ (re-estimation), (iv) tp0, and (v) tp+. All five estimates are close to the true values, showing that tw0 and tp0 are consistent without iteration. The second row, which gives the standard deviation of the estimates in the monte-carlo, shows that the \textsc{tw+} and \textsc{tp+} estimates are slightly less variable than \textsc{tw} and \textsc{tp}, showing efficiency gains. The third row gives the mean of the estimated asymptotic standard errors which are quite similar to row two, showing that the asymptotic approximation is quite accurate. The next two rows labeled $q_{05}$ and $q_{95}$ present the 5 and 95 percentage points of the empirical distribution of the standardized estimates. They are quite close to the quantiles from the normal distribution of $\pm 1.64$. The last row shows that coverage is generally satisfactory and reinforces the adequacy of the asymptotic normal approximation. In summary, the proposed \textsc{tp} estimates are already consistent but the convergence rate of $C_{it}$ depends on the position of all $(i,t)$ as given in Proposition (ref). The updated estimates make use of additional information and have improved statistical properties.
Having shown that the properties of $\tilde X$ are satisfactory, we proceed to evaluate the covariance estimates using economically meaningful benchmarks. From a $T \times N^*$ panel of complete data $X^*$ treated as stock returns, we first calculate the “true" $N^*\times N^*$ covariance matrix $\bf \Sigma_{X^*}$. Using the $N^*\times 1$ vector of portfolio weights defined from $\bf \Sigma_{X^*}$, we compute $r^*_{pt}= \mathbf R^{*'}_t w^*_{p} $ and treat it as the “true" portfolio return at time $t$. We then randomly select $N$ stocks and assume missing values in the south-west block of the returns matrix $\mathbf R^*$. Hence the data matrix for analysis is of dimension $(T,N)$.
For each method, the (in-sample) bias and root-mean-squared error are computed for the following performance measures over $B=1000$ replications.
Of these measures, the first two are based on equally weighted portfolios while the last three are based on returns data. Portfolio volatility is often used as a risk measure when the benchmark is cash. Portfolio value-at-risk is used as a measure of downside risk.\footnote{Portfolio VaR is generally quoted as a positive number (ie, as a loss).} For call, we price plain vanilla European call option prices written on SP500 stocks calculated using Black-Scholes formula. The call price, variance, and covariance measures are calculated only for returns with incomplete data. Since we have $N_m$ series with missing values, we compute $N_m$ call prices and variances, and $ N_oN_m + (N_m)(N_m-1)/2$ covariances at each iteration.
The following notation is used Tables (ref) to (ref).
\paragraph{Strict Factor Model:} $X^*$ is simulated from a strict factor model:
where $\mathbb{I}_r$ is the $r \times r$ identity matrix, $\sigma_F^2$ and $\sigma_{\lambda}^2$ are the variances of each loading and factor, respectively, and $R_i^2$ is the percentage contribution of the systematic component to total variance (i.e. coefficient of determination) for series $i$. Note that setting the variance of $e_{it}$ in the way above guarantees that the coefficient of determination of each series is exactly equal to $R_i^2$. Since $r$ is assumed known, we set $r=5$ without loss of generality. We set $R^2_i$ to 0.6, $\sigma^2_{\lambda}$ to 1, and $\sigma^2_F$ to $0.035$, respectively. These give average volatility of $9.6\%$, similar to monthly SP 500 returns of $9.4\%$. The DGP abstracts from time varying loadings because the covariance estimators and will be affected by omitted time variation in the same way. Here, we focus on estimating the covariance of $\tilde X$.
Table (ref) reports results over 1000 replications for $(T,N)=(339, 100)$ with 15% of the values missing (60% of rows and 40% of columns). Whether we use tw or tp seems to make little difference for smj estimation. This can be due to the fact that either way, the missing idiosyncratic errors are set to zero and does not make a big difference for covariance estimation. Re-estimation is, however, beneficial for sm estimation especially when the first step is based on tw. This makes sense because tw uses a more restricted sample in the first step estimation and stands to gain more from re-estimation. Of the four methods, \textsc{sm+2} and \textsc{sm+4} which resample from the non-zero $\tilde e_{it}$ of series $i$ alone have smaller bias and root-mean-squared error.
Turning to the factor-based estimators, since the idiosyncratic errors are mutually uncorrelated by design, we had expected sfa and sfa+ which impose the strict factor structure to give better estimates than the sm counterparts. While sfa+ is comparable to sm2 and sm4, it is inferior to \textsc{sm+2} and \textsc{sm+4}. This can be due to the fact that our implementation of \textsc{sm} is an average of $S$ imputed covariances which has a noise reducing property. It is also possible that the \textsc{sm} estimators impose sparsity through univariate re-sampling, and as mentioned above, residual cross-section dependence can only come from the non-zero $\tilde e_{it}$ associated with the non-missing data. As these are estimates of errors that are uncorrelated by design, the \textsc{sm} estimators should return sample covariances that are approximately diagonal.
As robustness check, the top panel of Table (ref) presents results for 30% missing (80% of rows and 60% of columns). The relative performance of the estimation methods are similar across levels of missingness, and as expected, estmation accuracy of a given model is higher when the missingness is lower. The results so far are presented for $N<T$. Results for $T=200, N=250$ are shown in the bottom panel Table (ref). Note that with this design, each $\tilde{\mathbf \Sigma}_{\tilde X}(s,j)$ is verified to be singular. However, the averaged estimator $\bar {\mathbf \Sigma}_{\tilde X}(j)$ is not. Similar to results in Table (ref), sm+2 and sm+4 perform well and often better than sf+a.
\paragraph{Monte Carlo Calibrated to SP500 Returns} The final exercise is a Monte-Carlo that takes $X^*$ to be SP500 monthly returns between 1990 and 2018 to capture the approximate factor structure. The first factor in $X^*$ of dimension $348\times 339$ explains over 26.2% of the variations, the second explains 4.1%. The next three factors explain 3.8, 2.7, and 2.1% of the variation, respectively. As in the previous monte carlo, we random select $N=100$ stocks in each replication and set some values in the south-west block to missing.
The results based on 1000 replications are reported in Table (ref). Imputation leads to some reduction in bias and variance, but the errors are significantly smaller than no adjustment upon comparing sm1-4 with sm0. Though double imputation is generally better than single imputation, the precise gain depends on the risk measure and the estimation sample. As in the strict factor case, tw and tp give similar risk-performance errors. Also as in the strict factor case, sfa+ is comparable to sm+1 and \textsc{sm+3}, but has larger errors than \textsc{sm+2} and \text{sm+4}. However, improvements from re-estimation both in terms of bias and RMSE are larger for this DGP that mimics the approximate factor structure.
Overall, we find that single imputation is strongly preferred over no imputation but is inferior to double imputation. Ultimately, the appropriate method depends on the data generating process and we offer several alternatives worthy of consideration. Averaging the covariances of imputed data overlaid with resampled errors has two desirable effects that we had not anticipated: it reduces noise, and the averaged estimator is full rank. This interesting finding and the role of resampling in covariance estimation is an interesting problem that warrants further investigation.
This paper provides three sets of results. The first is a tp algorithm that can consistently estimate the entire low rank matrix without iteration and a distribution theory for the estimates is provided. The second result makes precise the conditions under which factors estimated from incomplete data can be treated as known in factor-augmented regressions. The third pertains to estimation of covariances from incomplete data. We consider different schemes to compensate for an omitted error in the level estimates. Implications for using imputed factors in augmented regressions are also discussed.
\baselineskip=12.0pt
\baselineskip=15.0pt