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,351 characters · 13 sections · 77 citation commands
Arellano-Bond LASSO Estimator for Dynamic Linear Panel Models$^*$
Panel data involve observations collected for cross-sectional units ($i=1,\ldots,N$) over multiple time periods ($t=1,\ldots, T$). Models for panel data are commonly used in economics and other social sciences because they allow researchers to control for unobserved unit and time heterogeneity and account for unit-level dynamics. These models have multiple applications, including evaluating job training and minimum wage regulations in labor economics, studying household consumption and economic growth in macroeconomics, estimating demand models for products in microeconomics, and analyzing payout policies and investment decisions in corporate finance. See bond2002dynamic for a review of methods and applications of dynamic panel data models.
holtz1988estimating introduced instrumental variable (IV) estimation in dynamic panel data models without strict exogeneity. Following this line of research, the Arellano-Bond estimator (AB) has become the most widely used fundamental method for panel models arellano1991some,arellano1995another. It applies to dynamic linear models that include lagged dependent variables and predetermined covariates as explanatory variables and unobserved unit and time fixed effects. After taking first differences or forward orthogonal deviations to remove the unit fixed effects, AB constructs moment conditions using sufficiently lagged dependent variables and covariates as instruments and applies the Generalized Method of Moments (GMM) to estimate the model parameters. However, AB might be severely biased in long panels. The problem arises because the number of moment conditions grows with the square of $T$, $T^2$, leading to many instrument bias caused by the large degree of overidentification in the GMM problem newey2004higher. More precisely, AB has an asymptotic bias of order $T/N$, which might not be negligible compared to $1/\sqrt{NT}$, the size of the stochastic error, when the time dimension $T$ is sufficiently large relative to $N$ alvarez2003time. This problem causes bias in estimators and undercoverage of confidence intervals.\footnote{The original motivation of AB was specification testing. We abstract from this issue because our method does not naturally lead to a specification test. Readers interested in this topic are referred to the existing literature on specification testing in IV settings with many instruments, including anatolyev2011specification, hansen2008estimation, chao2014testing, and shi2025testing.}
We address the bias issue of AB within long panels through a two-step method. After removing the unobserved unit fixed effects by forward orthogonal deviations, we first select the most informative moment conditions, followed by applying a linear instrumental variable estimation using instruments derived from these conditions. Specifically, as the number of AB's moment conditions varies across time periods, we perform a moment selection procedure on cross-section data for each time period separately. We utilize the least absolute shrinkage and selection operator (LASSO) of tibshirani1996regression as our selector, given the naturally sparse structure of the moment conditions under appropriate weak temporal dependence conditions.
Several moment selection methods have been previously established in other contexts. For instance, donald2009choosing introduced an alternative method for selecting instruments based on asymptotic mean squared error calculations. belloni2012sparse described a similar approach to select optimal instruments using LASSO in cross-section instrumental variable models. Other methods for constructing optimal instruments through model averaging have been proposed by kuersteiner2010constructing and okui2011instrumental. luo2016selecting expanded the LASSO selector to nonlinear GMM settings with many potential moments, noting its computational advantages over donald2009choosing. Although some studies employ the AB estimator to underpin their analyses or conduct simulations newey2009generalized, none addresses the AB estimator directly, mainly because they assume independent and identically distributed data.
LASSO utilizes the $\ell_1$-norm to select the moment conditions. To show the validity of the selector, we build on theoretical results achieving near-oracle rates for LASSO and related estimators in bickel2009simultaneous, BCH2011 and belloni2013least. An additional complication comes from the high dimensionality of the problem. As we mentioned above, the number of moment conditions grows with $T^2$. Moreover, we allow for high-dimensional covariates including multiple lags of the dependent variable and other strictly exogenous covariates, which are becoming common in modern applications, and unit and time fixed effects, whose number grows with the two dimensions of the panel. While our method does not suffer from shrinkage and model selection biases because the moment conditions of the second step are Neyman-orthogonal with respect to the parameters estimated in the first step, it can still be subject to over-fitting bias, specially in the presence of high-dimensional covariates. We deal with this problem by combining the two steps of our procedure using cross-fitting DML. Thus, we partition the panel in two parts. We select the moment conditions in the first part and estimate the parameters in the second part. Then, we repeat the procedure reversing the roles of the two parts and aggregate the results by averaging the estimates of the model parameters from the two orderings. Cross-fitting does not change the large sample properties of the estimator in this case because the use of forward orthogonal deviations to remove the unit fixed effects attenuates the dependence between the two stages by spreading the transformed error over many future periods. We nevertheless find small sample improvements in various simulation settings.\footnote{In a previous version of the paper, we used first differences to remove the unit fixed effects. In that case, cross-fitting improved the rate condition for deriving the asymptotic distribution of the estimator by a factor of $\sqrt{T}$. The improvement arises because cross-fitting removes the dependence between the generated errors of the selected instruments and the transformed errors in the main regression by conditioning on the sub-sample used for moment selection.}
There is an extensive recent literature on panel data with large $T$, including dynamic linear models. alvarez2003time studied the properties of AB and other estimators in long panels. They showed that AB exhibits asymptotic bias when $T/N$ tends to a constant. moral2013likelihood, moral2019dynamic and alvarez2022robust developed likelihood-based alternatives to the AB estimator and chen2019mastering proposed a debiasing method based on applying the split-panel idea of dhaene2015split to the cross-section dimension of the panel, whereas nickell1980correcting, kiviet1995bias, hahn2002asymptotically and chudik2018half developed alternatives based on bias corrections of the fixed effects estimator. We compare our method to these alternatives in numerical simulations. okui2009optimal proposed a method to select instruments by characterizing the mean squared error of the one-step AB estimator that uses the matrix of second moments of the instruments as weighting matrix in models with homoskedastic errors and strictly exogenous regressors. In the same setting, carrasco2024regularized developed a version of the one-step AB estimator that regularizes the weighting matrix building on carrasco2012regularization. Our instrument selection method is different and applies to both one-step and two-step AB estimators. It also does not rely on homoskedastic errors and allows for weakly exogenous covariates. cheng2023weight have recently proposed a minimum distance estimator that combines graphical LASSO to regularize the weighting matrix with cross-fitting to reduce bias. While the steps of the procedure are similar to ours, the setting and regularization method are different.
We make four main theoretical contributions. First, we show that the moment conditions of AB exhibit an approximately sparse structure under suitable temporal dependence conditions (see, e.g., Proposition (ref)) and propose a LASSO version of AB that fruitfully exploits this structure. In particular, the effective dimension of the non-zero coefficients in the first step estimation of the selected instruments at each time period $t$ is the minimum of $\log N$ and $t$, which is very “low" relative to the cross-sectional size $N$. Second, we propose a cross-fitting procedure based on sample splitting that reduces finite-sample overfitting bias and opens the door to the use of machine learning methods other than LASSO. Unlike the existing high-dimensional IV literature, which typically focuses on cross-sectional settings, our framework is developed for dynamic panel models, where the relevant moment conditions may vary across time periods. We deal with this issue by performing moment selection separately for each time period and then aggregating the information over both $i$ and $t$ in the final instrumental variable estimator. The corresponding theory accommodates temporal dependence in the outcome and covariates, as well as the presence of time effects. Third, we show that first differences and forward orthogonal deviations have different properties in our setting. In particular, under suitable conditions, forward orthogonal differences yields more efficient estimators than first differences, and its large-sample properties are not affected by the use of cross-fitting. Fourth, we consider models with high-dimensional covariates that are becoming common in modern applications; e.g., klosin2022estimating,semenova2017estimation,gao2024robust,takeshima2025detecting. Here, we achieve moment selection by constructing orthogonal conditions and employing regularized GMM with a sparse weighting matrix via the Dantzig selector, in order to partial out the effect of high-dimensional nuisance parameters.
The finite sample properties of our method are illustrated through comprehensive Monte Carlo simulations, where we compare it with alternative approaches such as likelihood-based, debiased fixed effects, and debiased AB estimators. Lastly, we report results from an empirical application to the effect of the opening of K-12 schools and other policies on the spread of COVID-19 using a panel of $2,510$ US counties over $32$ weeks, extracted from the dataset used in chernozhukov2021association. We estimate a panel regression model with rich dynamics, incorporating four lags of the dependent variable and several predetermined covariates. Due to the large number of instruments ($m=3,375$), the small bias condition for AB of chen2019mastering, $m^2/(NT) \to 0$, fails for the AB estimator ($m^2/(NT) \approx 168$).\footnote{In this calculation $T=27$ because the first 5 observations are used as initial conditions of the dynamic model.} Compared to AB, our method finds that policies such as the opening of K-12 schools, stay-at-home orders, and banning gatherings have more muted absolute long-run effects, and college visits have substantially smaller effects in both the short and long run.
For a column vector $v=(v_1,\ldots,v_d)^\top\in \mathbb{R}^d$ and a constant $r\geq 1$, we denote $|v|_r=(\sum_{i=1}^d |v_i|^r)^{1/r}$ and $|v|_{\infty}=\max\limits_{1\leq i\leq d}|v_i|$. Define $|v|_0$ as the zero norm, i.e. the number of non-zero coordinates. For a matrix $A=(a_{ij})_{1\le i\le m, 1\le j\le n}$, we define $|A|_{\max}=\max\limits_{1\leq i\leq m,1\leq j\leq n}|a_{ij}|$, $|A|_{1}=\max\limits_{1\leq j\leq n}\sum_{i=1}^m|a_{ij}|$, $|A|_{\infty}=\max\limits_{1\leq i\leq m}\sum_{j=1}^n|a_{ij}|$, and $|A|_{1,1} = \sum_{i=1}^m\sum_{j=1}^n|a_{ij}|$. For a random variable $X_{it}$, we say $X_{it}\in\mathcal L^r$ if $\lVert X_{it}\rVert_r\stackrel{\mathrm{def}}{=}(\mathop{\mbox{\sf E}}|X_{it}|^r)^{1/r}<\infty$ for some $r>0$, and define the sub-Gaussian norm as $\|X_{it}\|_{\psi_{1/2}} = \inf\{s >0: \mathop{\mbox{\sf E}}\operatorname{exp}(X_{it}^2/s^2)\leq 2\}$, where $\mathop{\mbox{\sf E}}$ denotes the expectation conditional on the unit and time effects. We denote the limit cross-sectional average by $\bar\mathop{\mbox{\sf E}}(\cdot)$, that is $\bar\mathop{\mbox{\sf E}}(X_{it}) = \lim\limits_{N\to\infty}N^{-1}\sum_{i=1}^N\mathop{\mbox{\sf E}}(X_{it})$, provided that the limit exists.\footnote{Throughout the paper, we assume that the relevant limits exist, whenever we take expectations.} Given two sequences of positive numbers ${a_n}$ and ${b_n}$, we write $a_n\lesssim b_n$ (resp. $a_n\asymp b_n$) if there exists $C>0$, which does not depend on $n$, such that $a_n/b_n\le C$ (resp. $1/C\le a_n/b_n\le C$) for all large $n$. For a sequence of random variables ${X_n}$, we use the notation $X_n\lesssim_\mathrm{P} b_n$ to denote $X_n=\mathcal{O}_{\mathrm{P}}(b_n)$. For two real numbers, set $x\vee y=\max(x,y)$ and $x\wedge y=\min(x,y)$.
The rest of the paper is organized as follows. Section (ref) introduces the model and estimators. Section (ref) presents the main theoretical results. Sections (ref) and (ref) report the results of the simulation study and empirical application, respectively. Section (ref) contains some concluding remarks. An appendix collects the deferred proofs of the theoretical results, additional results for the analysis that uses first differences, instead of forward orthogonal deviations, to remove the unobserved unit effects, and supplementary results for the simulation study.
Let $\{(Y_{it},D_{it},C_{it}): 1 \leq i \leq N, 1 \leq t \leq T\}$ be a panel dataset, where $i$ and $t$ index unit and time period, respectively. $Y_{it}$ is a scalar outcome or response variable, $D_{it}$ is the policy variable or treatment of interest, and $C_{it}$ is a vector of covariates of fixed dimension including, for example, $Y_{i,t-1}$ and other treatments. To measure the effect of $D_{it}$ on $Y_{it}$, we consider a dynamic linear panel model:
where $\theta^0$ is the parameter of interest, $\alpha_i$ is an unobserved unit effect, $\gamma_t$ is an unobserved time effect, and $\varepsilon_{it}$ is an idiosyncratic error with zero mean and constant variance. We might also be interested in functions of $\theta^0$ such as long-run effects in dynamic models that include lags of the dependent variable as covariates. We refer to the empirical application in Section (ref) for an example.
In the theoretical analysis, we shall treat the unobserved unit and time effects as fixed parameters. This is equivalent to conditioning on the realization of all these effects.\footnote{Due to this conditioning, all probability statements should be qualified with almost surely. We shall omit this qualifier for notational convenience.} We assume that $\{(D_{it},C_{it},\varepsilon_{it}): 1 \leq t \leq T\}$ are independent over $i$, and $\varepsilon_{it}$ is an uncorrelated sequence over $t$. In addition, we assume that the treatment and covariates in $X_{it}$ are predetermined with respect to $\varepsilon_{it}$ in the sense that $$\mathop{\mbox{\sf E}}(X_{is}\varepsilon_{it}) = 0, \text{ for all } 1 \leq s \leq t \leq T.$$
We remove the unobserved effects by taking forward orthogonal deviations (FOD) over time and demeaning all the variables at the unit level, namely
where $\Delta \widetilde Z_{it} = \Delta Z_{it} - \sum_{j=1}^N \Delta Z_{jt}/N$, $\Delta Z_{it} = c_t\{Z_{it} - \sum_{s=1}^{T-t} Z_{i,t+s}/(T-t)\}$, and $c_t = \sqrt{(T-t)/(T-t+1)}$, for $Z_{it} \in \{Y_{it}, X_{it}, \varepsilon_{it}\}$.\footnote{In the previous version of the paper chernozhukov2024arellano, we used first differences instead of FOD to remove the unobserved effects. Appendix (ref) contains some of the results for the resulting estimators.} The transformed error, $\Delta \widetilde \varepsilon_{it}$, is an uncorrelated sequence over $t$ and satisfies the moment conditions $$ \mathop{\mbox{\sf E}} (X_{is}\Delta \widetilde \varepsilon_{it}) = 0, \text{ for all } 1\leq s \leq t \leq T-1. $$
AB uses these moment conditions to construct a GMM estimator of $\theta^0$.\footnote{More precisely, the version of AB that uses moment conditions in FOD is from arellano1995another.} It should be noted that AB is biased when $T$ is large due to the large number of moment conditions, i.e. $m = \mathcal{O}(T^2)$; see, e.g., newey2004higher for more discussion. We propose an alternative estimator that is computationally simple and has lower bias when $T$ is large. It is based on the application of LASSO to select the most informative moment conditions to estimate the parameters. Thus, the estimator has two stages. It first selects moment conditions using LASSO, and then estimates the parameters of interest by instrumental variables, with the predicted values of the endogenous regressors obtained from the selected moment conditions serving as instruments. We name the new estimator as AB-LASSO as a shorthand for Arellano-Bond LASSO estimator.
AB is an instrumental variable estimator. Its bias comes from overfitting because the same observations are used to project the endogenous regressors on the instruments and to estimate the parameters phillips1977bias,angrist1995split,angrist1999jackknife. The order of the bias is $m/n$, where $m$ is the number of moment conditions and $n$ is the sample size. In the case of AB, $m=\mathcal{O}(T^2)$ and $n=NT$, so that the order of the bias is $T/N$. The order of the sampling noise is $n^{-1/2} = (NT)^{-1/2}$, so that the small bias condition of chen2019mastering is $m/n^{1/2} \to 0$ or equivalently $m^2/n = T^3/N \to 0$. Our proposed AB-LASSO estimator reduces the overfitting bias by selecting moment conditions and by using the FOD transformation. Up to logarithmic terms, the small bias condition for AB-LASSO becomes $\max\limits_{1\leq t\leq T-1}\sqrt{s_t^*/N} \to 0$ ($s_t^*$ is the dimension of effective instruments for each $t$). When $s_t^*$ is moderately large relative to $N$, AB-LASSO might still exhibit small sample or higher order bias. To reduce this bias, we develop a sample-splitting procedure over the cross-section dimension following the idea of the split-sample IV estimator of angrist1995split. We name the version of AB-LASSO with sample splitting and cross-fitting as AB-LASSO-SS.
In many empirical panel applications, researchers augment dynamic specifications with a rich set of additional covariates to strengthen identification and improve robustness. Examples include dynamic treatment effect models with many policy indicators and firm-level panels with extensive macroeconomic controls. In modern applications, the number of such controls can be large relative to the sample size, especially when flexible specifications or many interaction terms are considered. To reflect this practice, we extend the basic model in (ref) by including additional covariates:
where $X_{2,it}$ is a possibly high-dimensional vector of covariates that is independent over $i$ and satisfies $\mathop{\mbox{\sf E}}(\Delta\widetilde X_{2,it} \Delta\widetilde\varepsilon_{it}) = 0$. The leading cases of such covariates are strictly exogenous variables with respect to $\varepsilon_{it}$. Denote the dimension of $X_{2,it}$ by $d_2$. The high-dimensional case arises when $d_2$ is large relative to the sample size $n=NT$ such that it is more appropriate to treat $d_2$ as increasing in the asymptotic analysis.
When the dimension of the additional covariates is small, their effects can be removed using standard projection arguments together with the unit and time fixed effects. However, when the number of covariates is large relative to the sample size, classical partialling-out is no longer feasible. In particular, direct estimation of the full parameter vector by least squares becomes ill-posed, and naive regularization (e.g., plugging in LASSO estimates) introduces shrinkage bias that contaminates inference on the low-dimensional parameters of interest. To address this issue, we treat the coefficients on the high-dimensional covariates as nuisance parameters and focus inference on a small set of structural parameters. Our strategy is to construct orthogonal moment functions that are insensitive, to first order, to regularization error in the estimation of the nuisance parameters. The key idea is to build time-specific orthogonalized instruments for the endogenous regressors that are strongly correlated with the components of interest and orthogonal to the high-dimensional covariates. This orthogonality ensures that small estimation errors in the nuisance component do not affect the asymptotic distribution of the estimator of the parameters of interest.
To formally explain this partialling-out procedure, it is convenient to rewrite the extended model (ref) as: $$Y_{it} = \alpha_i + \gamma_t + X_{1,it}^{\top}\theta^0_1 + X_{2,it}^{\top}\theta^0_2 + \varepsilon_{it}, \quad X_{1,it} := (D_{it}, C_{it}^\top)^{\top},$$ where $\theta_1^0\in\mathbb{R}^{d_1}$ and $\theta_2^0\in\mathbb{R}^{d_2}$. Denote $d=d_1+d_2$, where $d_1$ is fixed and $d_2$ is growing with $n$. Assume the sparsity assumption $|\theta_2^0|_0=\mbox{\tiny $\mathcal{O}$}(n)$. The moment functions are given by $$g_{it}(\theta_1,\theta_2) = \mathop{\mbox{\sf E}}\big\{(\Delta\widetilde Y_{it} - \Delta\widetilde X_{1,it}^\top\theta_1 - \Delta \widetilde X_{2,it}^\top\theta_2)U_{it}\big\},$$ where $U_{it}=(U_{it}^{0\top}, \Delta \widetilde X_{2,it}^\top)^\top$, and $U_{it}^0$ ($d_1\times1$) contains the most informative IVs for $\Delta\widetilde X_{1,it}$.
To construct orthogonalized instruments in the high-dimensional setting, we seek a time-specific weighting matrix $\mathcal W_t$ ($d\times d_1$) such that the transformed instruments $\mathcal W_t^\top U_{it}$ remain strongly correlated with $\Delta \widetilde X_{1,it}$, to preserve identification of $\theta_1^0$, and are orthogonal to $\Delta \widetilde X_{2,it}$, so that the influence of the high-dimensional nuisance component is removed. Formally, we aim to choose $\mathcal W_t$ so that $\mathop{\mbox{\sf E}} \big\{\Delta \widetilde X_{2,it} (\mathcal W_t^\top U_{it})^\top \big\} = 0$, while ensuring that $\mathop{\mbox{\sf E}} \big\{ \Delta \widetilde X_{1,it} (\mathcal W_t^\top U_{it})^\top \big\}$ has full rank $d_1$. This problem can be solved by Dantzig selector:
where we have replaced $U_{it}^0$ by the LASSO predictions, i.e. $\widehat U_{it}=(\widehat {\Delta \widetilde X_{1,it}^\top}, \Delta \widetilde X_{2,it}^\top)^\top$, and $\mathbf I_{d\times d_1}$ represents the $d\times d_1$ sub-matrix of the $d\times d$ identity matrix. Then, the instrument for $\Delta\widetilde X_{1,it}$ is $\widehat{\mathcal W}_{t}^{\top}\widehat{U}_{it}$ and the estimator of the parameters of interest becomes:\footnote{When $X_{2it}$ includes second or higher lags of the dependent variable the summation over $t$ in (ref) needs to be modified to include only the observed values of $\widehat U_{it}$. See Remark (ref) for a related discussion.}
Conceptually, in the low-dimensional case this orthogonalization reduces to a standard projection onto the orthogonal complement of the nuisance score. In the high-dimensional case, however, the projection cannot be computed exactly. We therefore approximate it by solving a constrained $\ell_1$-minimization problem that delivers a sparse weighting matrix. This step can be viewed as constructing an approximate Neyman-orthogonal score tailored to the panel structure of the model. The resulting estimator achieves valid inference for $\theta_1^0$ even when the dimension of the additional covariates grows with the sample size, provided the nuisance parameters are sufficiently sparse.
In this section, we present the theoretical foundation of the proposed estimator. We begin with the basic model (ref), which is a special case of the general model (ref), and establish the main results for this setting before turning to the extended model. Throughout this section, we impose the following conditions on the data generating processes.
We allow for overlap in the innovations $\xi_{it}$ and $\zeta_{it}$, as long as the exogeneity conditions specified in Section (ref) are satisfied, i.e., $\mathop{\mbox{\sf E}}(X_{is}\varepsilon_{it}) = 0$, for all $1 \leq s \leq t$. Note that Assumption (ref) allows for certain forms of non-stationarity in the process for $X_{it}$ such as unit-specific deterministic time trends and unrestricted time effects. For example, the policy variables of the empirical application in Section (ref) can be modeled as unit-specific deterministic time trends. It does rule out, however, other forms of stochastic non-stationarity such as unit roots. The following definition, along with Assumptions (ref) and (ref)(i) below, adapts the functional dependence measure proposed by wu2005nonlinear for stationary time series processes to heterogeneous panel data processes.
Example (ref) provides an example of a univariate heterogeneous linear process that satisfies Assumption (ref)(i). In Example (ref), we verify Assumptions (ref)--(ref) for a basic panel AR(1) model.
The m.d.s.\ condition in Assumption (ref)(ii) aligns with the standard large $T$ panel literature; see, for example, alvarez2003time and arellano2003panel. Using standard techniques such as $m$-dependent approximation and blockwise time series analysis, as in chen2022inference, the setting can be generalized. Additionally, for practitioners, we note that the sub-Gaussian conditions in Assumption (ref)(iii) rule out heavy-tail distributions for $X_{it}$ and $\varepsilon_{it}$. However, this assumption is not critical for our analysis. They can be relaxed to a polynomial tail conditions with more demanding rate assumptions.
We will first demonstrate the consistency property of the LASSO estimator $\widehat\Pi_t$, which is obtained in step 1 of AB-LASSO by (ref). For this purpose, a few definitions and assumptions are introduced as follows.
For each component $W_{it}$ of $\Delta\widetilde X_{it}$, denote the $N\times1$ vector containing $(W_{it})_{i=1}^N$ by $\boldsymbol{W}_t$. Recall that $V_{it} = (1,X_{i1}^\top,\ldots,X_{it}^\top)^\top$. Denote the dimension of $V_{it}$ as $m_t$, which is the number of instruments for each time period $t$. We further stack $V_{it}^\top$ by rows for all $i=1,\ldots,N$, to create the $N\times m_t$ matrix $\boldsymbol{V}_t$. For each $t=1,\ldots,T-1$, define the linear projection for step 1: $$\bm W_t = \bm V_t\Pi_t^{0}+ \bm\eta_t,$$ with $\Pi_t^{0} \stackrel{\mathrm{def}}{=} \bar\mathop{\mbox{\sf E}}(V_{it}V_{it}^{\top})^{-1}\bar\mathop{\mbox{\sf E}}(V_{it}W_{it})$, and $\boldsymbol{\eta}_t$ is an $N\times1$ vector of errors $(\eta_{it})_{i=1}^N$. Assumptions (ref)--(ref) guarantee that each component in $V_{it}$ satisfies the finite moment conditions $\mathop{\mbox{\sf E}}|V_{it,k}|^{2r}<\infty$ and $\mathop{\mbox{\sf E}}|V_{it,k}\eta_{it}|^r<\infty$, for some $r\geq2$, where $k=1,\ldots,m_t$, $i=1,\ldots,N$, $t=1,\ldots,T-1$.
In addition to the population OLS coefficients $\Pi_t^0$, we consider a sparse approximation $\Pi_t^*= \arg\min\limits_{|\Pi_t|_0\leq s_t^*}\bar\mathop{\mbox{\sf E}}|V_{it}^\top(\Pi_t- \Pi_t^0)|^2$, where $s_t^*$ (such that $0<s_t^*\leq m_t$) characterizes the sparsity of $\Pi_t^*$. This sparsity level is determined by the degree of temporal dependency in the data. For instance, using a special case of the panel AR(1) model from Example (ref), we demonstrate that the coefficients in $\Pi_t^0$ exhibit geometric decay governed by the autoregressive coefficient $\theta_1^0$, which leads to the approximately sparse structure of the moment conditions.
We illustrate the sparsity structure in Proposition (ref) with a special case for analytical tractability. In this example, the projection coefficients are zero for regressors beyond $X_{it}=(D_{it},Y_{i,t-1})^\top$ and are therefore exactly sparse. This exact sparsity arises from the serial independence of $D_{it}$. If $D_{it}$ is serially correlated, for example follows an ARMA(1,1) process with innovations $v_{it}$ correlated with $\varepsilon_{is}$ for $s<t$, then the projection coefficients tail off as the time distance increases and become approximately sparse.
To quantify the approximation error of $\Pi_t^*$ with respect to $\Pi_t^0$ empirically, we consider the prediction norm defined by
We shall express all the general rate of convergence results in terms of $C_{s_t^*}$. Under the approximate sparse condition, Theorem (ref) in Appendix (ref) shows the oracle order $s_t^*\asymp \log N \wedge t$, which represents the optimal solution to the risk minimization problem. Consequently, the corresponding oracle bound for the approximation error follows as $C_{s_t^*}\lesssim_\mathrm{P} \sqrt{(\log N \wedge t)/N}$.
Let $\Pi_{t,k}^*$ be the $k$-th element of $\Pi_t^*$, $k=1,\ldots,m_t$. Define the indices sets $J_t\stackrel{\mathrm{def}}{=}\{k\in\{1,\ldots,m_t\}:\Pi_{t,k}^*\neq0\}$ and $J_t^{c}\stackrel{\mathrm{def}}{=}\{k\in\{1,\ldots,m_t\}:\Pi_{t,k}^*=0\}$. For any $\delta_t\in\mathbb{R}^{m_t}$, let $J_{t,0}\subseteq\{1,\ldots,m_t\}$ be a set of indices with cardinality $|J_{t,0}|\leq s_t^*$, and let $J_{t,1}\subseteq\{1,\ldots,m_t\}$ be the set of indices corresponding to the $s_t^*$ largest in absolute value coordinates of $\delta_t$ outside of $J_{t,0}$. In the case of $s_t^*>m_t-|J_{t,0}|$, it corresponds to the $m_t-|J_{t,0}|$ largest absolute values. Let $J_{t,01}\stackrel{\mathrm{def}}{=} J_{t,0}\cup J_{t,1}$. Define $\delta_{t,J_t}$ as the sub-vector of $\delta_t$ corresponding to $J_t$, and define $\delta_{t,J_t^{c}}$ and $\delta_{t,J_{t,01}}$ similarly.
To show identification of $\Pi_t^*$, we consider two events (for each $t$) associated with the restricted eigenvalue (RE) conditions, as outlined in Section 3 of bickel2009simultaneous. For $c_0>0$, define
where $\kappa_t(c_0,s_t^*)$ represents positive constant that depends on $c_0$ and $s_t^*$. In Lemma (ref), we will prove that these events occur with probabilities approaching 1 as $N\to\infty$, for some $\kappa_t(\cdot)>0$ related to the RE condition in population, as per Assumption (ref).
Lemma (ref) shows that $\mathrm{P}(\mathcal A_{1t})\to1$ as $N\to\infty$ for each $t$, which is in line with Lemma 1 of belloni2013least. We can similarly verify that $\mathrm{P}(\mathcal A'_{1t})\to1$ as $N\to\infty$. These results ensure that the Gram matrix $V_{it}V_{it}^\top$ is well conditioned along the sparse directions and establish the identification of the sparse solution $\Pi_t^*$ within the subspace, provided such a solution exists.
Recall the LASSO estimator $\widehat\Pi_t$ obtained by (ref). To achieve good prediction performance of the estimator, properly chosen penalty tuning parameters and weights are necessary. For each $t=1,\ldots,T-1$, let $\bm\omega_t$ be an $m_t\times1$ vector, with the first element being 1 and the remaining elements collecting the non-negative penalty weights $(\omega_{ts}\mathbf{1}_s)_{s=1}^{t}$, where $\mathbf{1}_s$ represents a vector of ones with the same dimension as $X_{is}$.
Assumption (ref) is consistent with Eq. (2.3) in bickel2009simultaneous and serves as a necessary condition for $\Pi_t^0$ to be the minimizer of the LASSO problem expressed in (ref) with probability at least $1-\alpha$. This assumption is crucial for establishing the consistency of our estimator and implies that an ideal choice of the tuning parameter $\lambda_t$ is given by the ($1-\alpha$) quantile of the random variable $c|\bm V_t^\top\bm\eta_t\oslash\bm\omega_t|_\infty$.\footnote{Empirically, since $\bm\eta_t$ is unobserved, the tuning parameter $\lambda_t$ can be selected either based on quantiles of the standard normal distribution or through a more data-dependent approach using the multiplier bootstrap, as discussed in lasso2018. In our empirical analysis, to avoid over-fitting, we adopt a data-independent choice of $\lambda_t$, which is more conservative, and account for heteroskedasticity through the penalty weights $\bm\omega_t$. Further details regarding the practical selection of $\lambda_t$ and $\bm\omega_t$ are provided in Section (ref).} Under Assumptions (ref)--(ref) and (ref), we can apply Lemma (ref), which provides the maximal tail probability for the partial sum of the $m_t$-dimensional process $\varpi_{it}\stackrel{\mathrm{def}}{=} V_{it}\eta_{it}\oslash\bm\omega_t$, to derive an upper bound for the ideal choice of $\lambda_t$, which is of order $\sqrt{N\log m_t}$.
Lastly, to conclude the consistency of the LASSO estimators, we present the prediction performance bounds for $\delta_{\Pi,t}\stackrel{\mathrm{def}}{=}\widehat\Pi_t-\Pi_t^*$ in the following theorem. This will be combined with the prediction norm of the approximation error in (ref) to derive a performance bound for $\widehat\Pi_t-\Pi_{t}^0$ using the triangle inequality.
Based on Theorem (ref), Corollary (ref) provides the joint prediction performance bounds, where the $\ell_2$-norm bound is derived by following Theorem 7.2 of bickel2009simultaneous.
According to the oracle order of $s_t^\ast$ and the corresponding oracle bound for $C_{s_t^*}$, as shown in Appendix (ref), we obtain explicit rates for the general results in Theorem (ref) and Corollary (ref). Specifically, when $s_t^\ast \asymp \log N \wedge t$, $C_{s_t^\ast} \lesssim_\mathrm{P} \sqrt{(\log N \wedge t)/N}$, $\lambda_t \lesssim_\mathrm{P} \sqrt{N\log m_t}$, and there exists a positive constant $\underline\kappa$ such that $\kappa_t(3,s_t^\ast)\ge \underline\kappa>0$, we have $$|\delta_{\Pi,t}|_{2,N} \lesssim_\mathrm{P}\sqrt{(\log N \wedge t)\log m_t/N},\quad |\delta_{\Pi,t}|_{1} \lesssim_\mathrm{P}(\log N \wedge t)\sqrt{\log m_t/N}.$$ Moreover, the joint error bounds satisfy the following rates: $$\max_{1\leq t\leq T-1}|\delta_{\Pi,t}|_{1} \lesssim_\mathrm{P}(\log N \wedge T)\sqrt{\log T/N},\quad \max_{1\leq t\leq T-1}|\delta_{\Pi,t}|_{2} \lesssim_\mathrm{P}\sqrt{(\log N \wedge T)\log T/N}.$$
In this subsection, we establish the asymptotic normality of the final estimator for both AB-LASSO and AB-LASSO-SS, which will allow us to perform large sample inference on the parameters of interest and functions of them.
Define $\bm\Theta_t^0$ (resp. $\widehat{\bm\Theta}_t$) by stacking $\Pi_t^0$ (resp. $\widehat\Pi_t$) by rows for each component $W_{it}$ of $\Delta\widetilde X_{it}$. Specifically, when the number of components in $X_{it}$ is $d$, we have ${\bm\Theta}_t^0$ and $\widehat{{\bm\Theta}}_t$ with dimensions $d\times m_t$ for each $t=1,\ldots T-1$. Recall the definition $V_{it} = (1,X_{i1}^\top,\ldots,X_{it}^\top)^\top$. It follows that $\widehat{\Delta\widetilde X_{it}}=\widehat{\bm\Theta}_tV_{it}$, and the AB-LASSO estimator obtained in (ref) can be expressed by $$\widehat\theta - \theta^0 = \bigg(\sum_{i=1}^N \sum_{t=1}^{T-1} \widehat{\bm\Theta}_tV_{it}\Delta\widetilde X_{it}^{\top} \bigg)^{-1} \bigg(\sum_{i=1}^N \sum_{t=1}^{T-1} \widehat{\bm\Theta}_tV_{it} \Delta\widetilde\varepsilon_{it}\bigg).$$ The asymptotic variance of $\widehat\theta$ has the sandwich form. We impose the following assumption to derive the specific formula for it.
Note that Assumption (ref)(ii) implies $\mathop{\mbox{\sf E}}(\varepsilon_{it}\varepsilon_{is} \mid V_{it})=0$ for $1\leq s<t\leq T$ and $\mathop{\mbox{\sf E}}(\Delta\widetilde\varepsilon_{it}\Delta\widetilde\varepsilon_{i,t-\ell}\mid V_{it})\lesssim\frac{1}{T-t+\ell}\big(1-\frac{1}{N}\big)$ for any $\ell\geq 1$. It follows that $\lim\limits_{T\to\infty}T^{-1}\sum_{t=\ell+1}^{T-1}\mathop{\mbox{\sf E}}(\Delta\widetilde\varepsilon_{it}\Delta\widetilde\varepsilon_{i,t-\ell} \mid V_{it})=0$ for any $\ell\geq 1$ and $$\lim\limits_{N,T\to\infty}(NT)^{-1}\sum\limits_{i=1}^N\sum\limits_{1\leq s<t\leq T-1}\mathop{\mbox{\sf E}}(\Delta\widetilde\varepsilon_{it}\Delta\widetilde\varepsilon_{is} \mid V_{it})=0.$$ Therefore, $\Delta\widetilde\varepsilon_{it}$ can be treated as approximately a martingale difference sequence over $t$. By defining $\Sigma_{0,t}\stackrel{\mathrm{def}}{=}\lim\limits_{N\to\infty}N^{-1}\sum_{i=1}^N\mathop{\mbox{\sf E}}(V_{it}V_{it}^\top(\Delta\widetilde\varepsilon_{it})^2)$, we can express the asymptotic variance of $\widehat\theta$ in the form of
Accordingly, the empirical analog of $\Omega$ is
where $\widehat Q = \{N(T-1)\}^{-1}\sum_{i=1}^N\sum_{t=1}^{T-1}\widehat{\bm\Theta}_tV_{it}\Delta\widetilde X_{it}^{\top}$ and $\widehat \Sigma_{0,t} = N^{-1}\sum_{i=1}^N V_{it}V_{it}^\top(\widehat{\Delta\widetilde\varepsilon_{it}})^2$, with $\widehat{\Delta\widetilde{\varepsilon}_{it}}=\Delta\widetilde Y_{it}-\Delta\widetilde X_{it}^\top\widehat\theta$.
We now establish the formal asymptotic properties of the AB-LASSO and AB-LASSO-SS estimators. The key results are consistency and asymptotic normality, which together justify the use of these estimators for inference on $\theta^0$ in large panels.
A notable feature of Theorem (ref) is that both estimators share the same limiting distribution, despite AB-LASSO-SS employing sample splitting. This reflects the fact that, in our setting, the bias reduction from sample splitting does not affect the central limit theorem, and both estimators attain the same $\sqrt{NT}$-rate and limiting Gaussian distribution. This theoretical finding is consistent with our simulation evidence, where AB-LASSO and AB-LASSO-SS exhibit similar inferential performance in large samples.
We conclude this subsection with a supporting lemma showing that the variance estimator $\widehat{\Omega}$ is consistent, which is essential for constructing feasible confidence intervals in practice.
This lemma requires that the number of relevant instruments $s_t^*$ grows slowly relative to $N$, a mild condition that is standard in high-dimensional IV settings and aligns with the rate condition in Theorem (ref).
We now extend the asymptotic theory to the general model (ref), which allows for many exogenous covariates. We present the limiting distribution of the estimator $\widehat\theta_1$ obtained via (ref).
The key difference relative to the basic model is that the instrument vector now takes the form $U_{it}=(U_{it}^{0\top}, \Delta \widetilde X_{2,it}^\top)^\top$, combining the ideal IVs for the endogenous component $\Delta \widetilde X_{1,it}$\footnote{The ideal IVs for $\Delta\widetilde X_{1,it}$ are structured similarly to those for $\Delta\widetilde X_{it}$ in the basic model, expressed as ${\bm\Theta}_t^0V_{it}$. The covariates $V_{it}$ are expanded by the additional moments arising from the strictly exogenous covariates in $X_{2,it}$, as commented in Remark (ref).} with the exogenous component $\Delta \widetilde X_{2,it}$, which serve as their own instruments. The asymptotic variance of $\widehat\theta_1$ inherits this structure and takes the sandwich form $$ \Omega_1 = Q_1^{-1}\bm\Sigma_1(Q_1^{-1})^\top, \quad \bm\Sigma_1 = \lim_{N,T\to\infty}(NT)^{-1}\sum_{i=1}^N\sum_{t=1}^{T-1}\mathcal W_t^\top\mathop{\mbox{\sf E}}(U_{it}U_{it}^\top(\Delta\widetilde\varepsilon_{it})^2)\mathcal W_t, $$ where $Q_1=\lim\limits_{N,T\to\infty}(NT)^{-1}\sum_{i=1}^N\sum_{t=1}^{T-1} \mathcal W_t^\top\mathop{\mbox{\sf E}}(U_{it}\Delta\widetilde X_{1,it}^{\top})$ captures the relevance of the instruments for the endogenous component and is assumed to be nonsingular, and $\bm\Sigma_1$ is the outer-product variance of the moment conditions, weighted by $\mathcal W_t$. The weighting matrix $\mathcal W_t$ is estimated via a Dantzig selector in practice.
To make inference feasible, we define the moment matrix $M_t\stackrel{\mathrm{def}}{=}\bar\mathop{\sf E}\left\{
U_{it}^{\top}\right\}$ and its empirical counterpart $\widehat M_t\stackrel{\mathrm{def}}{=} N^{-1}\sum_{i=1}^N\left\{
\widehat U_{it}^{\top}\right\}$, where the ideal instruments $U_{it}^0$ in $U_{it}$ are replaced by their LASSO predictions, i.e.\ $\widehat U_{it}=(\widehat{\Delta\widetilde X_{1,it}^\top}, \Delta\widetilde X_{2,it}^\top)^\top$.
The following assumption collects the regularity conditions needed on the weighting matrix $\mathcal W_t$ and the tuning parameters $\ell_t$ under the general model.
Assumption (ref)(i) bounds the magnitude of the weighting matrix and the moment matrix, while Assumption (ref)(ii) translates these bounds into a rate condition on the sparsity of the first-stage LASSO. This condition is analogous in spirit to the rate condition in Theorem (ref), but adjusted for the additional complexity introduced by the high-dimensional exogenous covariates and the estimated weights. Assumption (ref)(iii) restricts the order of the tuning parameter $\ell_t$ used in the Dantzig selector, which controls the penalty level and the deviation $|\widehat{\mathcal W}_t - \mathcal W_t|_{\max}$.
The theorem confirms that the $\sqrt{NT}$-rate and Gaussian limiting distribution established for the basic model carry over intact to the general setting with many exogenous covariates, provided the additional regularity conditions in Assumption (ref) are satisfied. The estimator thus remains suitable for inference in the empirically relevant case. Notably, the crucial rate condition in Assumption (ref)(ii) remains unchanged if a sample-splitting procedure is employed. Specifically, $\widehat{\mathcal W}_t$ obtained via the Dantzig selector and the instrument selection by LASSO are constructed from a sub-sample that is cross-sectionally independent of the sub-sample used to compute the final estimator in (ref). Further details are provided at the end of the proof of Theorem (ref).
We illustrate the finite sample properties of the proposed AB-LASSO and AB-LASSO-SS estimators through numerical simulations, where we also compare them with other alternative methods. Following moral2013likelihood and moral2019dynamic, we use the following data generating process of bun2006effects:
where $\alpha_i\stackrel{\operatorname{i.i.d.}}{\sim}\operatorname{N}(0,2.96)$. We consider 2 cases, one with conditional homoskedasticity and the other with conditional heteroskedasticity. In the homoskedastic case, $\varepsilon_{it}$ and $v_{it}$ are i.i.d.\ $t(4)$ (Student's $t$-distribution with four degrees of freedom), whereas in the heteroskedastic case, $\varepsilon_{it}=\{1+0.5*\mathbf{1}(v_{it}>0)\}e_{it}$, with $e_{it}$ and $v_{it}$ being i.i.d.\ $t(4)$. The sequences ${e_{it}}$, ${v_{it}}$, and ${\alpha_i}$ are mutually independent. Here, $D_{it}$ is predetermined with respect to $\varepsilon_{it}$ but not strictly exogenous, with the parameter $\phi$ capturing the feedback from $Y_{i,t-1}$ to $D_{it}$. In addition, the fixed effects are correlated with both $D_{it}$ and the feedback. We set $\theta_1=0.75$, $\theta_2=0.25$, $\rho=0.5$, $\phi=-0.17$, and $\pi=0.67$. While our methods do not rely on stationarity, we initialize the process with the same starting values of $Y$ and $D$ as moral2013likelihood and moral2019dynamic, which ensure mean stationarity. In Appendix (ref), we conduct another numerical simulation where the data generating process is calibrated to the empirical study in Section (ref).
We assume that $Y_{it}$ is first observed at $t=1$ and use all the available lags of $Y_{it}$ and $D_{it}$ to construct the moment conditions for AB-LASSO and AB-LASSO-SS with FOD, that is, $$ \mathop{\mbox{\sf E}}(Z_{it} \Delta\widetilde{\varepsilon}_{it}) = 0, \quad Z_{it} = (Y_{i,t-1}, \ldots, Y_{i1}, D_{it},\ldots,D_{i1})^{\top}, \quad t = 2,\ldots,T-1, $$ and $\mathop{\mbox{\sf E}}(D_{i1} \Delta\widetilde{\varepsilon}_{i1}) = 0. $ See Remark (ref) for a description of how the AB-LASSO and AB-LASSO-SS are modified when we do not observe the initial condition $Y_{i0}$. The LASSO fitting is carried out using post-LASSO with penalty tuning parameters $\lambda_t$ that are independent of the design matrix. Specifically, for each $t=2,\ldots,T-1$, we set $\lambda_t=c\sqrt{N}\Phi^{-1}(1-0.1/(2m_t))$, where $c=1.1$ and $m_t$ denotes the number of instruments, i.e., the number of regressors in the LASSO regression for each time period $t$.\footnote{For the practical choice of the constant $c$, we recommend calibrating it on the data at hand. In our application-based calibration (Appendix (ref)), we find that $c=1.5$ performs well and adopt this value for the empirical analysis in Section (ref), while in this simulation study we find that $c=1.1$ performs well.} To account for heteroskedasticity, the penalty weights associated with each regressor $V_{it,k}$, for $k=1,\ldots,m_t$, are defined as $\sqrt{\sum_{i=1}^N V_{it,k}^2\hat\eta_{it}^2}\big/\sqrt{N}$, where $\hat\eta_{it}$ are preliminary estimates of the error terms in the first step regressions. This approach is adopted for its conservatism to prevent overfitting lasso2018. For the AB-LASSO-SS estimator, we implement a $K$-fold sample-splitting and cross-fitting procedure with $K \in \{2,5\}.$ According to the formula for the asymptotic variance estimation shown in (ref), we calculate the standard error using the residuals based on the final estimate obtained after aggregation. We also repeat the estimation 100 times using random sample splits and aggregate the results through the medians.
In addition to the conventional two-step AB, we compare the results with several alternatives. First, the debiased AB of chen2019mastering using one random split (i.e. 2 folds) along the cross-section dimension. We refer to this estimator as DAB-SS. Second, the analytically debiased fixed effects estimator (DFE-A), which corrects the bias of the FE estimator arising from the incidental parameter problem nickell1980correcting,kiviet1995bias,hahn2002asymptotically. Third, the likelihood-based estimator moral2013likelihood,moral2019dynamic. For all the comparison estimators, we use analytical standard error clustered at the individual level. We do not consider the half-panel jackknife bias corrected estimator of chudik2018half because it relies on unconditional stationarity of all the variables, which is restrictive for many applications. For example, most of the covariates in our empirical application of Section (ref) are staggered policy indicators that are not stationary over time.
For each estimator, we report the root mean square error (RMSE), standard deviation (SD) and bias in percent of the true parameter value, together with the length and empirical coverage of confidence intervals (CI) with a nominal confidence level of $95\%$. All of them are computed using 500 simulations. Tables (ref) and (ref) display the results for the treatment coefficient $\theta_2$, which is typically of primary interest, for $T \in \{20; 30; 40; 50; 60\}$ in the heteroskedastic case for $N=100$ and $N=200$, respectively. Similar results for the homoskedastic case can be found in Appendix (ref). These sample sizes lead to $m \in \{360; 840; 1,520; 2,400; 3,480\}$ moment conditions for AB, which are large relative to the corresponding sample sizes $n = NT \in \{2,000 ; 3,000; 4,000; 5,000; 6,000\}$ when $N=100$ and $n = NT \in \{4,000 ; 6,000; 8,000; 10,000; 12,000\}$ when $N=200$. Thus, the resulting orders of the small bias condition are $m^2/n \in \{64.8; 235.2; 577.6; 1,152; 2,018.4\}$ when $N=100$ and $m^2/n \in \{32.4; 117.6; 288.8; 576; 1,009.2\}$ when $N=200$, which are not negligible.
We find that AB-LASSO consistently outperforms AB and DAB across all evaluation criteria for all sample sizes. It is evident that using LASSO to select the most relevant moments significantly reduces bias and results in more accurate coverage rates than AB. The CI length of our methods is generally shorter than that for AB, suggesting that the bias reduction does not come at the expense of more dispersion. Comparing with DFE-A, our methods have similar SD, RMSE and CI length, but smaller absolute bias. Moreover, AB-LASSO and AB-LASSO-SS display coverages closer to the nominal level than DFE-A. The undercoverage of DFE-A is due to the bias being comparable to the standard deviation.
When the panel is relatively short, the performance of AB-LASSO in terms of RMSE and CI coverage is comparable to that of ML. Although ML attains a lower RMSE when $T=20$, this advantage vanishes for the other values of $T$ considered. We acknowledge that selecting lags through LASSO might not be essential when the number of instruments is modest (e.g., when $T=20$). However, ML is only feasible when the number of observations is large relative to the number of moment conditions. Specifically, with $N=100$, the maximum feasible value of $T$ is $50$. Moreover, for $N=100$ and $T=50$, the likelihood-based model fails to converge in roughly two-thirds of the simulations.\footnote{We implement the maximum likelihood estimator on dynamic panel models using the dpm function in the R package dpm dpm.} In such cases, parameter estimates are available (as initial values are returned), but standard errors are not, as measures of fit are not reported. Consequently, we do not report ML results for $T\in\{50,60\}$ when $N=100$. Our methods remain robust as $T$ increases, demonstrating the practical suitability for longer panels when $N$ is moderate. In practice, long panels with $T>50$ and moderate $N$ arise in applications such as modeling GDP growth rates for OECD (Organization for Economic Co-operation and Development) countries. For example, galvao2011quantile used annual data from 1948 to 2008 on 18 countries, that is $N=18$ and $T=61$, where our methods remain feasible but ML becomes impractical. In general, relatively long panels are common in macroeconomic applications that use annual or quarterly data on a set of countries or regions (e.g., garcia1987macroeconomic,hoogstrate2000pooling).
In this case sample splitting does not further improve the performance of AB-LASSO significantly, likely because FOD already mitigates much of the overfitting bias and there are only two right-hand-side observed variables in the equation for $Y_{it}$. Overall, the results for AB-LASSO-SS remain robust across different choices of $K$ in terms of estimation, while using $K=5$ split folds improves the inference precision in some cases by producing narrower CI with less conservative coverage rates. Comparing the results for $N=200$ and $T=20$ in Table (ref) with those for $N=100$ and $T=40$ in Table (ref), we find that AB-LASSO and AB-LASSO-SS exhibit superior performance in the latter case, showing lower bias and more accurate CI coverage. This indicates their robustness when the time dimension $T$ is large relative to $N$, while the sample size $n=NT$ remains fixed. Finally, we give an example comparing the complexity of AB and AB-LASSO in terms of computer memory usage. Implementing two-step AB with an efficient weighting matrix under heteroskedasticity for $N=200$ and $T=30$ using the pgmm function in the R package plm plm on a single sample (with $840$ moment conditions) consumes approximately $1.7$ GB of peak memory, whereas AB-LASSO-SS with a single 5-fold partition requires only $45$ MB.\footnote{We track the RAM usage in R using the package peakRAM peakRAM.}
We apply AB-LASSO to study the effect of K-12 schools opening and other policies on the spread of COVID-19 in the U.S. We use a balanced panel of 2,510 counties over 32 weeks between April 1st and December 2nd, 2020. This panel was extracted from chernozhukov2021association, which constructed an unbalanced panel of U.S. counties including 7-day moving averages of daily observations for the same period. We aggregate the observations at the week level to avoid spurious serial correlation coming from the moving averages.
In this application, $Y_{it}$ is the logarithm of the number of reported COVID-19 cases in county $i$ at week $t$, $D_{it}$ is a measure of visits to K-12 schools from SafeGraph foot traffic data, $C_{it}$ contains other treatments and control variables, $\alpha_i$ is a county fixed effect and $\gamma_t$ is a week fixed effect. We estimate the model:
where $C_{1it}$ includes a measure of visits to colleges and policy indicators on mask mandates, stay-at-home orders and the ban on gatherings of more than 50 persons, and $C_{2it}$ includes a measure of the weekly growth rate in the number of tests. We assume that there is no serial correlation in $\varepsilon_{it}$ over $t$. The variables in $D_{it}$ and $C_{1it}$ enter the model lagged one week to account for the time lag between infection and case confirmation. Additionally, we assume that $D_{it}$, $C_{1it}$ and $C_{2it}$ are predetermined with respect to $\varepsilon_{it}$. Accordingly, we use $Y_{i,t-1},\ldots, Y_{i1}$, $D_{it},\ldots, D_{i1}$, $C_{1it},\ldots, C_{1i1}$, and $C_{2it},\ldots, C_{2i1}$ to construct moment conditions at each $t=1,\ldots,T-1$, for AB-LASSO with FOD. This yields $m=3,375$ moment conditions and $n=NT = 2,510 \times 27 = 67,770$ observations.\footnote{The first 5 weeks are used as initial conditions.} AB is likely to be biased in this case because $m^2/n \approx 168$.
Table (ref) presents estimates and standard errors for AB-LASSO, AB-LASSO-SS with $K=2$, AB, DAB-SS, and DFE-A.\footnote{For conciseness, we do not report the results for the autoregressive coefficients, as they are not the policy parameters of interest, but they are available from authors upon request.} The likelihood-based estimator is excluded because applying the dpm function in R that we use in the numerical simulations results in a singular design error. For the AB-LASSO and AB-LASSO-SS estimators, the penalty parameters are chosen based on the calibrated simulation in Appendix (ref). The coefficients of the model measure the short run effects of the covariates. In addition to these coefficients, we report results for the long-run effects obtained as $\theta_k/(1-\sum_{j=1}^4\beta_j)$ where $\theta_k$ is the coefficient of the covariate of interest and $\beta_1,\ldots,\beta_4$ are the coefficients of $Y_{i,t-1},\ldots,Y_{i,t-4},$ respectively.
All the methods reveal positive and significant effects of K-12 school openings and college visits on the spread of COVID-19. In particular, allowing K-12 school openings and college visits is associated with a higher number of cases in both the short and long run. The estimated effects of mask mandates, banning gatherings, and stay-at-home orders are negative and significant both in the short and long run. The positive coefficient of weekly test growth is found to be marginally significant.
Compared with AB, AB-LASSO and AB-LASSO-SS produce considerably smaller absolute estimates of the long-run effects, and smaller estimates for college visits in both the short and long run. Our methods also provide more precise estimates of both effects than AB, particularly for school opening and college visits. Compared with DFE-A, AB-LASSO and AB-LASSO-SS predict smaller absolute short-run effects but larger absolute long-run effects for school opening and college visits, as well as substantially larger absolute effects for long-run impacts of mask mandates, stay-at-home orders and banning gatherings.\footnote{In the calibrated simulation of Appendix (ref), we find that AB-LASSO and DFE-A have biases with reverse signs for the short run effects of K12 school openings, college visits, mask mandates and stay-at-home orders.} In general, the results of AB-LASSO-SS are not sensitive to the number of folds $K$, so we report only $K=2$. Debiasing the standard AB through half-splitting the panel does not substantially change most short-run effect estimates.
According to AB-LASSO, we conclude that the opening of K-12 schools one week is associated with an increase in the number of COVID-19 cases of about $50\%$ the week after and has a compounded long run increase of nearly $300\%$. Allowing college visits one week is associated with approximately an $80\%$ increase in cases in the following week and about a $500\%$ increase in the long run. Mask mandates, stay-at-home orders, and banning gatherings have more modest effects. The reduction in the number of cases is $10\%$, $7\%$, and $7\%$ after one week and roughly $65\%$, $40\%$, and $45\%$ in the long run, respectively. These effects are all statistically and economically relevant.
We propose a LASSO and cross-fitting based estimator of dynamic linear panel models. This estimator shows better large sample properties and finite-sample performance in simulations than the classical AB estimator in long panels. In an empirical application, our estimator finds that policies such as the closure of K-12 schools, mask mandates and stay-at-home orders reduced the spread of COVID-19 both in the short and long run. Our estimates of the long run effects, however, are less optimistic than the estimates obtained with AB.
A potential avenue for future research is to analyze the performance of our method under weak identification. This situation arises, for example, in dynamic models when the process of the outcome $Y_{it}$ is very persistent such that lagged values of $Y_{it}$ might not be strongly correlated with the transformed outcome $\Delta Y_{it}$. Here, we can follow blundell1998initial and employ a system GMM estimator that exploits additional moment conditions for the equation in levels, if we make stationarity assumptions about the initial conditions of the process. The number of these moment conditions grows fast with $T$ and may introduce additional many moment conditions bias. In preliminary numerical simulations, we find that LASSO does not do a good job selecting among these moment conditions as they do not have a natural sparse structure. We leave a more detailed treatment of the system GMM estimator with many moment conditions to future research. Alternatively, newey2009generalized developed a jackknife GMM estimator to deal with many weak instruments for cross-sectional data, which could be extended to incorporate cross-fitting in our panel setting. Here, we also conjecture that standard methods for dealing with weak instruments for cross-section data, such as the use of the Anderson-Rubin statistics, can be fruitfully applied to our setting anderson1949estimation. We also leave this analysis to future research.
\vskip 2em \centerline{ \bf Appendix} \vskip -1em \setcounter{subsection}{0} \vskip 2em