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.
102,815 characters · 22 sections · 108 citation commands
The Factor-Lasso and K-Step Bootstrap Approach for Inference in High-Dimensional Economic Applications
\address{Booth School of Business, University of Chicago, Chicago, IL 60637} \email{[email removed]}
\address{Department of Economics, Rutgers University, New Brunswick, NJ 08901} \email{[email removed]}
\allowdisplaybreaks{
JEL classification: C33, C38.
Keywords: panel data, treatment effects, high-dimensional
Data in which there are many observable variables available for each observation, i.e. “high-dimensional data,” are increasingly common and available for use in empirical applications. Having rich high-dimensional data offers many opportunities for empirical researchers but also poses statistical challenges in that regularization or dimension reduction will generally be needed for informative data analysis. The success of regularized estimation for either forecasting or inference using high-dimensional data relies on using a regularization device that is appropriate for the type of data at hand. Effective regularization imposes substantive restrictions in estimation, and the resulting estimates can perform very poorly, for example suffering from large biases and missing important explanatory power, when the restrictions provide poor approximations to the underlying data generating mechanism. It is thus important to employ regularized estimators that accommodate sensible beliefs about the structure of an underlying econometric model.
Two structures which are common in the econometrics literature are sparse structures and factor structures. To fix ideas, consider the linear regression model
where $i \leq n$ indexes individual observations, $y_i$ is the observed outcome of interest, $x_i$ is a $p \times 1$ vector of observed predictor variables with $p \gg n$ allowed, and $\varepsilon_i$ is a regression disturbance. A sparse structure essentially imposes that the number of non-zero elements in $\theta$ is small. Intuitively, the sparse structure relies on the belief that the majority of the explanatory power in the observed predictor variables concentrates within a small number of the available variables. Estimators that are appropriate for sparse models, such as the lasso or variable selection procedures, may perform very poorly when the true model is “dense” in the sense that there are many non-zero elements in $\beta$ that are moderate in magnitude.
A commonly employed version of a linear factor model employs a different structure where
Here $f_i$ denotes a latent $K\times 1$ vector of factors with $K \ll n$ that are important in determining both the observed outcome of interest, $y_i$, and the observed $p \times 1$, with $p \gg n$, vector of observed predictor variables $x_i$. Within this structure, one may obtain estimates of the latent factors and build a model for the outcome given the extracted factors; see, e.g. bai03, BN02, SW02 and fan2016sufficient. The basic factor model differs markedly from the sparse linear model ((ref)). Importantly, data generated from model ((ref))-((ref)) would generally result in a dense coefficient vector $\theta$ in the regression of $y_i$ onto $x_i$, and sparsity based estimation strategies would tend to perform poorly as a result. Of course, if the data generated by the sparse model ((ref)), common factors will generally not capture the explanatory power, which loads on a small number of the raw regressors, and pure factor-based estimation will perform poorly.
In this paper, we propose a simple model that nests both the sparsity-based and factor-based structures. The model allows for the observed predictors to have a factor structure but then allows both the factors and the factor residuals, the $U_i$ in equation ((ref)), to load in the outcome equation. That is, we replace ((ref)) with
and impose that $\theta$ is sparse. This model allows for the fact that all of the relevant explanatory power in the predictors may not be captured entirely by the factors but imposes that any predictive power not captured by the factors concentrates on only a few elements of the high-dimensional covariate vector. ((ref)) clearly reduces to ((ref)) when there is no factor structure in $x$ and reduces to ((ref)) when $\theta = 0$. We note that this model shares much in common with factor augmented regression models, e.g. BN06 and BBE:FAVAR, with the key points of departure being that we do not assume the identity of the additional variables to include in the model is known and that $U$ is not observable. HMC:PFM consider a model that shares the essential structure of ((ref)) and ((ref)) from a Bayesian standpoint. They show that forecasts obtained from their Bayesian estimator of this model tend to outperform forecasts obtained based on either pure sparsity or pure factor based models.
The first key contribution of the present paper is offering a practical estimation and inference procedure that is appropriate for inference in a panel generalization of the model given by equations ((ref)) and ((ref)) and providing a formal treatment of the procedure's theoretical properties. Specifically, we proceed by first running a factor extraction step and taking residuals from regressing each observed variable on the estimated factors. Using these residuals, we then follow the lasso-based estimation and inference procedures of BCHK:FE. We show that the resulting estimator of parameters of interest specified ex ante by the researcher is asymptotically normal with readily estimated asymptotic variance under sensible conditions. These conditions allow for errors in selection of the elements of the covariate vector that load after controlling for the factors but maintain sufficiently strong conditions to allow oracle selection of the number of factors. The theoretical analysis is substantially complicated by the fact that factors and factor-residuals are not observed and must be extracted from the data. The estimation error in this extraction then enters the second step nonlinear and non-smooth lasso problem. Due to this complication, the theoretical results in this paper make use of arguments that, to our knowledge, are not implied by results existing in the current factor modeling literature or the current lasso literature. These results may be of some interest outside of the context of establishing the properties of our proposed inferential procedure.
By addressing estimation and inference in an interesting high-dimensional factor augmented regression model appropriate for panel data, our paper contributes to the rapidly growing literature dealing with obtaining valid inferential statements following regularized estimation. See, for example, belloni2012sparse,BCK-LAD,BCY-honest,BCFH:Policy,belloni2014inference,BCHK:FE, BerkEtAl13, DoubleML, dezeure2016high, fan2001variable, FL11, Farrell:JMP, GautierTsybakovHDIV, AdaptiveGraphLasso, FST14, JM:ConfidenceIntervals, Kozbur:FS, LeeTaylorScreening, LeeEtAlLasso, LockhartEtAlLasso, LoftusStepwise, TaylorEtAlAdaptive, vdGBRD:AsymptoticConfidenceSets, WagerAthey:RandForTE, and ZhangZhang:CI for approaches to obtaining valid inferential statements in a variety of different high-dimensional settings.
As a second main contribution, we offer a new, computationally convenient bootstrap method for inference. Specifically, we consider a bootstrap where we apply our main procedure, including extraction of factors and lasso estimation steps, within each bootstrap replication. As computation of the lasso estimator within each bootstrap sample may be demanding, we explicitly consider a k-step bootstrap following andrews2002higher where we start at the lasso solution from the full sample and then iterate a numeric solution algorithm for the lasso estimator for k-steps. We make use of solution algorithms for which the updates are available in closed form which leads to fast computation. We provide high-level conditions under which the procedure provides asymptotically valid inference for parameters of interest and provide specific examples with lower level conditions. The k-step bootstrap we propose complements other bootstrap procedures that have been proposed for lasso-based inference, for example, BCFH:Policy, chatterjee2011bootstrapping, CCK:AOS13, and dezeure2016high. In particular, the approach we take is something of a middle ground between CCK:AOS13, which uses resampling of model scores to avoid recomputation of the lasso estimator, and dezeure2016high which fully recompute the lasso solution within each bootstrap replication. The former approach is extremely computationally convenient and asymptotically valid but does not capture any finite sample uncertainty introduced in the lasso selection, while the latter may be computationally cumbersome due to fully recomputing the lasso solution within each iteration. We note that the bootstrap procedure could be easily applied outside of the specific model considered in this paper and that the technical analysis here is new and may be of interest outside of the present context.
The remainder of this paper is organized as follows. In Section (ref), we describe the panel factor-lasso model and outline the basic algorithm we will employ for inference. We present formal results for the proposed procedure in Section (ref), providing regularity conditions under which the estimator of parameters of interest is asymptotically normal and valid confidence statements may be obtained. Section (ref) describes the k-step bootstrap approach in detail and provides a formal analysis establishing the validity of the resulting bootstrap inference. Section (ref) discusses the factor extraction part of the problem in more detail and provides examples with accompanying low-level conditions that are sufficient for the high-level conditions stated in Section (ref). We then provide simulation and empirical examples that motivate the model we consider and illustrate the use of the estimation procedure in Section (ref). Key proofs are collected in an appendix with additional results provided in a supplementary appendix.
Throughout the paper, we use $\|\beta\|_1$ and $\|\beta\|_2$ to respectively denote the $\ell_1$- and $\ell_2$- norms of a vector $\beta$; use $\|A\|$ and $\|A\|_F$ to respectively denote the spectral and Frobenius norms of a matrix $A$. In addition, denote by $|J|_0$ as the cardinality of a finite set $J$. Finally, for two positive sequences $a_n, b_n$, we write $a_n\asymp b_n$ if $a_n=O(b_n) $ and $b_n=O(a_n).$
Consider the linear panel model defined by
where $i \leq n$ indexes cross-sectional observations, $t \leq T$ indexes time series observations, $X_{it}$ are observed potentially confounding variables, and $d_{it}$ is an a priori specified “treatment” variable of interest.\footnote{Our results will immediately apply to the case where $d_{it}$ is an $r \times 1$ vector with $r$ fixed. The analysis could also be extended to handle unbalanced panels where observations are missing at random. We omit both cases for convenience.} $f_i$ is a $K \times 1$ vector of latent factors with time-varying $K \times 1$ factor loading vectors $\xi_t$, $\delta_{dt}$ and $p \times K$ dimensional factor-loading matrix $\Lambda_t$. We will take asymptotics where $\dim(X_{it})=p \rightarrow \infty$, $n \rightarrow \infty$, and $T$ is either fixed or growing slowly relative to $n$ and $p$ when stating our formal results, and we explicitly allow for scenarios where $p \gg nT$. $K$ is assumed fixed throughout the paper. Our object of interest is the parameter $\alpha$ on the variable of interest $d_{it}$. Following HMC:PFM, we refer to the model ((ref))-((ref)) as the “panel partial factor model” (PPFM).\footnote{HMC:PFM consider a similar structure to ((ref))-((ref)) which excludes the individual and time effects and imposes that the $\epsilon_{it}$ are i.i.d. Gaussian innovations. They refer to this model as a partial factor model.}
In each equation, we also allow for additive unobserved individual effects, $(g_i,\zeta_i,w_i')$, and time specific effects, $(\nu_t,\mu_t,\rho_t')$, where $g_i$, $\zeta_i$, $\nu_t$, and $\mu_t$ are scalars and $w_i$ and $\rho_t$ are $p \times 1$ vectors. We do not impose structure over the individual or time specific effects and thus treat them as fixed effects. This treatment differentiates the common factors, $f_i$, from the additive heterogeneity $(g_i,\zeta_i,w_i')$ and $(\nu_t,\mu_t,\rho_t')$ as we impose that the $f_i$ are common to each observed series with common, time-varying loadings. Term $U_{it}$ represents the part of the observed $X_{it}$ that is orthogonal to the factors and unobserved time and individual specific heterogeneity. We allow $U_{it}$ to be correlated to both the outcome and variable of interest after controlling for the factors and individual and time fixed effects. Because $p \gg nT$, we assume that $\theta$ and $\gamma_d$ are approximately sparse vectors. We assume that observed right-hand side variables are strictly exogenous so that ${\mathrm{E}}[\eta_{it}|X_{i1},...,X_{iT}] = 0$ and ${\mathrm{E}}[\epsilon_{it}|X_{i1},...,X_{iT},d_{i1},...,d_{iT}] = 0$. We will assume that data are iid across $i$ but allow for dependence across time periods, $t$. Finally, we note that while we treat the PPFM defined in ((ref))-((ref)) in the formal analysis, the results clearly apply to models without additive fixed effects or to a single cross-section.\footnote{We consider a cross-sectional instrumental variables version of the model in both a simulation and an empirical example.}
As noted in the Introduction, the PPFM generalizes the high-dimensional sparse fixed effects model examined in BCHK:FE and conventional large-dimensional factor models and factor augmented regression models; e.g. BN06. The PPFM is also related to, but distinct from, interactive fixed effects models as in, for example, bai,bai2014theory, MW10,MW11, pesaran and Su13.\footnote{See also bonhomme:manresa for a distinct but related approach based on a grouped fixed effects model.} A simple version of the interactive fixed effects model analogous to ((ref)) is
In this model, $z_{it}$ represents a known, low-dimensional set of variables that must be controlled for in addition to the factors in $f_i$. There appear to be three key distinctions between the high-dimensional PPFM and interactive fixed effects approaches. First, we relax the assumption that one knows the exact identity of the variables that should appear in the model, $z_{it}$, by allowing for a high-dimensional set of observed potential confounds in $X_{it}$. Second, we allow for the fact that the relevant explanatory power in the predictors may not be captured entirely by the factors, but impose that any predictive power not captured by the factors concentrates on only a few elements of the high-dimensional vector $U$. Third, we directly extract estimates of the factors and $U$ from $X$ which can proceed even when $T$ is small. Approaches to estimating the interactive fixed effects structure rely on having a large number of observations in both the time series and cross-sectional dimensions. We thus view the PPFM and interactive fixed effects approaches as complementary where one may prefer one or the other depending on the nature of the data at hand.
To estimate $\alpha$, we begin by taking the within transformation of all observed variables to remove the fixed effects. To this end, let $$ \tilde z_{it}=z_{it}-\bar z_{\cdot t}- \bar z_{i \cdot} +\bar{\bar z} $$ for any variable $z_{it}$ where $\bar z_{\cdot t} = \frac{1}{n}\sum_{i=1}^{n} z_{it}$, $\bar z_{i \cdot} = \frac{1}{T}\sum_{t=1}^{T} z_{it}$, and $\bar{\bar z} = \frac{1}{nT}\sum_{i = 1, t = 1}^{n,T} z_{it}$. We can then define a demeaned model as
After removing the additive unobserved heterogeneity, we estimate the (demeaned) latent factors as well as the (demeaned) idiosyncratic components from the model $\tilde X_{it}=\tilde\Lambda_t'\tilde f_i+\tilde U_{it}$.\footnote{We note that recovering the untransformed $f_i$ and $U_{it}$ would only be possible with large $n$ and $T$ due to the presence of the unrestricted fixed effects. Fortunately, recovering these quantities is unnecessary within the model with common coefficients $\theta$, $\gamma_d$, and $\alpha$ as only $\tilde f_i$ and $\tilde U_{it}$ appear in the equations of interest. This simplification would not generally occur if we allowed heterogeneity in $\theta$, $\gamma_d$, or $\alpha$ over time or across individuals, and we would need to consider incidental parameters bias introduced by removing the additive fixed effects. We leave exploration of this issue to future research.} Let $\widehat F=(\widehat f_1,...,\widehat f_n)'$ be the $n\times K$ matrix of estimated factors. We shall discuss some examples of $\widehat F$ in Section (ref). Given $\widehat F$, we estimate $\tilde\Lambda_t$ and $\tilde U_{it}$ by least squares:
Substituting ((ref)) to ((ref)), we obtain
Now let $\tilde Y_t = (\tilde y_{1t},...,\tilde y_{nt})'$ and $\tilde D_t = (\tilde d_{1t},...,\tilde d_{nt})'$ denote the vectors of outcome and treatment variable within each time period $t$. We next regress $\tilde Y_t$ and $\tilde D_t$ onto the extracted factors $\widehat F$ time period by time period to obtain $\{\widehat \delta_{yt}\}_{t=1}^{T}$ and $\{\widehat \delta_{dt}\}_{t=1}^{T}$ for
We then run the lasso with the residuals from each of these factor regressions as dependent variable and the estimated factor disturbances $\widehat U_{it}$ as predictors. That is, we obtain
where the tuning parameter $\kappa_n$ is chosen as, for some $c_0>1$ and $q_n \rightarrow 0$,
and $\widehat\Psi^y$ and $\widehat\Psi^d$ are diagonal penalty loading matrices. Given the fixed effects panel structure, we use the clustered penalty loadings of BCHK:FE which have diagonal elements defined as
where $\widehat e_{it}$ is an estimator of $\tilde e_{it} = \tilde y_{it} - \tilde \delta_{yt}'\tilde f_i- \tilde U_{it}'\gamma_y$ and $\widehat \eta_{it}$ is an estimator of $\tilde\eta_{it} = \tilde d_{it} - \tilde \delta_{dt}'\tilde f_i - \tilde U_{it}'\gamma_d$. \footnote{We obtain $\widehat e_{it}$ and $\widehat \eta_{it}$ through an iterative algorithm similar to that of belloni2014inference, which starts from a preliminary estimate. In addition, we use $c_0 = 1.1$ and $q_n = .1/\log(n)$ in the simulation and empirical examples. }
For the final step, we adopt the post-double-selection procedure of belloni2014inference. Let $\widehat J = \{j \leq p: \tilde \gamma_{y,j} \neq 0 \} \cup \{j \leq p: \tilde \gamma_{d,j} \neq 0\}$, and let $\widehat U_{it, \widehat J}$ be a subvector of $\widehat U_{it}$ whose elements are $\{\widehat U_{it, j}: j \in \widehat J\}$. We then run the regression of $\tilde y_{it}-\widehat\delta_{yt}'\widehat f_i$ on $\widehat U_{it, \widehat J}$ and $\tilde d_{it} - \widehat\delta_{dt}'\widehat f_i$ on $\widehat U_{it, \widehat J}$ and obtain
The final estimator of $\alpha$ is then given by
where $\widehat e_{it} = \tilde y_{it}- \widehat\delta_{yt}'\widehat f_i-\widehat U_{it,\widehat J}'\widehat\gamma_y$ and $\widehat \eta_{it} = \tilde d_{it}- \widehat\delta_{dt}'\widehat f_i-\widehat U_{it,\widehat J}'\widehat\gamma_d$ are the residuals from the regressions specified in ((ref)) and ((ref)).
The estimator $\widehat\alpha$ can be expressed more compactly in matrix form. Write
In addition, for a matrix $A$, define $M_A=I- A(A'A)^{-}A'$, where $(A'A)^{-}$ represents a generalized inverse of $A'A$. Then it is straightforward to verify that $$\widehat e = M_{\widehat U_{\widehat J}}(I_T\otimes M_{\widehat F})\tilde Y, \ \text{and} \ \widehat\eta = M_{\widehat U_{\widehat J}}(I_T\otimes M_{\widehat F})\tilde D $$ are the estimated residuals $(\tilde e_1',...,\tilde e_T')'$ and $(\tilde \eta_1',...,\tilde \eta_T')'$ defined above. Then $$ \widehat\alpha=(\widehat\eta'\widehat\eta)^{-1}\widehat\eta'\widehat e. $$
Note that the estimator $\widehat\alpha$ is numerically equivalent to the coefficient on $\tilde d_{it}$ in the regression of $\tilde y_{it}$ on $\tilde d_{it}$, $\widehat f_{i}$ interacted with time dummy variables, and $\widehat U_{it, \widehat J}$. In Theorem (ref) of the next section, we verify that inference for $\widehat\alpha$ can proceed using the output from this OLS regression as long as clustered standard errors (e.g. arellano:feinf, bdm:cluster, and hansen:cluster) are used.
The following algorithm summarizes the estimation strategy detailed above.
In this section, we present a set of sufficient conditions under which we establish asymptotic normality of $\widehat\alpha$ and provide a consistent estimator of its asymptotic variance. Throughout we consider sequences of data generating processes (DGPs) where $p$ increases as $n$ and $T$ increase and where model parameters are allowed to depend on $n$ and $T$. We suppress this dependence for notational simplicity. We use the term “absolute constants” to mean given constants that do not depend on the DGP.
Write $\epsilon_t=(\epsilon_{1t},...,\epsilon_{nt})'$, $\eta_t=(\eta_{1t},...,\eta_{nt})'$, and $U_{t}=(U_{1t}',...,U_{nt}')'.$ Similarly, let $\epsilon_i=(\epsilon_{i1},...,\epsilon_{iT})'$ and $\eta_i=(\eta_{i1},...,\eta_{iT})'$, $U_i=(U_{i1}',...,U_{iT}')'$.
Our first two conditions collect various restrictions on dependence, tail behavior, and moments of the unobserved features of the model. We assume there are positive absolute constants $C_1,C_2$ and $C_3$ such that the following assumption holds.
Assumption (ref) collects reasonably standard regularity conditions that restrict the dependence across observations and tail behavior of random variables. These conditions impose that the unobserved variables in the model are cross-sectionally independent, are weakly dependent and stationary in the time series, and have sub-Gaussian tails. Assumption (ref).(iii) further imposes weak conditional dependence in the factor residuals, $U_{it}$. In the simple case where $U_{it}$ is independent of $f_i$, $\eta_i$, and $\epsilon_i$ for all $t$, this condition reduces to weak intertemporal correlation and no strong dependence among the columns of $U_{it}$. Importantly, it does not imply that all correlation among the observed $X_{it}$ is captured by factors but allows for the presence of a rich covariance structure in the part of $X_{it}$ that is not linearly explained by the factors. The condition also allows for some dependence between “control” variables $U_{it}$ and structural unobservables $\eta_i$ and $\epsilon_i$ but restricts the magnitude of any such dependence so that it is asymptotically negligible. Finally, condition (v) requires that all the low dimensional parameters are well bounded.
Recall that $e_{it}=\alpha\eta_{it}+\epsilon_{it}$.
Assumption (ref) collects additional high-level moment bounds. The bounds on moments of normalized sums in Condition (i) could be established under a variety of sufficient lower level conditions. Condition (ii) places restrictions on the dependence between $\{U_{it},e_{it}\}_{i=1,t=1}^{n,T}$ and $\{f_i\}_{i=1}^{n}$.
Before stating the next assumption, we decompose the high dimensional coefficients as $$ \gamma_y=\underbrace{\gamma_y^0}_{\text{exactly sparse}}+\underbrace{R_y}_{\text{remainder}} \quad \text{and} \quad \gamma_d=\underbrace{\gamma_d^0}_{\text{exactly sparse}}+\underbrace{R_d}_{\text{remainder}} $$ where $\gamma_y^0$ and $\gamma_d^0$ are sparse vectors that approximate the potentially dense true coefficient vectors $\gamma_y$ and $\gamma_d$ and $R_y$ and $R_d$ represent approximation errors. Let $ J=\{j\leq p: \gamma_{y,j}^0\neq 0\}\cup\{j\leq p: \gamma_{d,j}^0\neq 0\} $ be the union of the support of the exactly sparse components.
Assumption (ref) collects restrictions on the quality of the approximation provided by $\gamma_y^0$ and $\gamma_d^0$ and rates of growth of model complexity as measured by $J$ and $p$ and sample sizes in the cross-sectional and time series dimension. Condition (iii) imposes the somewhat nonstandard requirement that $T$ be much smaller than $n$. The need for this condition arises from the fact that we need to obtain high-quality estimates of the idiosyncratic term in the factor equation, $U_{it}$, which depends on accurately estimating both the unknown factors and the loadings. Estimating the loading matrix $\Lambda_t$ well for any given $t$ requires a relatively large $n$, and we thus require $T$ to be smaller than $n$ as the number of unknown loading matrices $\{\Lambda_t\}_{t\leq T}$ is $O(T)$.
Our next assumption restricts the covariance matrix of the within-transformed factor residuals $\tilde U_{it}$.
{Maintaining Assumptions (ref)-(ref),} a simple sufficient condition for Assumption (ref) is that all the eigenvalues of $\frac{1}{nT}\sum_{i}\sum_tE( U_{it}-\bar U_{i,\cdot})(U_{it}-\bar U_{i\cdot})'$ are well bounded. This is a typical condition in high-dimensional approximate factor models (e.g., bai03, SW02). It ensures that the idiosyncratic components are weakly dependent and therefore the decomposition $\tilde X_{it}=\tilde \Lambda_t'\tilde f_i+\tilde U_{it}$ is asymptotically identified (as $p\to\infty$).
Finally, we present high-level conditions on the accuracy of $\widehat F$ in Assumption (ref). The high-level conditions potentially allow for many estimators of the factors, and we verify that these conditions hold under more primitive assumptions for the case of estimating the factors using PCA in Appendix (ref).
One of the major technical tasks of this paper is to show that the effects of estimating the latent factor and idiosyncratic terms are stochastically dominated by the plug-in tuning parameter $\kappa_n$ in ((ref)). Since $\kappa_n\asymp\sqrt{\frac{\log p}{nT}}$, this is a strong requirement, and gives rise to Assumption (ref) (and Assumption (ref) below for the bootstrap sample). Technically, existing results in the literature on estimating factors models are not directly applicable to verify these conditions. In Appendix (ref), we show that
when $\widehat f_t$ is estimated via PCA. While this result is essentially standard and allows conditions involving $\Delta_F$ to be directly verified, however, it does not imply the uniform convergence condition $\max_{t\leq T}\|\widehat f_t-H'\tilde f_t\|_2$. Nor is this result sufficient to verify the other stated conditions because other terms, e.g. $\Delta_{eg}, \Delta_{fum}, \Delta_{fe}$, involve “weighted averages" of $\{\widehat f_i- H'\tilde f_i\}$ whose rates of convergence can be derived and shown to be faster than that of $\Delta_F=\frac{1}{pT}+\frac{1}{n^2}+\frac{1}{nT^2}$. For instance, if we use a simple Cauchy-Schwarz inequality to bound $\Delta_{ud}$, we would have $$ \max_{m, k\leq p}\|\frac{1}{nT} \sum_{i=1}^n \sum_{t=1}^T(\widehat f_i -H'\tilde f_i) \tilde U_{it,m} \tilde h_{tk}' \|_F^2\leq \frac{1}{n}\sum_{i=1}^n\|\widehat f_i-H'\tilde f_i\|_2^2 \max_{m,k\leq p}\frac{1}{n}\sum_{i=1}^n\|\frac{1}{T}\sum_{t=1}^T \tilde U_{it,m} \tilde h_{tk}\|_2^2. $$ It can be shown that $\max_{m,k\leq p}\frac{1}{n}\sum_{i=1}^n\|\frac{1}{T}\sum_{t=1}^T \tilde U_{it,m} \tilde h_{tk}\|_2^2=O_P(\frac{\log p}{T})$, so this crude bound gives us $\Delta_{ud}=\Delta_F\sqrt{\frac{\log p}{T}}$. Unfortunately, this bound is not sharp enough to verify the condition $ \Delta_{ud}=o(\sqrt{\frac{\log p}{nT}})$ unless $n=o(pT)$. In the special case that $T$ is fixed, requiring $n=o(p)$ is a restrictive condition. Rather than relying on these crude bounds, we achieve sharper bounds by directly deriving the rate of convergence for each required term in Appendix (ref) which relies on some novel technical work. These conditions only require $n=o(p^2T)$ which provides much more freedom on the ratio $n/p$.
The asymptotic variance of $\widehat\alpha$ will depend on the quantities $$ \sigma_{\eta\epsilon}=\operatorname{Var} \left(\frac{1}{\sqrt{nT}}\sum_{i=1}^n\sum_{t=1}^T( \eta_{it}-\bar \eta_{i\cdot})(\epsilon_{it}-\bar\epsilon_{i\cdot}) \right) \quad \text{and} \quad \sigma_{\eta}^2=\frac{1}{nT}\sum_{i=1}^n\sum_{t=1}^T\operatorname{Var}(\eta_{it}-\bar \eta_{i\cdot}) $$ for which $$ \widehat \sigma_{\eta\epsilon}=\frac{1}{nT}\sum_{i=1}^n \left(\sum_{t=1}^T\widehat\eta_{it}\widehat\epsilon_{it}\right)^2 \quad \text{and} \quad \widehat \sigma_{\eta}^2=\frac{1}{nT}\sum_{i=1}^n\sum_{t=1}^T\widehat\eta_{it}^2 $$ are natural estimators. Note that $\widehat\sigma_{\eta\epsilon}$ is just the usual clustered covariance estimator with clustering at the individual level.
The main implication of Theorem (ref) and Corollary (ref) is that $\widehat\alpha$ converges at a $\sqrt{nT}$ rate and that inference may proceed using standard asymptotic confidence intervals and hypothesis tests. Importantly, the inferential results hold uniformly across a large class of approximately sparse models which include cases where perfect selection over which elements of $\tilde U_{it}$ enter the model is impossible even in the limit. It is also important to highlight that the conditions on estimation of the factors do rule out the presence of weak factors, and the inferential results do not hold uniformly over sequences of models in which perfect selection of the number of factors and fast convergence of the factors and factor loadings do not hold. The difficulty with handling weak factors arises due to the entry of the estimation errors of the factors in the cluster-lasso problems ((ref)) and ((ref)) and the non-smooth and highly nonlinear nature of this problem. Extending the results to accommodate the presence of weak factors and imperfect selection of the number of factors would be an interesting direction for further research.
In this section, we present a computationally tractable bootstrap procedure that can be used in lieu of the plug-in asymptotic inference formally presented in Theorem (ref) and Corollary (ref). While well-developed in low-dimensional settings, there are relatively few formal treatments of bootstrap procedures in high-dimensional settings, though see chatterjee2011bootstrapping, CCK:AOS13, BCFH:Policy, and dezeure2016high for important existing treatments. In the following, we consider a bootstrap procedure which only approximately solves the cluster-lasso problem within each bootstrap replication and thus may remain computationally convenient while also intuitively capturing the sampling variation introduced in the lasso selection.
Let $D^*=\{\tilde y_{it}^*, \tilde d_{it}^*, \tilde X_{it}^*\}_{i\leq n, t\leq T}$ denote a sample of bootstrap data, and let $\widehat\alpha^*$ be the estimator obtained by applying the factor-lasso estimator with data $D^*$. Let $B$ denote the number of bootstrap repetitions.
A potential computational problem with bootstrap procedures for lasso estimation is that one needs to solve $B$ lasso problems where $B$ will typically be fairly large. To circumvent this problem, we adopt the approach of andrews2002higher by using the fact that the complete lasso estimator based on the original data, denoted by $\tilde \gamma_{lasso}$, should be close to the complete lasso estimator based on bootstrapped data $D^*$, denoted by $\tilde\gamma^*_{lasso}$. Hence, within each bootstrap replication, we can use $\tilde \gamma_{lasso}$ as the initial value for solving the lasso problem and iteratively update the estimator for $k$ steps. Denote the resulting $k$-step bootstrap lasso estimator by $\tilde\gamma^*$. We simply use $\tilde\gamma^*$ in place of $\tilde\gamma^*_{lasso}$ wherever the solution to a lasso problem shows up in the factor-lasso problem. The main result of this section is showing that the $k$-step bootstrap procedure is first-order valid for statistical inference about $\alpha$ as long as the minimization error after $k$ steps is less than the statistical error (i.e. $o_{P^*}(({nT})^{-1/2})$.
The substantive difference between the present context and andrews2002higher is that andrews2002higher makes use of Newton-Raphson updates for the k-steps while face a regularized optimization problem at each iteration. Tractability relies on the fact that there are a variety of procedures for updating within the lasso problem that are available in closed form. Using these analytic updates greatly reduces the overall computational task and makes a k-step bootstrap procedure attractive within the lasso context.
Specifically, consider the following lasso problems on the bootstrap data. Let
where
The definitions of $\{\tilde y_{it}^*, \tilde d_{it}^*, \widehat\delta_{yt}^{*}, \widehat\delta_{dt}^{*}, \widehat f_i^*, \widehat U_{it}^{*} \}_{i\leq n, t\leq T}$ will be formally given below. Let $\tilde\gamma_{y}$ and $\tilde \gamma_d$ be the lasso solutions obtained from the original data. Also, note that we fix the value of $\kappa_n$ and of the penalty loadings $\widehat\Psi^y$ and $\widehat\Psi^d$ to the same values as used to obtain the solutions $\tilde\gamma_y$ and $\tilde\gamma_d$ in the original data.
Within each bootstrap replication, we then approximately solve the lasso problems ((ref)) by applying the following procedure. The maximum number of steps $k$ to be taken should be determined on a case-by-case basis according to the available computational capacity.\footnote{In applications where obtaining the full lasso solution is not too burdensome, one may simply iterate to convergence.}
There are a variety of iteration mappings that can be used in Step (A2) of the k-step lasso problem. A commonly used and simple mapping is the “coordinate descent method,” also known as the “shooting method,” studied by fu1998penalized.\footnote{Another commonly used iterative scheme that could readily be applied in the present setting is the “composite gradient method” (e.g. nesterov2007gradient and agarwal2012fast). We choose to focus on the coordinate descent method as our concrete example as it does not rely on additional tuning parameters and performed well numerically in preliminary simulation experiments. In addition, coordinate descent requires weaker regularity conditions than the composite gradient method for our theoretical analysis.} For solving problem ((ref)), write the solution after the $l^{\text{th}}$ iteration as $\gamma_{y,l}=(\gamma_{y,l,1},...,\gamma_{y,l,p})'$. The coordinate descent method updates $\gamma_{y,l+1}$ by iteratively cycling through all coordinates. Specifically, we solve the following one-dimensional optimization problem for $m=1,...,p$,
Here $m^-=\{j: j< m\}$; and $\gamma_{y,l+1,m^-}$ and $\widehat U_{it, m^-}^{*} $ are $\mathbb{R}^{m-1}$ dimensional vectors whose components are respectively those of $\{\gamma_{y,l+1, j}: j<m\}$ and $\{\widehat U^*_{it, j}: j<m\}$. Similarly, $m^+=\{j: j>m\}$; and $\gamma_{y,l,m^+}$ and $\widehat U_{it,m^+}^{*} $ are $\mathbb{R}^{p-m}$ dimensional vectors whose components are respectively those of $\{\gamma_{y,l, j}: j>m\}$ and $\{\widehat U^*_{it, j}: j>m\}$. When $m=1$, $m^-$ is empty; and when $m=p$, $m^+$ is empty. In these cases, the corresponding subvectors, $\gamma_{y,l+1,m^-}$ and $\widehat U_{it, m^-}^{*} $ or $\gamma_{y,l,m^+}$, and $\widehat U_{it,m^+}^{*} $, are defined as zero. Note that when $\gamma_{y, l+1,m}$ is being updated the previous $m-1$ elements have already been updated, while the remaining $p-m$ elements are yet to be updated. Thus, $\gamma_{y, l+1,m^-}$ is a subvector of $\gamma_{y, l+1}$, but $\gamma_{y, l,m^+}$ is a subvector of $\gamma_{y, l}$. Denote by $\gamma_{y, l+1}^{(m)}:=(\gamma_{y, l+1, m^-}, \gamma_{y, l+1,m}, \gamma_{y, l, m^+})'$ the vector that results immediately after the $m^{\text{th}}$ coordinate has been updated during the $(l+1)^{\text{th}}$ iteration. When $m=p$, all the components have been updated; and we obtain $\gamma_{y, l+1}:=\gamma_{y, l+1}^{(p)}.$
Importantly, ((ref)) is a one-dimensional $\ell_1$-penalized quadratic problem which has an analytical solution given by the soft thresholding operation:
where $ Z_{it, l,m}^*:=\tilde y_{it}^*-\widehat\delta_{yt}^{*'}\widehat f_i^*-\widehat U_{it,m^-}^{*'}\gamma_{y,l+1, m^-}-\widehat U_{it,m^+}^{*'}\gamma_{y,l, m^+}, $ $(x)_+=\max\{x,0\}$, and $sgn(x)$ takes the sign of $x$. Therefore, the mappings in ((ref)) are given by $$ \mathcal{S}_y(\gamma_{y,l})=(\gamma_{y, l+1, 1},...,\gamma_{y, l+1, p})', \quad \text{where each $\gamma_{y, l+1, m}$ is given in (\ref{eq4.3}).} $$ $\mathcal{S}_d(\gamma_{d,l})$ is obviously defined similarly.
With the k-step lasso program defined, we now state the complete algorithm for the proposed k-step bootstrap procedure. We make use of a wild residual bootstrap to generate the data at each bootstrap replication.
In the following, we present conditions under which we verify that the bootstrap confidence intervals are asymptotically valid: $$ P\left(\alpha\in \left[\widehat\alpha\pm \frac{q_{\tau}^*}{\sqrt{nT}}\right]\right)\to 1-\tau. $$
The first assumption imposes high-level conditions that will admit the use of general updating rules in ((ref)) of the $k$-Step Lasso Iteration. The assumption provides high-level conditions on the computational properties and sparsity of the solution resulting after taking $k$ iterations in the solution of the lasso problem. Recall that $ \tilde\gamma_y^*=\gamma_{y,k}$ and $\tilde\gamma_d^*=\gamma_{d,k}.$
Condition (i) requires that the minimization error should be negligible compared to the statistical error after $k$ iteration steps. Condition (ii) guarantees the sparsity of the iterated solutions. As a concrete example, we verify both conditions for the coordinate descent method. We note that, to the best of our knowledge, showing the $|J|_0$-sparsity of the $k$-step iterated coordinate descent estimator has not been done previously when $p$ is potentially much larger than $n$ and may be of some independent interest.
We next impose a fairly standard notion of regularity on the high-dimensional component $\tilde U_{it}$ which shows up in the infeasible lasso problem with known factors.
This assumption has been discussed by many authors, and various sufficient conditions have been provided (e.g., Raskutti:2010 and loh2015regularized). The following lemma provides a simple sufficient condition for both this assumption and the restricted/sparse eigenvalue assumption.
As we described earlier, even if $p/n\to\infty$, requiring that the eigenvalues of the population covariance matrix are well-bounded is not a stringent condition. Note that this condition is imposed only on the factor-residuals, $U_{it}$, and that similar conditions on the population covariance matrix of factor residuals are typically imposed in the formal analysis of large approximate factor models.
The following conditions are imposed on the bootstrap weights.
The sub-exponential condition for the bootstrap weights enables us to bound many stochastic processes uniformly in $m\leq p$ and $t\leq T$. In our numerical studies, we follow mammen1993:bootstrap and use $w_i^x=\zeta_{1,i}^x/\sqrt{2}+((\zeta_{2,i}^x)^2-1)/2$ where $\zeta_{1,i}^x$ and $\zeta_{2,i}^x$ are independent standard normals and $x \in \{U,Y,D\}$.
Finally, we impose further regularity on the quality of estimation of the factors in the bootstrap data.
As with Assumption (ref), we show that
when $\widehat f_t^*$ is estimated using PCA in Appendix (ref) which allows direct verification of conditions involving $\Delta_F^*$. We handle the remaining terms by directly deriving the rates of convergences for each required term in Appendix (ref).
Under these additional conditions, we are able to verify that the confidence interval resulting from application the $k$-step bootstrap procedure has asymptotically correct coverage.
In this section, we discuss estimation of factors and factor residuals using principal components (PC).\footnote{We choose to focus on the PC estimator as a concrete example because it is relatively simple and is free of tuning parameters. One could consider other options which would also satisfy our assumed high-level conditions. For example, the weighted PC estimator (e.g. choi, BL13), can be more efficient than the standard PC estimator but requires additional tuning parameters for practical application.} We also provide low-level conditions under which the high-level conditions used in establishing Theorem (ref) are satisfied for PC.
Let $$ \tilde X=
_{pT\times n}, \quad \tilde \Lambda=
_{pT\times K},\quad \tilde F=
_{n\times K}, $$ and define $\tilde U$ similarly. The matrix form of the factor model is then $$ \tilde X=\tilde\Lambda \tilde F'+\tilde U, $$ where the individual and time effects have already been removed.
One of the most commonly used factor estimators is based on the PC of the $n\times n$ matrix $\tilde X'\tilde X$. Let $\widehat F$ denote the $n\times K$ matrix of the estimated factors. The columns of $\widehat F/\sqrt{n}$ are the eigenvectors of the first $K$ eigenvalues of $\tilde X'\tilde X/(npT)$. Let $V$ be the $K$ by $K$ diagonal matrix consisting of the first $K$ eigenvalues. Then the PC estimator estimates $\tilde F$ up to a $K\times K$ rotation matrix (e.g., SW02 and bai03) $H$ defined by $$ H=\frac{1}{npT}\tilde\Lambda'\tilde\Lambda \tilde F'\widehat FV^{-1}. $$ The factor estimator in the bootstrap sampling space is defined similarly with $\widehat F^*$ denoting the $n\times K$ matrix of the estimated factors whose columns are $\sqrt{n}$ times the first $K$ eigenvectors of $\tilde X^{*'}\tilde X^*/(npT)$. Finally, it is important to note that we do not need to estimate the factors but only need to estimate the space spanned by the factors for Theorem (ref) to hold.
We now present additional regularity conditions which are sufficient to verify that the PC estimator satisfies the conditions given in Assumption (ref). These conditions are standard for high-dimensional approximate factor models.
Assumption (ref) effectively implies that the factors do not load on a small number of series but rather are related to a large number of the available $X$-variables. The use of this assumption in high-dimensional factor models provides part of the motivation for the factor-lasso approach where at least some forms of association between factors that are not pervasive but instead load on only a few elements in $X$ and an outcome can be captured through the presence of the factor residuals in the equations of interest.
The left hand side of the third condition equals $\max_i \operatorname{Var}(\frac{1}{\sqrt{p}} \sum_{m=1}^p(\sqrt{T}\bar U_{i\cdot, m})^2)$. In addition, if we ignore the absolute value, then the left hand side of the fourth condition equals $\max_{im} \operatorname{Var} (\frac{1}{\sqrt{p}} (\sqrt{T}\bar U_{i\cdot ,k})(\sqrt{T}\bar U_{i\cdot, m}))$. Hence the third condition means that the variance of the standardized squared average should be bounded, and the fourth condition is slightly stronger than requiring $\max_{im} \operatorname{Var} (\frac{1}{\sqrt{p}} (\sqrt{T}\bar U_{i\cdot ,k})(\sqrt{T}\bar U_{i\cdot, m}))<C$. In the special case when $\{U_t\}$ is serially independent across $t$, all of the four conditions in Assumption (ref) can be directly verified under various notions of weak cross-sectional dependence.
The following proposition show that the high-level Assumptions (ref) and (ref) are satisfied by the PC estimator.
The conditions $|J|_0^4n=o(p^2T)$ and $|J|_0^2\log n=o(p)$ require lower bounds on the growth of $p$. These conditions differ from those used in the literature on inference in purely sparse high-dimensional, e.g. belloni2014inference, in that lower bounds on $p$ are not required in the purely sparse setting. These lower bounds arise since accurately estimating the unknown factors using PCA requires a large number of observed series. Indeed, the “average rate of convergence" is $$ \frac{1}{n}\sum_{i=1}^n\|\widehat f_i-H'\tilde f_i\|_2^2=O_P(\frac{1}{n^2}+ \frac{1}{nT^2}+\frac{1}{pT}), $$ where the product $pT$ is the dimension of $\tilde X_i$. In the special case $|J|_0=O(1)$, these conditions require $$ T \ll n \ll p^2T, \ \log ^{3}p=O(n), \ \text{and} \ \log n=o(p). $$ The results developed in this paper will thus be inappropriate in settings where $p$ is quite small relative to $n$. Of course, in the setting with $p$ small relative to $n$, a simple approach is to just use all of the available variables without dimension reduction.
Finally, though we have been assuming the number of factors, $K$, is known a priori, our procedure admits data-dependent methods (e.g., BN02 or AH) for selecting $K$. Under mild conditions such as those employed in BN02 or AH, $\widehat K$, the estimator for $K$, is consistent. All the preceding results hold can then be shown to hold following first-step estimation of $K$ by conducting the theoretical analysis conditional upon the event that $\widehat K = K$ and then arguing that the results asymptotically hold unconditionally as $P(\widehat K=K) \to 1$.
We now present simulation and empirical results in support of the formal analysis presented in the previous sections. The first simulation example is based directly on the PPFM given in ((ref))-((ref)). The second simulation example is based on a purely cross-sectional model that allows for instrumental variables estimation of the parameter on an endogenous variable in the presence of a low-dimensional set of instrumental variables (IVs) and a large number of potential control variables.\footnote{The formal development in the IV case with a small number of instruments is a notationally burdensome but straightforward extension of the results developed in this paper.} Following the simulation experiments, we then present results from two empirical applications. In the first, we apply the developed procedure to estimate the effects of gun prevalence on crime following cook:ludwig:guns using the data from BCHK:FE. In the second example, we apply the instrumental variables strategy of acemoglu:colonial to try to estimate the effect of institutions on growth.
In our first set of simulations, we report results for estimation and inference on $\alpha$ with data generated according to
with $n = 100$, $T = 10$, $K = 3$, and $p = 100$. We take $\epsilon_{it} \sim N(0,1)$, $\eta_{it} \sim N(0,1)$, and $U_{it} \sim N(0_{p},\Sigma_U)$ where $0_{p}$ is a $p \times 1$ vector of zeros, $\Sigma_U$ has $(r,s)$ element given by $[\Sigma_{U}]_{[r,s]} = .7^{|r-s|}$, and $\epsilon_{it}$, $\eta_{it}$, and $U_{it}$ are i.i.d. over $i$ and $t$ and jointly independent of each other. We generate unobserved individual-specific and time-specific heterogeneity by taking $n$ i.i.d. draws, one for each individual, $(g_i,\zeta_i,w_i) \sim N(0_{p+2},I_{p+2})$ where $I_{p+2}$ is a $(p+2) \times (p+2)$ identity matrix and taking $T$ i.i.d. draws, one for each time period, $(\nu_t,\mu_t,\rho_t) \sim N(0_{p+2},I_{p+2})$. The latent factors, $f_i$, are generated as i.i.d. draws from $ N(0_{K},I_K)$. The factor loading vectors $\xi_t$ and $\delta_{dt}$ and factor loading matrix $\Lambda_t$ are drawn independently over time with each entry generated as an independent draw from a standard normal random variable. The individual-specific, time-specific heterogeneity terms and factor loadings are drawn once, and the same values are used in each simulation replication.
We set the $j^{\text{th}}$ entry of $\theta$ and $\gamma_d$ as $\theta_j = \gamma_{d,j} = \frac{1}{j^2}$. $c_{\Lambda}$, $c_{\delta}$, $c_{\gamma}$, $c_{\xi}$, and $c_{\theta}$ are scalars that are set to alter the relative strength of $f_i$ and $U_{it}$ in each equation. We choose $c_{\Lambda}$ so that the average $R^2$ from the $p$ regressions of $X_{it,j}$ on $f_i$ is 0.5. We choose $(c_{\delta},c_{\gamma})$ so that the $R^2$ of the infeasible regression of $d_{it} - \zeta_i - \mu_t$ on $(c_{\delta}\delta_{dt})'f_i + U_{it}'(c_{\gamma} \gamma_d)$ is 0.7 and the factors account for 0%, 25%, 50%, 75%, or 100% of the explanatory power in this regression. We similarly choose $(c_{\xi},c_{\theta})$ so that the $R^2$ of the infeasible regression of $y_{it} - \alpha d_{it} - g_i - \nu_t$ on $(c_{\xi}\xi_{t})'f_i + U_{it}'(c_{\theta} \theta)$ is 0.7 and the factors account for 0%, 25%, 50%, 75%, or 100% of the explanatory power in this regression. Finally, we set $\alpha = 1$.
We compare the performance of the procedure developed in this paper to several benchmarks. Because we consider a design with $p < nT$, ordinary least squares of $y_{it}$ on $d_{it}$, $X_{it}$ and a full set of individual and time dummy variables is feasible (OLS). We also consider estimating $\alpha$ based on the assumption that confounding is entirely captured by latent factors. To implement this procedure, we extract factors, $\widehat f_i$, from $\tilde X_{it}$ by PCA as discussed in Section (ref). We then regress $y_{it}$ on $d_{it}$, $\widehat f_i$ interacted with a complete set of time dummy variables, and a full set of individual and time dummy variables to obtain the estimator for $\alpha$ (Factor). For our third procedure, we directly apply the fixed effects double-selection procedure of BCHK:FE which is appropriate for a sparse high-dimensional model with fixed effects (Double Selection). We then consider two ad hoc variants of the double-selection approach. In the first, we extract the first 20 principal components and interact these with a full set of time dummies. We then apply the fixed effects double-selection procedure of BCHK:FE to the data $(Y,D,X^*)$ where $X^*$ denotes the original $X$ variables augmented to include the interactions of principal components with time dummies (Double Selection F). The second ad hoc procedure extracts factors from $\tilde X_{it}$ by PCA. We then obtain estimates $\widehat U_{it}$ as in ((ref)) and apply the fixed effects double-selection procedure of BCHK:FE to the data $(Y,D,\hat{U}^*)$ where $\hat{U}^*$ denotes the matrix formed by combining $\widehat U$ with the interactions of principal components with time dummies (Double Selection U). Finally, we directly apply the factor-lasso approach outlined in this paper (Factor Lasso). We use the AH procedure to select the number of factors to use in obtaining the Factor, Double Selection U, and Factor Lasso results.
Figure (ref) gives simulation RMSEs for the estimator of $\alpha$ resulting from applying each procedure. The RMSEs are truncated at 0.1 for readability of the figure. The most striking feature of Figure (ref) is that only the proposed factor lasso procedure delivers uniformly good performance regardless of the relative strength of the factors and factor residuals in this simulation design. Each of the other procedures exhibits behavior that depends strongly on the exact strength of the factors in the different equations. In terms of RMSE, the factor-lasso procedure uniformly dominates regular OLS, Double Selection ignoring the factor structure, and the ad hoc procedure Double Selection F within the design considered. The factor-lasso estimator of $\alpha$ is outperformed by the pure factor model in the case where all of the explanatory power in the outcome equation is contained in the factors, which corresponds to the case where the pure factor model is correctly specified and there is no additional confounding based on the factor residuals, and the Double Selection U procedure when the factors have no explanatory power in the treatment (D) equation but all explanatory power in the Y equation. It is also important to note that the performance loss is small in these few cases where the factor lasso is outperformed. A final interesting point to note is that the conventional lasso-based double selection procedure is outperformed by the factor lasso even when the factors do not load in either the treatment or outcome equation. It seems likely that the loss in this case is due to the presence of the factors in the observed explanatory variables which leads to strong correlation among these variables. This strong correlation among the $X$'s is well-known to pose challenges for lasso-type estimators.
We report size of 5% level tests based on standard asymptotic approximations for each of the six procedures considered in Figure (ref) where the sizes are truncated at 0.3 for readability of the figure. In each panel, we report the rejection frequency of the standard t-test of the null hypothesis that $\alpha = 1$ with standard errors clustered at the individual level. The most striking feature of the figure is again the uniformly good performance of tests based on the proposed factor lasso procedure. Tests based on the factor-lasso procedure effectively control size, with size ranging between 3.3% and 5.3% across the design parameters considered in the simulation. This behavior is in sharp contrast to the other procedures considered which may have large size distortions depending upon exactly how large the relative contribution of the factors is in the $D$ and $Y$ equations. Importantly, this good behavior does not come at the cost of using an inferior estimator as evidenced by the RMSE results.
We conclude this discussion by looking at the performance of the k-step bootstrap. In Figure (ref), we report size of 5% level tests using the factor-lasso estimator and the asymptotic approximation provided in Theorem (ref), the k-step bootstrap, and a score bootstrap based on BCFH:Policy. The k-step bootstrap and asymptotic approximation have similar performance that keeps size close to the promised level. Interestingly, the score-based bootstrap that does not reestimate the factors or the lasso parts of the model exhibits mild size distortions across all of the design settings in this example.
We supplement the simulation results from the PPFM with additional simulations in a cross-sectional version of the model generalized to allow for an endogenous variable. Specifically, we generate data from the model
with $n = 100$, $K = 2$, and $p = 100$. Within this model, $d_i$ is an endogenous variable with coefficient of interest $\alpha$ and $z_i$ is an instrumental variable. We generate $\epsilon_{i} \sim N(0,1)$ and $\eta_{i} \sim N(0,1)$ with ${\mathrm{E}}[\epsilon_i\eta_i] = .8$ i.i.d. across $i$ and independent of all other random variables. We generate i.i.d. draws for $U_{i}$ as before, and $v_i \sim N(0,1)$ independently from $U_i$. We also generate $(\nu,\mu,\rho, \xi,\delta_d,\Lambda, c_{\Lambda}, c_{\gamma_z}, c_{\delta_d}, c_{\xi}, c_{\theta})$ as before. We set $\theta = \gamma_d = \gamma_z$ to be vectors with $j^{\text{th}}$ entry given by $\theta_j = \gamma_{d,j} = \gamma_{z,j} = \frac{1}{j^2}$. To control the strength of the instrument, we choose $(c_{\delta_z},c_{\gamma_z})$ so that the $R^2$ of the infeasible regression of $z_{i} - \zeta$ on $(c_{\delta_z}\delta_{z})'f_i + U_{i}'(c_{\gamma_z} \gamma_z)$ is 0.7 and the factors account for 50% of the explanatory power in this regression. We set $\pi$ so that the fraction of variation accounted for by $z_i$ in the regression of $d_i$ on $z_i$, $f_i$ and $U_i$ is 25%. Finally, we set $\alpha = 1$.
We again estimate $\alpha$ using six different IV procedures similar to those implemented in the previous simulation with one exception. As the number of features is equal to the sample size in these simulations, we consider an infeasible “oracle” estimator that estimates $\alpha$ from IV regression of $y_{i} - (c_{\xi} \xi)'f_i - U_{i}'(c_{\theta} \theta) - \nu$ on $d_{i} - (c_{\delta_d}\delta_{d})'f_i - U_{i}'(c_{\gamma_d} \gamma_d) - \mu$ using $z_{i} - (c_{\delta_z}\delta_{z})'f_i - U_{i}'(c_{\gamma_z} \gamma_z) - \zeta$ as instrument (Oracle). This estimator provides a type of best-case benchmark and allows us to ascertain that instruments are strong enough that the usual asymptotic approximation provides a reasonable approximation in the idealized scenario where one is able to perfectly remove the effect of confounding from all variables.
Figure (ref) gives simulation RMSEs for the estimator of $\alpha$ resulting from applying each procedure. The RMSEs are truncated at 0.1 for readability of the figure.\footnote{Theoretically, the MSE of the IV estimator does not exist in this context. We report root mean truncated squared error with a truncation point of 1.} Again, we see that the factor lasso procedure delivers good performance regardless of the relative strength of the factors and factor residuals in this simulation design. Each of the other procedures exhibits behavior that depends strongly on the exact strength of the factors in the different equations. It might be noted that the dominance of the factor-lasso estimator, in terms of RMSE, over the “Oracle” procedure is due to the definition of the oracle that we use which fully removes the variation in each variable due to factors and factor residuals even in situations in which some of these variables produce no confounding. For example, one need not remove the variation in the instruments due to the factors in cases where the factors have zero loadings in the outcome equation, but this variation is always removed due to the way we have defined the oracle model.
We report size of 5% level tests based on standard asymptotic approximations for each of the six procedures considered in Figure (ref) with size truncated at 0.3 for readability of the figure. In each panel, we report the rejection frequency of the standard t-test of the null hypothesis that $\alpha = 1$ using heteroscedasticity robust standard errors. Here, we see that the only procedure that uniformly controls size is the infeasible oracle. Among the feasible procedures, the proposed factor lasso approach performs relatively well in keeping size distortions small across the majority of combinations of relative strengths of the factors. In this case, we do see that the factor-lasso procedure suffers from reasonably large size distortions when the factors account for all of the confounding in the outcome equation and a moderate amount of counfounding in the treatment equation. We also see that the pure factor model controls size well in this case, but performs very poorly once all variation in the outcome equation is not due to the factors.
We again conclude by looking at the performance of the k-step bootstrap in Figure (ref). We see that there is a modest, but clearly visible, improvement from using the k-step bootstrap relative to the asymptotic approximation. The score based bootstrap, on the other hand, lines up reasonably well with the asymptotic approximation.
Overall, the results from the two simulation experiments are supportive of the asymptotic theory. We see that the factor-lasso approach delivers estimators with good properties relative to other feasible procedures that leverage either a pure factor structure or a pure sparse structure in partial factor model settings. We see that both point estimation properties, measured in terms of RMSE, and inferential quality, as measured by size of tests, are competitive or much better than the other procedures considered in our simulation design. The results also suggest that the proposed k-step bootstrap procedure works relatively well and may offer some gains relative to the asymptotic Gaussian approximation.
In this example, we follow BCHK:FE who build upon the work of cook:ludwig:guns and attempt to estimate the effect of gun prevalence on crime in a setting with a high-dimensional set of potential controls. As in BCHK:FE, we focus exclusively on trying to measure the effect of gun prevalence on homicide rates. An important difficulty with estimating the effect of gun prevalence in the United States is that exact gun-ownership numbers are difficult to obtain. Due to this difficulty, cook:ludwig:guns use the fraction of suicides committed with a firearm (abbreviated FSS) within a county to proxy for county-level gun ownership rates. cook:ludwig:guns provide a series of arguments and evidence from secondary data sources supporting the claim that FSS provides a useful proxy for gun ownership. For the analysis in this paper, we simply take it as given that estimating a causal effect of FSS on crime measures is worthwhile and abstract from any further measurement or data issues surrounding the use of this proxy.
Both cook:ludwig:guns and BCHK:FE estimate linear fixed effects models of the form
where $g_i$ and $\nu_t$ are treated as parameters to be estimated, $X_{it}$ are control variables, and $Y_{it}$ is one of three dependent variables: the overall homicide rate within county $i$ in year $t$, the firearm homicide rate within county $i$ in year $t$, or the non-firearm homicide rate within county $i$ in year $t$. cook:ludwig:guns use the four variables percent African American, percent of households with female head, nonviolent crime rates, and percent of the population that lived in the same house five years earlier as their set of controls $X_{it}$. BCHK:FE maintain the assumption of approximate sparsity and employ their variable selection approach using a much larger set of potential controls generated by taking variables compiled by the US Census Bureau as $X_{it}$. Their variables include county-level measures of demographics, the age distribution, the income distribution, crime rates, federal spending, home ownership rates, house prices, educational attainment, voting patterns, employment statistics, and migration rates along with interactions of the initial (1980) values of all control variables with a linear, quadratic, and cubic term in time.
Rather than adopt the approximately sparse model in ((ref)), we employ the PPFM, ((ref))-((ref)), and factor-lasso approach to estimate $\alpha$ using 909 variables in $X_{it}$ constructed as in BCHK:FE.\footnote{The exact identities of the variables are available upon request. The data is from the U.S. Census Bureau USA Counties Database, http://www.census.gov/support/USACdataDownloads.html.} The PPFM model seems very appropriate for this data as it directly incorporates a mechanism to accommodate the concern that there are features of counties that are not directly observed, the $f_i$, but are related to the evolution of the outcome and treatment variable of interest, which is captured by the time-varying factor loadings. Obviously, exclusion of these factors would then lead to omitted variables bias in any estimator of $\alpha$ that fails to capture them. Concern about the existence of such factors is common in empirical applications involving aggregate panel data.
The key assumption that we leverage to allow us to simply accommodate these latent factors is that the same correlated unobserved factors that lead to confounding are related to the evolution of other observed county-level aggregates and that we have access to a large number of these auxiliary aggregates. While this key assumption is strong, the PPFM also naturally provides some robustness to the presence of shocks ($U_{it}$) that are related to movements of the observed $X_{it}$ series as well as movements in the variable of interest and outcome. Such shocks may be motivated, for example, by the factor structure being misspecified, by the presence of variables that are not strongly related to factors but are confounded with the treatment and outcome, and simply by the presence of local shocks not captured by the factors that are related to the observed series.
We present estimation results in Table (ref) with results for each dependent variable presented across the columns and rows corresponding to different estimation approaches. As a baseline, we report numbers taken directly from the first row of Table 3 in cook:ludwig:guns in the first row of Table 1 (“Cook and Ludwig (2006) Baseline”). We report results obtained from our data the remaining rows.\footnote{All results are based on weighted regression where we weight by the within-county average population over 1980-1999.} For these results, we first report the point estimate and estimate of the asymptotic standard error obtained by clustering by county. Immediately below these results, we report the 95% confidence interval obtained from applying the k-step bootstrap procedure in brackets. The rows labeled “Post Double Selection” apply the procedure of BCHK:FE. The rows labeled “Factor” are based on a pure factor model; the rows labeled “Factor-Lasso” use the proposed factor-lasso procedure. All factors are estimated using PCA and the number of factors is selected using AH.
We see that the estimates and inferential statements produced for the firearm homicide rate (“Gun”) and the non-firearm homicide rate (“non-Gun”) are broadly consistent with each other. In all cases, there is a fairly large positive point estimate for the effect on the firearm homicide rate with corresponding 95% confidence intervals that exclude zero, suggesting positive association between the used measure of gun prevalence and gun homicides. For the non-firearm homicide rate, all point estimates are negative and modest and confidence intervals include both positive and negative values. The broad results for the overall homicide rate (“Overall”) are slightly more mixed. The baseline results for cook:ludwig:guns and results from a pure factor model suggest a strongly significant, positive effect of gun prevalence on the overall homicide rate. Assuming sparsity and applying BCHK:FE yields a positive estimate of the effect which is statistically insignificant at the 5% level. Finally, the factor-lasso estimator is similar in magnitude to the sparsity-based estimator but borderline significant at the 5% level using the bootstrap confidence interval.
A more interesting comparison can be made by looking more closely and considering the variable and factor selection results. The “Post Double Selection” procedure ends up selecting three variables for estimating the effect on overall homicide rates, three variables for gun homicide rates, and two variables for non-gun homicide rates. The pure factor model uses one factor. The factor-lasso approach then uses one factor in all cases but selects eight additional variables for estimating the effect on the overall homicide rate, eight additional variables for the gun homicide rate, and five additional variables for the non-gun homicide rate. These results suggest that the “Post Double Selection” and “Factor” results may be based on models that fail to adequately capture the effect of potential confounds. We also see that the “Factor” estimates are substantially shifted away from the “Factor Lasso” estimates relative to standard errors and that the factor-lasso estimates are the most precise in the sense of having the shortest confidence intervals. Both findings are consistent with the asymptotic theory and with the simulation results.
We revisit the example considered in acemoglu:colonial. acemoglu:colonial are interested in the parameter $\alpha$ in a structural model of the form
based on aggregate country level data where “Protection from Expropriation” is a measure of the strength of individual property rights that is used as a proxy for the strength of institutions and $x_i$ is a set of variables that are meant to control for geography. acemoglu:colonial adopt an IV strategy where they instrument for institution quality using early European settler mortality to estimate $\alpha$ as institutions are clearly potentially endogenous. They point out that their instrument would be invalid if there were other factors that are highly persistent and related to the development of institutions within a country and to the country's GDP. A leading candidate for such a factor that they discuss is geography. To address this possibility, acemoglu:colonial control for the distance from the equator in their baseline specifications and consider different sets of geographic controls such as continent dummies within their robustness checks.\footnote{E.g. acemoglu:colonial Table 4.}
There are, of course, many other ways to measure geography besides distance to the equator or continent where a country is found. Rather than ex ante choose a small number of variables to proxy for geography, we put a large number of variables that potentially capture geography in $x_i$ and then use the data to reduce dimension. Specifically, we consider dummies for Africa, Asia, North America, and South America as well as longitude, renewable water, land boundary, land area, amount of coastline, territorial seas, amount of arable land, average temperature, average high temperature, average low temperature, average precipitation, elevation of highest point, elevation of lowest point, fraction of area that is low-lying, latitude, and spherical distance from London.
We adapt the analysis of acemoglu:colonial to the present setting by considering estimation of a partial factor instrumental variables model
using our 20 geography measures as $x_i$ and the 64 countries from the original acemoglu:colonial data. The factor-lasso approach seems quite sensible in this setting. Each of the observed geography measures could reasonably be taken as a noisy proxy for a country's geography. This relationship is likely to be complicated and uneven with the chief features leading to association between the geography proxies plausibly being only weakly related to the notions of geography that are important predictors of mortality and institutions. The factor-lasso approach, by allowing a small number of elements of $U_i$ to enter the equation of interest in addition to any common geography factors, readily accommodates this latter possibility in a parsimonious, data-dependent way.
We report estimation results for the first stage coefficient on the instrument in Table (ref). We report results from the factor-lasso approach in the row “Factor-Lasso”. For comparison, we also report results from a few natural alternative models. The row labeled “Latitude” uses the single variable distance from the equator to control for geography as in the baseline results from acemoglu:colonial. We report results from applying OLS using all 20 available geographic controls without dimension reduction in the row labeled “All Controls.” We apply the double selection approach of belloni2014inference which would be appropriate if the relationship between geographic controls and the variables of interest were well-approximated by a sparse linear model in “Double Selection.” Finally, “Factor” reduces dimension through positing a conventional factor model. All factors are estimated using PCA with number of factors selected by applying the procedure from AH.
The first-stage results using only the latitude control suggest there is a fairly strong relationship between the instrument and endogenous variable if latitude is a sufficient control for geography. The first stage F-statistic using just latitude is 10.9 which many would take to indicate that the instrument is sufficiently strong to identify the effect of interest.\footnote{A benchmark that is commonly used in the applied literature to assess whether there is sufficient variation in the instrument to identify the effect of interest is to compare the first stage F-statistic to 10, with smaller values indicating weak identification.} The results change in a potentially substantive way after allowing for the possibility that geography is not adequately captured by latitude. For each of the remaining approaches considered, the first-stage F-statistic drops substantially below 10, with all methods besides applying the pure factor model returning first-stage coefficients that are statistically insignificant at the 5% level.
One might dismiss the lack of significance after including all controls without dimension reduction as it seems likely that a model with 20 covariates in addition to the variables of interest and only 64 observations is overfit. The next strongest result is from the pure factor model which makes use of a single extracted component and produces a first-stage F-statistic of 7.5. As evidenced in the simulation example, inference results based on a pure factor model may be highly misleading when elements of $U_i$ also have explanatory power. It is then interesting that the double-selection approach and the factor-lasso approach deliver almost identical results indicating a weak association between the endogenous variable and instrument after controlling parsimoniously for geography. The double-selection procedure selects four variables\footnote{These variables are the Africa dummy, average temperature, average high temperature, and amount of arable land.}, and the factor-lasso approach uses one factor and two additional variables.\footnote{The two selected variables in addition to the factor are the Africa and Asia dummies.} One might take this to mean that the four variables selected in the double-selection procedure approximately capture the same information as the single factor and two variables used in the factor-lasso results. In either case, the results suggest that, at best, identification of the structural effect of institutions as measured by “Protection from Expropriation” using settler mortality as instrument is weak after geography is controlled for in a parsimonious, data-dependent way. Given this apparent weak identification, we do not report second stage estimates of the structural effect.\footnote{We note that it would be straightforward to adapt the weak-identification robust procedure of ch:WeakId to the present setting. We do not pursue this extension for brevity.}
We believe the two empirical examples illustrate the potential applicability of partial factor models and the associated factor-lasso approach in applied economics. The model provides a natural generalization to standard factor models and sparse high-dimensional models and seems appropriate for many economic applications, especially those that make use of aggregate panel or cross-sectional data. The results in the first example based on cook:ludwig:guns roughly line up with the original results, though they demonstrate the potential for efficiency gains from adopting the methods developed in this paper. In the second example, we draw substantively different conclusions about the strength of identification than one would draw following the approach in acemoglu:colonial due to the ability to control more flexibly for the leading candidate for confounding. Overall, the results suggest that application of the proposed methods may usefully complement the sensitivity analyses performed in empirical economics and also have the potential to strengthen the plausibility of any conclusions drawn.