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.
104,338 characters · 10 sections · 78 citation commands
Inference in High-Dimensional Panel Models: Two-Way Dependence and Unobserved Heterogeneity
\setcounter{page}{1} \thispagestyle{empty} \pagestyle{plain}
In economic research, high dimensionality typically refers to the large number of unknown parameters relative to the sample size, under which traditional estimations are either infeasible or tend to yield noisy estimates. The issue of high dimensionality becomes more relevant as data availability grows and economic modeling involves more flexibility. Commonly, the problem of high dimensionality appears at least in three scenarios as follows:
Particularly, the modeling of unobserved heterogeneity in panel models makes high dimensionality more of a practical issue rather than just a theoretical concern. As a motivating example, let's consider a panel model where all three sources of high dimensionality are involved:
where ${D}_{it}$ are low-dimensional treatments or policy variables. ${X}_{it}$ are high-dimensional controls. $D_{it}$ can also contain some higher-order effects and interactive effects with a subset of the controls to allow for nonlinear and heterogeneous effects in a parametric way. $g(.)$ is an unknown function, e.g., an infinite-dimensional parameter; $c_i$ and $d_t$ are unobserved heterogeneous effects. The interest is in the inference on the low-dimensional parameters ${\theta}_0$.
Without considering the features of panel data and the unobserved heterogeneity, it is a classic partial linear model that has been well-studied in the semiparametric literature. With high dimensionality, sparse approximation and regularization approaches have been widely employed to reduce the dimensionality. Essentially, regularization, also known as the machine learning approach, trades off bias for smaller variance to achieve desirable rates of convergence. However, due to the bias introduced by regularization and overfitting, inference can be challenging. Typically, some bias-correction procedures are involved to obtain estimators with better statistical properties and to conduct valid inference.
In the case of panel data, three challenges remain with the existing high-dimensional approaches. First of all, the statistical properties of many regularized estimators remain unknown with panel data where dependence exists across space/unit and time. Secondly, some bias-correction procedures for inference, such as sample-splitting/cross-fitting, are particular about sampling assumptions, and existing approaches are not valid under two-way dependence in the panel. Thirdly, the unobserved individual and time effects may appear in a flexible way, which further complicates estimation and inference.
Although the standard LASSO that assumes the Gaussian error condition does not restrict the dependence structure of data as long as the (long-run) variance is finite, the validity of this approach largely depends on whether a theoretically required penalty level can be found in practice. The common approach for the standard LASSO to decide the penalty level is the K-fold cross-validation, and its validity has been established in chetverikov2021cross, but proper modifications and validity for dependent data remain as open questions. In another strand of literature, gao2024robust establish new Fuk-Nagaev type concentration inequality, by which they show two variants of LASSO continue to work with two-way dependence and non-Gaussianity, with a rate requirement on the penalty level. However, to determine the practical penalty level, it is unclear if the modified BIC approach continues to work under two-way dependence. Furthermore, the validity of this approach under cluster dependence where the correlation may be strong within the cluster remains unknown. chernozhukov2021lasso deal with a system of (time series) equations, and they propose to run LASSO equation-by-equation while determining the penalty level jointly through a multiplier bootstrap which accounts for cross-sectional dependence. However, their performance bound for LASSO only works for each time series and is not compatible with a pooled panel. Also, its validity under cluster dependence is not unknown. I proposed a variant of LASSO that uses regressor-specific penalty weights robust to two-way cluster dependence and weak temporal dependence across clusters. Such a LASSO approach is labeled as the two-way cluster-LASSO, corresponding to the heteroskedasticity-robust LASSO in belloni2012sparse and the cluster-LASSO in belloni2016inference. This approach theoretically derives the common penalty level $\lambda$ up to a constant and a small-order sequence that does not vary across different data-generating processes. Therefore, data-driven tuning, such as information criterion and cross-validation, is not needed, which makes it more computationally efficient and avoids non-trivial theories that take data-driven tuning into account. In Table (ref), a comparison of aforementioned LASSO-type estimators are presented.
According to the rate result, the proposed two-way cluster-LASSO is consistent as both $N,T$ diverge jointly, but the convergence rate is not as fast as the common rates for LASSO under independence or weak dependence\footnote{The intuition for the slow rate of convergence is illustrated in \hyperref[app_0]{Appendix}.}. When $N,T$ diverge at a comparable rate, the rate of convergence is slower than $(N\wedge T)^{-1/2}=(NT)^{-1/4}$ , which is a common rate requirement for inferential theory. This is where the second challenge arises: if a faster rate of convergence is not achievable due to the two-way cluster dependence, some bias-correction approaches are needed to relax the rate requirement for valid inference. There are many bias-correction methods for high-dimensional models. The orthogonalization of moment functions in Chernozhukov2018 provides a general way for constructing estimators that features multiplicative error terms, by which the rate requirement on the nuisance parameter estimation can be relaxed. Combining with a cross-fitting procedure to further control the overfitting bias, they obtain valid inferential results for high-dimensional regression models. However, cross-fitting is sensitive to the sampling assumption. Building upon recent development of cross-fitting approaches for dependent data (Chiang2022,Semenova2023), I propose a clustered-panel cross-fitting scheme and establish its validity. Effectively, this inferential procedure extends the double/debiased machine learning (DML, hereafter) approach by Chernozhukov2018 to panel data models, and it is labeled as the panel DML. Asymptotic normality for the panel DML estimator and the consistency for the variance estimator are established. It is shown that the crude requirement on the rate of convergence can be relaxed to $(N\wedge T)^{-1/4}$, which admits the first-step estimation through the two-way cluster-LASSO.
For the third challenge caused by the unobserved heterogeneity, existing approaches assume that $(c_i,d_t)$ are either additive (belloni2016inference,kock2019uniform,clarke2025double) or interactive (vogt2022cce). To allow for flexible function forms while remaining tractable, I propose to model $(c_i,d_t)$ as correlated random effects through a generalized Mundlak device . In that way, a very rich form of heterogeneity is permitted. A closely related idea has been implemented in wooldridge2020inference and clarke2025double, what's different in this paper is that both unit and temporal unobserved effects are considered, and they are not separable from observable covariates. This is made possible by exploring the sparsity condition on $g_0$. Furthermore, a subtle issue of cross-fitting is discussed with the presence of unobserved heterogeneity. Although valid inference remains challenging for high-dimensional models without cross-fitting in general, I show that inferential theory can be established in model ((ref)) using the full sample with a slightly stronger sparsity condition.
The simulation results are quite revealing. When the model is mildly high-dimensional, non-robust methods, i.e. those designed for i.i.d. data or one-way cluster data, are comparable to the proposed two-way cluster-robust method, and the unregularized method remains valid. This is true for both i.i.d and two-way clustered panel. However, when the nuisance parameters are truly high-dimensional, the pattern changes significantly: (1) with i.i.d data, the unregularized method completely fails the task, and the proposed two-way cluster-robust methods has comparable performance to the non-robust methods in term of both estimation and inference; (2) with two-way clustered panel data, both unregularized methods and non-robust high-dimensional methods severely over-reject in the tests, and the proposed two-way cluster-robust method dominates in terms of estimation bias, SD, RMSE, and, particularly, inference coverage.
In the empirical application, I re-examine the effect of government spending on the output of an open economy following the framework of nakamura2014fiscal, a well-cited empirical macro paper. The baseline model is not concerned with the high dimensionality: a linear panel model with a small number of covariates and additive unobserved heterogeneous effects; the identification is through the instrumental variable. However, even in a conventionally low-dimensional setting, high dimensionality could be hidden because the true model could be highly nonlinear in the covariates and the unobserved heterogeneity. To avoid the endogeneity caused by the potential misspecification in the function form, I consider extending the baseline model in a flexible way as in (ref). The proposed dependence-robust estimation and inference for high-dimensional models can be leveraged, and the results can be used for a robustness check. It is shown that the estimates are consistent with the baseline results, which indicates that the nonlinear and interactive effects may not be very relevant in this model. However, existing approaches that are not robust to high dimensionality or two-way cluster dependence tend to over-fit, bringing noisy estimates and inaccurate inference results.
The rest of the paper is outlined as follows: Section (ref) presents the two-way cluster-LASSO estimator and the investigation of its statistical properties under two-way cluster dependence. Section (ref) introduces the clustered-panel cross-fitting for inference. Section (ref) studies the partial linear model with unobserved heterogeneity as a leading example. Simulation evidence is given in Section (ref). Section (ref) presents an empirical estimation of the government spending multiplier as an illustration of hidden high dimensionality and the application of the proposed toolkit. Section (ref) concludes the paper with empirical recommendations.
Here is a collection of frequently used notations in this paper. Some extra notations are defined along with the context. ${\rm E}$ and ${\rm P}$ are as generic expectation and probability operators. $\mathcal{P}_{NT}$ is an expanding collection of all data-generating processes ${P}$ that satisfy certain conditions. ${P}_{NT}$ is a sequence of probability laws such that ${P}_{NT} \in \mathcal{P}_{NT}$ for each $(N,T)$. The dependence on $(N,T)$ and ${P}_{NT}$ will be suppressed whenever clear in the context. $\Vert . \Vert$ is the Euclidean (Frobenius) norm for a matrix. Let $\mathbf{x}$ be a generic $k\times 1$ real vector, then the $l^q$ norm is denoted as $\Vert \mathbf{x} \Vert_q := \left(\sum_{j=1}^k x_{j}^q \right)^{1/q}$ for $1\leq q < \infty$; $\Vert \mathbf{x} \Vert_{\infty} := \max_{1\leq j\leq k} |x_j|$. The $L^q({P})$ norm is denoted as $\Vert f \Vert_{{P},q} := \left(\int \Vert f(\omega)\Vert^q d {P}(\omega)\right)^{1/q} $ where $f$ is a random element with probability law ${P}$. I denote the empirical average of $f_{it}$ over $i=1,...,N$ and $t=1,...,T$ as $\mathbb{E}_{NT}[f_{it}] = \frac{1}{NT} \sum_{i=1}^N \sum_{t=1}^T f_{it}$ and the empirical $L^2$ norm as $ \left\Vert f_{it} \right\Vert_{NT,2} = \left(\frac{1}{NT} \sum_{i=1}^N\sum_{t=1}^T \Vert f_{it}\Vert ^2\right)^{1/2}$. Correspondingly, I denote the empirical average of $f_{it}$ over the sub-sample $i\in I_k$ and $t\in S_l$ as $\mathbb{E}_{kl}[f_{it}] = \frac{1}{N_k T_l} \sum_{i\in I_k, t\in S_l} f_{it}$ and the empirical $L^2$ norm over the subsample as $ \left\Vert f_{it} \right\Vert_{kl,2} = \left(\frac{1}{N_kT_l} \sum_{i\in N_l}\sum_{t\in T_l} \Vert f_{it}\Vert ^2\right)^{1/2}$, where $I_k,S_l$ are sub-sample index sets and $N_k,T_l$ are sub-sample sizes that will be introduced next section.
In the existing literature, not much is known in terms of statistical properties for high-dimensional methods under cluster dependence in both cross-section and time. In this section, a variant of the $l1$-regularization methods, also known as the LASSO, is proposed and examined. To focus on the LASSO approach under two-way dependence, I consider a simple conditional expectation model of a scalar outcome given a potentially high-dimensional vector of covariates. Let $(Y_{it},X_{it})$ be a sample with $i=1,...,N$ and $t=1,...,T$. The conditional expectation model can be expressed as follows:
where $f(X_{it}):= {\rm E}[Y_{it}|X_{it}]$ is an unknown conditional expectation function of potentially high-dimensional covariates $X_{it}$; $V_{it}$ is the associated stochastic error.
To characterize the two-way cluster dependence in the panel, I assume the random elements $W_{it}:=(Y_{it},X_{it},V_{it})$ are generated by the following process: \begin{assumption_ahk}[Aldous-Hoover-Kallenberg Component Structure Characterization]
where $\mu = {\rm E}[W_{it}]$, $f$ is some unknown measurable function; $(\alpha_{i})_{i\geq 1}$, $(\gamma_{t})_{t\geq 1}$, and $(\varepsilon_{it})_{i\geq 1, t\geq 1}$ are mutually independent sequences, $\alpha_{i}$ is i.i.d across $i$, $\varepsilon_{it}$ is i.i.d across $i$ and $t$. $\{\gamma_t\}_{t\geq 1}$ is strictly stationary and is beta-mixing at a geometric rate:
for some constants $\kappa>0$ and $c_{\kappa}\geq 0$, where $\beta(X,Y) = \frac{1}{2}\Vert {P}_{X,Y}-{P}_{X}\times {P}_{Y} \Vert_{TV}$ and $\Vert . \Vert_{TV}$ denotes the total variation norm. \end{assumption_ahk}
Assumption \hyperref[ahk]{AHK} is motivated by a representation theorem for an exchangeable array, named after Aldous-Hoover-Kallenberg (AHK, hereafter), which states that if an array of random variables $(X_{ij})_{i\geq 1, j\geq 1}$ is separately or jointly exchangeable\footnote{An array $(X_{ij})_{i\geq 1, j\geq 1}$ is separately exchangeable if $\left(X_{\pi(i),\pi'(j)}\right)\overset{d}{=} \left(X_{ij}\right)$, and jointly exchangeable if the same condition holds with $\pi=\pi'$.}, then $X_{ij} = f(\xi_i,\zeta_{j},\iota_{ij})$ where $(\xi_i)_{i\geq 1},(\zeta_j)_{j\geq 1},(\iota_{ij})_{i\geq 1, j\geq 1}$ are mutually independent, uniformly distributed i.i.d. random variables. However, the exchangeability is not likely to hold for arrays with the presence of a temporal dimension since it is naturally ordered. In macroeconomics, for instance, we can interpret the time components $(\gamma_{t})_{t\geq 1}$ as unobserved common time shocks, which are naturally correlated over time, implying that the exchangeability is violated. Therefore, by allowing $\gamma_t$ to be correlated, it introduces temporal dependence across all clusters, making the characterization more sensible in the panel data setting. The beta-mixing condition restricts the temporal dependence of the common time effects to decay at an exponential rate, which is common in literature. The relaxation of the independence condition on $(\gamma_{t})_{t\geq 1}$ can be viewed as a generalization of the component structure representation, as argued by CHS_Restat. It is important to note that the components in (ref) simply characterize the dependence in panel data. Differing from factor models or models with unobserved heterogeneity, they do not affect the identification of the regression model.
Due to the potential high dimensionality in $X$, traditional nonparametric methods are not appropriate for estimating the unknown function $f$ due to the curse of dimensionality. A common approach to reduce the dimensionality is by taking into account the sparsity information in models through regularization. Although the unknown function $f$ is an infinite-dimensional parameter, which is not exactly sparse, we can view it from a sparse approximation perspective, following belloni2012sparse:
\begin{assumption_asm}[Approximate Sparse Model] The unknown function $f$ can be well-approximated by a dictionary of transformations $f_{it} = F(X_{it})$ where $f_{it}$ is a $p\times 1$ vector and $F$ is a measurable map, such that
where the coefficients $\zeta_0$ and the approximation error $r_{it}$ satisfy
\end{assumption_asm} Assumption \hyperref[asm]{ASM} views the high-dimensional linear regression as an approximation. It requires a subset of the parameters $\zeta_0$ to be zero while controlling the size of the approximation error. Compared to the sparsity condition in previous literature, here it imposes a slower rate of growth restriction on the non-zero slope coefficients. For example, $s = o(NT)$ corresponds to the case of heteroskedasticity-robust LASSO under i.i.d data in belloni2012sparse; $s=(Nl_T)$ corresponds to the cluster-robust LASSO under temporal dependence panel data in belloni2016inference where $l_T\in [1,T]$ is an information index that equals T when there is no temporal dependence and equals 1 when there is cross-sectional independence and perfect temporal dependence. In other words, the underlying component structure restricts the growth of nonzero slope coefficients of the model in a way similar to the perfect temporal-dependence case.
Under Assumption \hyperref[asm]{ASM}, we can rewrite the model (ref) as
We then apply $l1$ regularization in the least squared error problem under the model (ref). Let $\lambda$ be some non-negative common penalty level and $\omega$ be some non-negative $p\times p$ diagonal matrix of regressor-specific penalty weights. Consider the following generic weighted LASSO estimator:
We note that the level of penalty term is determined by both $\lambda$ and $\omega$. A large penalty level can reduce avoid noisy estimation due to overfitting but meanwhile too large a penalty level can cause under selection which introduces missing variable bias. This is clearly a tradeoff between overfitting variance and regularization bias. Therefore, to obtain the desirable property of LASSO estimation, $\lambda$ and $\omega$ need to be determined jointly both in theory and in practice. A common choice of $\omega$ is $\mathbb{E}_{NT}[\dot{f}_{it}'\dot{f}_{it}]$ where $\dot{f}_{it}$ be the demeaned $f_{it}$ using the sample mean\footnote{The demeaning is done because of the inclusion of the intercept term which is not penalized.}. With this choice of $\omega$, regressors are standardized, so the model selection is not affected by the scale of the regressors. In theory, given $\omega$, $\lambda$ needs to be chosen in a way that the following event happens with high probability:
where $C_\lambda>1$ is some unknown constant. Condition (ref) is referred to as the “regularization event” in the literature. If the error term $V_{it}$ is conditionally Gaussian or sub-Gaussian, then the Gaussian tail inequality implies an asymptotic order for $\lambda$. To enable applications in broader settings such as asymmetry and heavy tails of the error term, recent literature has investigated in new tools that ensures validity of LASSO with weaker conditions. Particularly, for panel data models, babii2023machine derives a Fuk-Nagaev type concentration inequality for panel data which ensures the validity of LASSO under non-Gaussianity and weak temporal dependence. Under functional dependence measure, {gao2024robust establish new Fuk-Nagaev type concentration inequality and chernozhukov2021lasso extend Gaussian approximation results to ensure LASSO validity with the presence of temporal and cross-sectional dependence. In practice, however, all approaches mentioned above require further estimation of data-driven tuning for the penalty level using cross-validation, information criterion, or bootstrap. These tuning methods are usually computationally costly and hard to justify in theory, and they are further complicated by the temporal and cross-sectional dependence in panel \footnote{chernozhukov2021lasso do provide theoretical justification for choosing $\lambda$ using the multiplier block bootstrap method.}. Furthermore, the restriction imposed for the functional dependence measure in chernozhukov2021lasso and gao2024robust excludes the two-way cluster dependence considered here as we allow the dependence to not decay over cross-sectional or time.
belloni2012sparse view $\omega$ as a self-normalizer and leverage moderate deviation theorems for the self-normalized sums to choose the penalty level that ensures Condition (ref). While this approach does not require extra data tuning once $\omega$ is chosen properly according to the theory, the validity of the approach is restricted: existing moderate deviation theorems only work for independent or weakly dependent random variables. For panel data with cross-sectional independence, one can cluster within each cross-sectional unit or construct temporal blocks that are approximately independent, but there is no existing moderate deviation theorem with self-normalizer that works for two-way dependence. Instead, I utilize the component structure characterization of the dependence and consider a Hoeffding-type decomposition of the high-dimensional mean-zero error term $f_{it}'V_{it}$: $$a_i = {\rm E}[f_{it}'V_{it}|\alpha_i], \ \ g_t = {\rm E}[f_{it}'V_{it}|\gamma_t], \ \ e_{it} = f_{it}'V_{it} - a_i - g_t.$$ The goal is to design penalty weights that account for the randomness from all three components and, importantly, being adaptive for both non-degenerate and degenerate cases.
To illustrate, let's first define a generic three-term regressor-specific penalty weight as follows:
where $\omega_{a,j},\ \omega_{g,j}$, and $\omega_{e,j} $ are non-negative weights correspond to the three components. Let $a_{i,j}$ be the entries of $a_i$ for $j=1,...,p$; $g_{t,j}$ and $e_{it,j}$ are defined similarly. To pin down an appropriate choice of penalty weights, we note that the following inequality holds, for each $j=1,...,p$,
where $c_\omega:= \frac{\omega_{a,j}+\omega_{g,j}+\omega_{e,j}}{\sqrt{\omega_{a,j}^2+\omega_{g,j}^2+\omega_{e,j}^2}}\in \left[1,\sqrt{3}\right]$ is a scaling constant. To see why (ref) holds, we note that for any positive numbers $A,B,C$, $a,b,c$ and $d=\frac{a+b+c}{\sqrt{a^2+b^2+c^2}}$, we have $$\left\{\frac{A+B+C}{\sqrt{a^2+b^2+c^2}}>z\right\} \subseteq \left\{\frac{A}{a}> \frac{z}{d} \right\} \cup \left\{\frac{B}{b}>\frac{z}{d}\right\} \cup \left\{\frac{C}{c}>\frac{z}{d}\right\}.$$ To prove it, suppose $\neg \left\{\left\{\frac{A}{a}>\frac{z}{d}\right\} \cup \left\{\frac{B}{b}>\frac{z}{d}\right\} \cup \left\{\frac{C}{c}>\frac{z}{d}\right\}\right\}$, which implies $\left\{\frac{A}{a}\leq \frac{z}{d}\right\} \cap \left\{\frac{B}{b}\leq \frac{z}{d}\right\} \cap \left\{\frac{C}{c}\leq \frac{z}{d}\right\}$. Then, we have $\frac{A+B+C}{\sqrt{a+b+c}} = \frac{A}{a} \frac{a}{\sqrt{a+b+c}} + \frac{B}{b} \frac{b}{\sqrt{a+b+c}} + \frac{C}{c} \frac{c}{\sqrt{a+b+c}} \leq \frac{z}{d} d =z $, which is $\neg \left\{\frac{A+B+C}{\sqrt{a+b+c}}>z\right\} $, so it is proved by contrapositive.
Inequality (ref) shows that it suffices to consider the tail probability of each component when choosing the penalty level. This greatly simplifies the problem since each of the components possesses much more tractable statistical properties compared to the original error term. It is shown in Appendix that, $a_{i,j}$ is independent over $i$, $\omega_{g,j}$ is weakly dependent over $t$, $e_{it,j}$ is independent conditional on $\{\gamma_t\}$, along with other useful properties. With the observation above, the infeasible regressor-specific penalty weights are proposed as follows:
where $B$ is the number of clusters/blocks, $h$ is the block length and $H_b$ is the associated index set. We note that $\omega_{a,j}$ can also be seen as a sample variance estimator (without centering around the sample mean). Furthermore, $\omega_{e,j}$ can be seen as a cluster variance estimator without estimation error, and $\omega_{g,j}$ can be seen as a correlated-cluster variance estimator (e.g., bester2008inference). It may seem natural to use their sample analogs as the feasible penalty weights, by estimating the components and then plugging in. However, it turns out that a direct sample analog of the infeasible penalty weights is not a good idea in certain settings. To see that, let's first focus on the estimation of the component without considering the estimation error in $V_{it}$, since the latter one is not the main source of the problem. For what follows, we consider the component estimators given $V_{it}$: let $v_{it,j}:= f_{it,j}V_{it}$ and define
Let $\tilde\omega_{a,j}^2$, $\tilde\omega_{g,j}^2$, and $\tilde\omega_{e,j}^2$ be defined as (ref) with the components replaced by those in (ref), and let $\tilde{\omega}_j^2 = \tilde\omega_{a,j}^2 + \tilde\omega_{g,j}^2 + \tilde\omega_{e,j}^2$. We also define the variances $\sigma_{a,j}^2 = E[a_{i,j}^2]$, $\sigma_{g,j}^2 = E[g_{j,t}^2]$, $\sigma_{e,j}^2 = E[e_{it,j}^2]$ and the long-run variances $\Sigma_{g,j} = \sum_{l=-\infty}^\infty E[g_{t,j} g_{t+l,j}]$, $\Sigma_{e,j} = \sum_{l=-\infty}^\infty E[e_{it.j} e_{it+l,j}]$. It can be shown that\footnote{More details can be found in Lemma (ref) in Appendix}, when either $a_{i,j}$ or $g_{i,j}$ is non-degenerate,
In this case, a direct sample analog of the infeasible one in (ref) could be valid. However, when both $a_{i,j}$ and $g_{i,j}$ are degenerate, e.g. $W_{it}$ is i.i.d. over $i$ and $t$, then $T\tilde{\omega}_{a,j}$ and $T\tilde{\omega}_{g,j}$ are exactly zero, and it can be shown that
Here we rescale the penalty weights so that their asymptotic limits are not degenerate. This implies that the rescaled penalty weights would converge to the limit that is at least twice larger as the limit of the infeasible penalty weights, making $\tilde{\omega}_j$ too conservative under the degeneracy. Furthermore, we note that the sample analog of ${\omega}_{e,j}^2$ is mechanically downward biased:
$ \tilde{\omega}_{e,j}^2$ is actually mechanically zero if we recenter the components when constructing the penalty weights. For similar reasons, the finite sample performance of $\tilde{\omega}_{g,j}^2$ would not be satisfying when the size of each correlated cluster $H_b$ is chosen too large.
Due to the aforementioned issues, I propose the following feasible penalty weights: Let $\hat{a}_{i,j}, \hat{g}_{t,j}$, and $\hat{e}_{it,j}$ be defined the same way as (ref) with $v_{it,j}$ replaced by $\hat v_{it,j} = f_{it,j}\hat V_{it}$ and $\hat V_{it}$ is the residual from some initial estimate. Then, we define
where $ k\left(\frac{|t-s|}{M}\right)$ is a Bartlett kernel and $M$ is a bandwidth with a requirement $M=o(T^{1/4})$. In menzel2021bootstrap and hounyo2025, similar adjustments are used to determine presence of the components in order for bootstrapping the components. Here these are used as self-normalizers adaptive to both non-degeneracy and degeneracy scenarios. We also note that $\hat{\omega}_{a,j}^2$ has been used as the feasible penalty weights in the cluster-LASSO method of belloni2016inference but they don't need the extra subtraction adjustment because $\hat{\omega}_{a,j}^2$ itself is in some sense adaptive under one-way cluster dependence, as hinted by Lemma (ref) in Appendix. Additionally, we observe that $\hat{\omega}_{g,j}^2$ is effectively the Driscoll-Kraay variance estimator which accounts for cross-sectional cluster dependence and weak serial dependence, and $\hat{\omega}_{e,j}^2$ is the "average of HACs" estimator that accounts for the serial correlation remaining in $e$.
Since the error term $v$ is unobserved, we need to first obtain the residual $\hat{v}$ from some initial estimation and then iterate the estimation. The detailed implementation is given as Algorithm 1 in \hyperref[app_0]{Appendix}. The validity of this proposal relies on whether the feasible penalty weights under the iterative estimation converges to the same limit of the infeasible penalty weights. We establish the validity formally in Appendix. For the main theorem below, we maintain a high-level assumption on the feasible penalty weights following belloni2012sparse,belloni2016inference: Let $\widehat{\omega}$ be the feasible diagonal weights. Suppose there exists $0<1/c_1<l\leq 1$ and $1\leq u<\infty$ such that $l\to 1$ and
where $\{\omega_j\}$ and $\{\widehat{\omega}_{j}\}$ are diagonal entries of $\omega$ and $\widehat{\omega}$, respectively.
Given this choice of penalty weights above, the common penalty level $\lambda$ is pinned down as
where $C_\lambda = c_1c_\omega>1$; $\gamma$ is a small order sequence. Other than $C_\lambda$ and $\gamma$, there are no unknown tuning parameters in this weighted LASSO estimation. In practice, $C_\lambda$ is chosen as a constant close to 1 and $\gamma = \alpha/\log( N \vee T)$ with $\alpha$ taken as a significance level\footnote{For example, $\alpha=0.1$ is more liberal than $\alpha=0.05$. The choices of $C_\lambda$ and $\gamma$ around the proposed practical choices do not matter much in finite sample.}. The order of $\gamma$ affects the convergence rate of the LASSO estimator: the theory only requires $\gamma=o(1)$ for LASSO to be consistent, but a faster rate of decay of $\gamma$ will result in a slower convergence rate of LASSO.
The key identification condition is given as follows. In the low-dimensional case, the identification in a linear regression is given by the non-singularity of ${\rm E}[f_{it}' f_{it}]$, which implies its empirical counterpart $\mathbb{E}_{NT} [f_{it}'f_{it}]$ is non-singular with high probability. However, as the dimension of $f_{it}$ grows larger than the sample size, $\mathbb{E}_{NT} [f_{it}'f_{it}]$ is singular almost surely. Fortunately, it turns out that under sparsity and $L1$ regularization, we only need certain sub-matrices to be well-behaved for identification. Define
where $\Delta(m)=\{\delta: \Vert\delta\Vert_0 = m, \Vert\delta\Vert_2=1\}$ and $ M_f = \mathbb{E}_{NT} [f_{it}'f_{it}]$.
\begin{assumption_se}[Sparse Eigenvalues] For any $C>0$, there exists constants $0<\kappa_1<\kappa_2<\infty$ such that with probability approaching one, as $(N,T)\to\infty$ jointly, $\kappa_1\leq \phi_{\rm min}(Cs)(M_f)<\phi_{\rm max}(Cs)(M_f) \leq \kappa_2$. \end{assumption_se} The sparse eigenvalue assumption follows from belloni2012sparse. It implies a restricted eigenvalue condition, which represents a modulus of continuity between the prediction norm and the norm of $\delta$ within a restricted set. More primitive sufficient conditions are discussed in bickel2009simultaneous and belloni2012sparse.
\begin{assumption_reg}[Regularity Conditions] (i) $\log (p/\gamma) = o\left(T^{1/6}/(\log T)^2\right)$. (ii) For some $\mu>1,\delta>0$, $\max_{j\leq p}\emph{E}[\lvert f_{it,j} \rvert^{8(\mu+\delta)}]<\infty$, $\emph{E}[\lvert V_{it} \rvert^{8\mu+\delta)}]<\infty$. (iii) $\emph{E}\left[\left(\sum_{t=1}^Te_{it,j}\right)^2|\{\gamma_t\}_{t=1}^T\right] >0$ almost surely. (iv) For each $j$, either (1) $\emph{E}\left(a_{i,j}^2\right) + \emph{E}\left(g_{t,j}^2\right) >\epsilon$ for some $\epsilon>0$ or (2) $f_{it,j}V_{it}$ is i.i.d over $i,t$. (v) $1\le \max_{j\leq p} \omega_j / \min_{j\leq p} \omega_j =O(1) $. \end{assumption_reg}
Assumption \hyperref[regcon]{REG}(i) imposes a restriction on the dimension of $f_{it}$, $p$, while allowing it to be greater than the sample size. The moment conditions in Assumption \hyperref[regcon]{REG}(ii) are common in the literature. \hyperref[regcon]{REG}(iii) is allows for both non-degeneracy and i.i.d case (degeneracy)\footnote{Here the i.i.d is taken as an special case of degeneracy. By Hoeffding decomposition, $e_{it}$ can be further decomposed into a term featuring $E[f_{it}V_{it}|\alpha_i,\gamma_t]$ and a residual term. This is mainly relevant when the both $a_i$ and $g_t$ are degenerate and $E[f_{it}V_{it}|\alpha_i,\gamma_t]$ is not, referred to as the non-Gaussian degenerate case discussed in menzel2021bootstrap. Also see recently davezies2025analytic and hounyo2025. This case still remains as a challenge and is a ongoing topic of research in the literature. We will focus on the Gaussian non-degenerate and degenerate cases.}.
A common way to mitigate the shrinkage bias of LASSO is to apply least square estimation based on the selected model by LASSO, which is named Post-LASSO. We denote the index set of selected regressors as $\widehat{\Gamma} = \{j\in {1,...,p}: |\widehat{\zeta}_j|>0\}$ where $\widehat{\zeta}_j$ are two-way LASSO estimates. The next theorem gives convergence rate results for both two-way cluster-LASSO and its associated Post-LASSO.
Theorem (ref) establishes convergence rates in terms of the prediction, $l1$, and $l2$ norms for the (post) two-way cluster-LASSO estimator in an approximately sparse model. These results are the first that give convergence rates for a LASSO-based estimator allowing for two-way cluster dependence. It is shown that under the two-way cluster dependence, the two-way cluster-LASSO is consistent but, under the non-degenerate case, has a convergence rate slower than those of LASSO-based methods under the random sampling condition or weak dependence. Without loss of generality, let $N = N \wedge T$, then by choosing $\gamma$ according to $\log(1/\gamma) \simeq \log(p\vee N)$, we have $\left\Vert \widehat{\zeta} - \zeta_0 \right\Vert_2= O_P\left(\sqrt{\frac{s\log (p\vee N)}{N} }\right)$. As a comparison, the rate of convergence in terms of the $l2$ norm is $O_P\left(\sqrt{\frac{s\log p }{NT} }\right)$ under the random sampling and the homoskedasticity Gaussian error assumptions in bickel2009simultaneous or the heteroskedasticity Gaussian error in Theorem 19.3 of hansen2022econometrics, $O_P\left(\sqrt{\frac{s\log (p\vee N) }{NT} }\right)$ under random sampling in belloni2012sparse, and $O_P\left(\sqrt{\frac{s\log (p\vee N) }{N I_T} }\right)$ under cross-sectional independence in belloni2016inference where the information index $I_T = 1$ when there is perfect dependence within the cross-sectional cluster.
As illustrated in the Introduction, the slow rate of convergence is due to the underlying factor structure. It is unclear if valid inference is possible under the rate of convergence results in Theorem (ref) or if it is possible to relax the requirement through a cross-fitting procedure. These questions are addressed in the next section.
In this section, I will propose a cross-fitting scheme suitable for clustered panel data. The idea of sample splitting is to split the sample in a proper way and use the sub-samples separately to estimate nuisance parameters and main parameters of interest. If the sub-samples are independent of each other, then the first-step estimation of the nuisance parameters will be independent of the sample used for the second-step estimation. With this property, the error term that causes the bias can vanish with a less stringent rate requirement on the nuisance estimation. Intuitively, the dependence between the two steps is eliminated so that a potentially over-fitted nuisance estimate from the first step does not pollute the second step as much as it would otherwise do.
Therefore, the goal of the cross-fitting scheme is to split the sample in a proper way so that the two resulting sub-samples are independent or, at least, \textquotedblleft approximately\textquotedblright \ independent. Under the AHK characterization in Assumption \hyperref[ahk]{AHK}, $W_{it}$ are cluster-dependent over both cross-section and time. Importantly, the cluster dependence does not vanish as the distance between observations (if there is any ordering) increases. If $\gamma_t$ is weakly dependent, which is the focus of this paper, then the dependence between observations that don't share the same cluster in either dimension dies out as the temporal distance grows. In that case, intuitively, one can split the sample so that the sub-samples do not share the same cluster and are far apart in temporal distance. This is exactly how this scheme works:
Later on, we also use $I_{-k}$ and $S_{-l}$ to denote the index sets for the auxiliary sample $W(-k,-l)$. Similarly, we denote $N_{-k}$ and $T_{-l}$ as the cross-sectional and temporal sample sizes for the auxiliary sample $W(-k,-l)$. Figure 1 illustrates the cross-fitting with $K=4$ and $L=8$.
{ Figure 1: Clustered-Panel cross-fitting with $K=4$ and $L=8$. Three graphs from left to right correspond to the main and auxiliary sample constructions with $(k,l)=(1,1)$, $(k,l)=(2,2)$, $(k,l)=(3,3)$. For a simple illustration, observations in the main sample are all adjacent in the cross-sectional dimension but it is not necessary in practice; the same applies to the auxiliary sample.}
Since the sub-samples $W(k,l)$ and $W(-k,-l)$ do not share any cluster, they are free from cluster dependence, and what's left is the weak dependence over time. Unless imposing $m-$dependence, the sub-samples above will not be independent. However, under certain regularity conditions regarding the weak dependence, it can be shown through the coupling technique that as long as the temporal distance between the sub-samples diverges at a certain rate, there exist coupling sub-samples that are independent of each other while having the same marginal distributions as the constructed sub-samples with probability converging to 1. Such a result is provided in Lemma (ref).
The proof of Lemma (ref) is given in \hyperref[app_0]{Appendix}. Lemma (ref) shows that the main and auxiliary samples from the proposed clustered-panel cross-fitting scheme are approximately independent as $N,T$ diverge. Note that the hypothetical samples $\tilde{W}(k,l)$ and $\tilde{W}(-k,-l)$ only serve as technical tools and do not matter in practice. The proof of Lemma (ref) is based on independence coupling results introduced in Semenova2023.
It turns out such properties of the cross-fitting scheme are very useful for inference in a high-dimensional panel model. Such procedure is used to remove the dependence between the first and the second steps of estimation so as to relax the rate requirement for the first-step estimation of the high-dimensional nuisance parameters. For what follows, I will define an inferential procedure for high-dimensional panel with a two-step estimator and cluster-robust variance estimator using the proposed cross-fitting scheme in the context of a semi-parametric moment restriction model, extending the DML approach in Chernozhukov2018 to a panel setting.
Let $\varphi(W_{it};\theta,\eta)$ be some identifying moment functions where $\theta$ is a low-dimensional vector of parameters of interest and $\eta$ are nuisance functions. For example, $\eta = g_0$ in (ref). Let $\psi(W_{it};\theta,\eta)$ be the orthogonalized moment function with the following properties:
Essentially, $\psi(W_{it};\theta,\eta)$ is adjusted for the fact that $\eta_0$ needs to be estimated and, as a result, the nuisance functions have no first-order effect locally on the orthogonalized moment conditions, based on which the estimation of $\theta_0$ is therefore robust to the plug-in of noisy estimates of $\gamma_0$. In contrast, the original identifying moment conditions do not possess such a property. In model (ref), $\varphi(W_{it};\theta,\eta) = D_{it} U_{it}$ and $\psi(W_{it};\theta,\eta) = \left(D_{it} - {\rm E}[D_{it}|X_{it},c_i,d_t]\right) \left(Y_{it}-D_{it}\theta - g(X_{it},c_i,d_t) \right)$. In the treatment effect model with unconfoundedness conditional on covariates and unobserved heterogeneous effects, $\varphi(W_{it};\theta,\eta) = {\rm E}[Y_{it}|D_{it}=1,X_{it},c_i,d_t]-{\rm E}[Y_{it}|D_{it}=0,X_{it},c_i,d_t] - \theta^{\rm ATE}$ and $\psi(W_{it};\theta,\eta)$ is the moment function corresponding to the well-known augmented inverse probability weighting estimator, which is doubly robust.
Due to the presence of unit and time components, I consider a two-way cluster robust variance estimator similar to CHS_Restat (CHS estimator) with adjustment due to cross-fitting. The variance estimator is motivated under arbitrary dependence in panel data and is shown to be robust to two-way clustering with correlated time effects in linear panel models. chen2024fixed show that such variance estimator can be written as an affine combination of three well-known robust variance estimators: Liang-Zeger-Arellano estimator, Driscoll-Kraay estimator, and the "average of HACs" estimator, and they propose a two-term variant named as, DKA estimator, for better finite sample performance\footnote{chen2024fixed also consider fixed-b asymptotic approximation and bias-correction as a function of the bandwidth $M$. It is not considered in this paper since the statistical property of such adjustment is not clear in this setting.}. We define the CHS and DKA variance estimator with cross-fitting adjustment as follows:
where $ \widehat{\bar{A}}:= \frac{1}{KL} \sum_{k=1}^K\sum_{l=1}^L \frac{1}{N_k T_l} \sum_{i\in I_k,s\in S_l}\psi^a({W_{it};\widehat{\eta}_{kl}})$ and, with $k\left( \frac{m}{M_l}\right):= 1-\frac{m}{M_l}$ for $m=0,1,...,M_l-1$ and 0 otherwise, for some $M_l\in [1,T_l]$,
The next definition summarizes the panel DML estimation and inference procedures for a semiparametric moment restriction model:
A formal study of the estimation and inference procedure given in Algorithm (ref) is given in the Online Supplementary Material. It is shown that the two-step panel-DML estimator of the fixed-dimensional parameter $\theta_0$ is asymptotically normal and the CHS- and DKA-type variance estimators are consistent for the asymptotic variance, under a non-degenerate condition of the underlying components. However, due to the cross-fitting adjustment, there is a subtle issue when the underlying components are degenerate: Since $(K,L)$ are assumed to be fixed, and the sample sizes in each block diverge, the asymptotic approximation of the variance estimators happens in the sub-sample level. If the components that drives the cluster dependence are degenerate, the rate of convergence would change, but the change of scaling for correct asymptotics of the variance estimator may not match with that of the DML estimator. It can be further shown that the proposed variance estimator would be conservative in the i.i.d case: In this case, the asymptotic variance is not degenerate when the DML estimator is scaled by $\sqrt{NT}$. The variance estimator in ((ref)) is designed for the asymptotic variance of the DML estimator scaled by $\sqrt{N}$. Therefore, to be adaptive to the i.i.d case, we need $T\hat{V}_{\rm CHS}$ to be consistent for the asymptotic variance. However, the correct asymptotics requires an adjustment less than $T$, which implies that current specification would be too conservative. This problem is not unique for the method in this paper but rather an intrinsic issue when cross-fitting is implemented in more than one-dimension. For example, Chiang2022 imposes a non-degenerate condition for their main results.
Instead of analytical variance estimators, one could consider resampling method. However, there is one more subtle issue when unobserved heterogeneities are present in the panel model, as revealed in the next section. For the latter issue, it is dealt by a strong assumption on the unobserved heterogeneities or DML without cross-fitting in a partial linear model as a special case. As an implication of these two subtle issues, approaches and sufficient conditions that allow for DML estimation and inference in general settings without cross-fitting would be favorable, and it is the focus of the ongoing research.
In this section, we study the partial linear model with unobserved heterogeneous effects as (ref) as a special case, and an excludable instrumental variable $Z_{it}$ is considered for identification:
For clearer presentation, $D_{it}$ is treated as a scalar variable. If the lags or leads of $D_{it}$ are exogenous, they can also be included in $X_{it}$. Doing so would not change the theory for estimation and inference, but could change the interpretation of $\theta_0$.
To deal with the unobserved heterogeneous effects that cause endogeneity, I take a correlated random-effects approach through the generalized Mundlak device: \begin{assumption_gmd}[Generalized Mundlak Device] For each $i=1,...,N$ and $t=1,...,T$,
where $\bar{F}_i = \frac{1}{T}\sum_{t=1}^T F_{it}$, $\bar{F}_{t} = \frac{1}{N}\sum_{i=1}^N F_{it}$, $F_{it}:= (D_{it},X_{it}')'$; $h_c$ and $h_d$ are some unknown measurable functions; $(\epsilon^c_i,\epsilon^d_t)$ are independent shocks; and $(c_i,d_t)$ are independent of $U_{it}$. \end{assumption_gmd}
To justify this assumption, we shall recall the idea of the conventional Mundlak device. To explicitly model the correlation between the random effects and the covariates, mundlak1978pooling proposes an auxiliary regression between the random effects and the cross-sectional sample average, and shows that if the random effects enter the model linearly then the resulting GLS estimator is equivalent to the common within-estimator. wooldridge2021two further shows that the equivalence relations exist among the POLS estimators resulting from the Mundlak device, within-transformation, and the fixed-effects dummies. Therefore, if aforementioned approaches are sensible for dealing with unobserved heterogeneity, then allowing the Mundlak device to have a more flexible function form should also be reasonable and more robust. Compared to similar ideas implemented in wooldridge2020inference and clarke2025double, the approach here is more flexible, while, as a tradeoff, relying more on the sparsity condition on the nuisance functions.
Before utilizing the proposed panel DML approach for inference, we note that there is a subtle issue: the Mundlak device uses the full history of the covariates which potentially generates dependence across the cross-fitting sub-samples. Similar issues also appear in a simple linear panel model with additive unobserved effects where within-transformation also introduces sample-averages. Therefore, the cross-fitting may not be compatible with the presence of unobserved heterogeneity. Without cross-fitting, it is challenging to establish an inferential theory with growing dimensionality in unknown parameters in general. Nevertheless, as is shown below, it is possible to establishing the asymptotic normality of the panel DML estimator using the full sample with a strengthened sparsity condition.
Under model (ref), $g(X_{it},c_i,d_t) = {\rm E}[Y_{it} - D_{it}\theta_0 | X_{it},c_i,d_t]$. We can rewrite (ref) as follows:
where $g_D(X_{it},c_i,d_t):= {\rm E}[ D_{it} | X_{it},c_i,d_t] $ and $g_Y(X_{it},c_i,d_t):= {\rm E}[ Y_{it} | X_{it},c_i,d_t] $. Under Assumption \hyperref[crecon]{GMD}, $g_D(X_{it},c_i,d_t)$ and $g_Y(X_{it},c_i,d_t)$ can be rewritten as compound functions, which are assumed to be well-approximated by a linear combination of a $\tau-$th order polynomial transformation $L^{\tau}$ as follows:
where $(\eta_D,\eta_Y)$ are slope coefficients and $(r_{it}^D,r_{it}^Y)$ are the approximation errors. Furthermore, we can define a vector of transformed regressors as $L_{1,it} =L^{\tau}(X_{it},\bar{F}_{i}, \bar{F}_{t})$ and a vector of unobserved regressors as $L_{2,it} = L^{\tau}(X_{it},\bar{F}_{i}, \bar{F}_{t},\epsilon_i^c,\epsilon_t^d) \backslash L^{\tau}(X_{it},\bar{F}_{i}, \bar{F}_{t}) $. Let $(\eta_{D,1},\eta_{D,2})$ be such that
And $(\eta_{Y,1},\eta_{Y,2})$ are defined in the same way. Under the sparse approximation and Assumption \hyperref[crecon]{GMD}, we can rewrite model (ref) as follows:
By defining a new error term $V_{it}^g := \left( L_{2,it} -{\rm E} [L_{2,it}]\right)\left(\eta_{Y,2}-\eta_{D,2}\theta_0\right)+ U_{it}$, a new approximation error $r_{it} = r_{it}^Y + r_{it}^D\theta_0$, the vector of observables $f_{it}:= \left(L_{1,it},1\right)$ with dimension denoted by $p$, and the nuisance vectors $\beta_0 := \left(\eta_{Y,1}, {\rm E} [L_{2,it}]\eta_{Y,2}\right)$, $\pi_0 := \left(\eta_{D,1}, {\rm E} [L_{2,it}]\eta_{D,2}\right)$, we can rewrite the model above as
Noticeably, in this case, the parameters associated with the unobservables $L_{2,it}$ can be arbitrarily non-sparse.
Given ${\rm E}[Z_{it}U_{it}]$ and the independence between $Z_{it}$ and $(\epsilon_i^c,\epsilon_t^d)$, we have the identifying moment condition ${\rm E}[Z_{it}V_{it}^g] = 0$. Let $\zeta_0$ be the linear projection parameter of $Z_{it}$ onto $f_{it}$ and let $V^Z_{it}$ be the corresponding linear projection errors. By Eq. (2.18) of Chernozhukov2018, the near-Neyman orthogonal moment function is given by:
where we denote $\eta_0=(\zeta_0,\beta_0,\pi_0)$. Under the sparse approximation, we can also rewrite the conditional expectation models for $Y$ and $D$ as
where ${V}^Y_{it} = \left( L_{2,it} -{\rm E} [L_{2,it}]\right)\eta_{Y,2}+ U_{it}^Y$ and ${V}^D_{it} = \left( L_{2,it} -{\rm E} [L_{2,it}]\right)\eta_{D,2}+ U_{it}^D$. For $l=Z,Y,D$, let ${\omega}_l$ be the infeasible penalty weights for the two-way cluster-LASSO estimation of $(\zeta_0,\beta_0,\pi_0)$, as defined in (ref) with $V_{it}$ replaced by $V^l_{it}$. Correspondingly, let $\widehat{V}^l$ be the residuals, and $\widehat{\omega}_l$ be the feasible penalty weights. The two-step debiased estimator $\widehat{\theta}$ for $\theta_0$ using the full-sample is defined as the solution of $\mathbb{E}_{NT}[\psi_{it}(\theta_,\widehat{\eta})]=0$ where $\widehat{\eta}$ are the (post) two-way cluster-LASSO estimators for $\eta_0$ obtained in the first step using the full-sample.
For statistical analysis, the following notations are used: $ {a}_i = {\rm E}[ V_{it}^ZV_{it}^g|\alpha_i]$, ${g}_t = {\rm E}[V_{it}^ZV_{it}^g|\gamma_t]$, ${\Sigma}_a = {\rm E}[{a}_i{a}_i']$, ${\Sigma}_g = \sum_{l=-\infty}^{\infty}{ E}[{g}_t{g}_{t+l}']$; ${a}_{i,j,l} = {\rm E}[f_{it,j}V^l_{it}|\alpha_i]$, ${g}_{t,j,l} = {\rm E}[f_{it,j}V^l_{it}|\gamma_t]$, $e_{it,j,l} =f_{it,j}V^l_{it}-{a}_{i,j,l} -{g}_{t,j,l}$ for $l =Z,Y,D$; ${A}_0 = { E}[V_{it}^Z{V}^D_{it}]$, ${\Omega}_0 = {\Sigma}_a+ c{\Sigma}_g$.
\begin{assumption_plm}[Regularity Conditions for the Partial Linear Model] \ \ \ \ \ \ \ \ \ \ \
\end{assumption_plm}
This set of regularity conditions follow from the assumptions for two-way cluster-LASSO. Here we focus on the non-degenerate case. The only extra condition is Assumption \hyperref[regcon_app]{REG-P}(ii) which is a smoothness condition that ensures the exogeneity properties of $\bar{F}_i$ and $\bar{F}_t$ inherited from $(c_i,\epsilon_i)$ and $(d_t,\epsilon_t)$.
Theorem (ref) establishes the validity of the proposed inference procedure using the full sample. Note that the sparsity condition and the condition of the approximation errors are stronger than the ones needed for two-way LASSO estimation itself. To estimate the asymptotic variance, the following variance estimators are adapted from CHS_Restat and chen2024fixed using the full sample:
where $ \widetilde{A}_{NT}:= \frac{1}{NT} \sum_{i=1}^N\sum_{t=1}^T (Z_{it}-f_{it}'\widetilde{\zeta})(D_{it}-f_{it}'\widetilde{\pi})$ and
For simplicity, we deliver the consistency results of variance estimators assuming the approximation is exact. Allowing for approximation errors does not change the main idea but only requires more regularity conditions on the approximation error and lengthier derivations.
Theorems (ref) and (ref) together justify the panel DML inference without cross-fitting by exploiting the partial linear structure and the sparsity condition \footnote{For more general models, recent literature resort to extra restrictions to avoid cross-fitting(e.g., chen2022debiased and cao2025neighborhood).}.
In this section, we examine the finite sample performance of proposed two-way cluster-LASSO estimation and debiased inference with or without cross-fitting. To focus on the main challenges caused by two-way clustering and high-dimensional nuisance estimation, the simulation study is anchored on linear approximately sparse models without unobserved heterogeneities. In the previous version of working paper, simulation results for linear model with exact sparsity and partial linear model with nuisance functions are also reported. For partial linear model, the performance is largely determined by whether the nuisance functions are well-approximated by basis functions. For the case of exact sparsity, the performance of LASSO-based method is sensitive to the size of the slope coefficients: a larger slope coefficient leads to better selection in the first stage across all LASSO methods. Since it is hard to justify proper choices of nuisance functions and sizes of slope coefficients in the DGP, we will consider a linear model with the sequence of slope coefficients decay in a polynomial rate, and LASSO-type approaches decides between “signal" and “noise” based on particular penalty choices\footnote{See Shen2025 for a critique on LASSO-type approaches for learning weak signals. They consider a scenario where signals are individually weak but jointly significant, which is implicitly excluded by the approximate sparsity condition considered in this paper. }.
The linear model with high-dimensional covariates, approximate sparsity, and two-way clustering dependence is specified as follows:
where $\theta_0=1$ is the true parameter of interest. $\beta_0 =\pi_0= \left(1, \frac{1}{2^2},..., \frac{1}{p^2}\right)'$ are $p$-dimensional nuisance parameters. $\alpha_{i}^u, \alpha_{i}^v, \varepsilon_{it}^u, \varepsilon_{it}^v$ are each random draws from $N(0,1)$; $(\alpha_{i,1},...,\alpha_{i,p})'$ and $(\varepsilon_{it,1},...,\varepsilon_{it,p})'$ are each random draws from a joint normal distribution with mean zero and variance-covariance matrix equal to $0.5^{|j-k|}$ in the $(j,k)$'s entry; $\gamma_t^u,\gamma_t^v, \gamma_{t,1},...,\gamma_{t,p}$ each follows an AR(1) process with the coefficient equal to $0.5$ and the initial values randomly drawn from $N(0,0.75)$. The weights $(w_1,w_2,w_3) = (1/3, 1/3, 1/3)$ control the relative importance of the components.
The simulation study examines the Monte Carlo bias (Bias), standard deviation (SD), root mean square error (RMSE), and coverage probability of estimators for $\theta_0$. All estimations are based on the orthogonal moment condition given by (ref) with $Z_{it} = D_{it}$ and $f_{it} = X_{it}$. The comparison will be among procedures with and without cross-fitting. The first-step estimations will be based on the POLS estimator (if feasible), the post heteroskedasticity-robust LASSO from belloni2012sparse, the post cluster-robust LASSO from belloni2016inference, and the post two-way cluster-LASSO. The CHS-type and DKA-type variance estimators, given in Section (ref) with cross-fitting and Section (ref) without cross-fitting, are used for inference. \footnote{In some unreported simulations, I also compare CHS/DKA type variance estimators with Eicker-Huber-White type estimators in Chernozhukov2018 for random sampling data and Cameron-Galbach-Miller type estimator from Chiang2022 for multiway clustered data. Since it is well-known that inference based on variance estimators that do not sufficiently account for the dependence would cause over-rejection, it is omitted here}.
The simulation results are based on 1000 Monte Carlo replications\footnote{For replication, the code is included in the supplementary material and also publicly available at \nolinkurl{http://kaichengchen.github.io/twlasso_paneldml_replication.zip}.}. It is a relatively small number of replications but it is necessitated by the high computational cost of multiple high-dimensional estimation and inference procedures, particularly with cross-fitting. For variance estimation, bandwidth parameters $M$ of the Bartlett kernel are required. I use the min-MSE rule from Andrews1991 for both purposes. For a generic scalar score $v_{it}$, the formula is given as $\hat{M}=1.8171\left( \frac{\widehat{\rho}^{2}}{\left( 1-\widehat{\rho }^{2}\right) ^{2}}\right) ^{1/3}T^{1/3} + 1$ where $\widehat{\rho}$ is the OLS estimator from the regression $\bar {\hat{v}}_{t}=\rho\bar{\hat{v}}_{t-1}+\eta_{t}$ where $\bar{\hat{v}}_{t}=\frac{1}{N}\sum_{i=1}^{N}\hat{v}_{it}$ and $\hat{v}_{it} = \hat{U}_{it}\hat{V}_{it}$.
Table \hyperref[table1]{5.1} presents a set of baseline results with a mild dimension of covariates, $p=200$. The number of covariates is much larger than either cross-sectional or temporal dimensions. In the first step, model selections are done using different LASSO approaches reported in the second column. Comparing the results obtained without using cross-fitting, it is shown that when the number of regressors is not extremely large relative to the sample size, the POLS estimator dominates the sparse methods through different LASSOs in terms of Monte Carlo bias, standard deviation, and coverage probability obtained using DKA standard error, even though the true model is approximately sparse. Among the sparse methods, cluster-robust methods including the proposed two-way cluster-LASSO exhibit the smallest biases and best coverages, with slightly larger standard deviation. In terms of selection, the proposed method selects the number of regressors closest to the true number of relevant regressors while other sparse methods over-select to different extents.
When cross-fitting is employed, all methods have witnessed a significant improvement in terms of sample coverage. This is particularly true for LASSO-based method not robust for dependent data. This is not too surprising because those non-robust sparse methods tend to over-select, and the cross-fitting is designed to remove the overfitting bias and to restore asymptotic normality. As a cost of cross-fitting, the Monte Carlo standard deviation increased, indicating the efficiency loss due to the exclusion of sub-samples in the first-step estimation. It is also worth emphasizing that the CHS- and DKA-type variance estimators designed for cross-fitting approaches play an important role in the desirable sample coverage. In some unreported simulations, it is shown that inference based on the cross-fitting variance estimators proposed in Chernozhukov2018 and Chiang2022 suffers from severe under-coverage. It implies that while two-way dependence potentially affects both estimation and inference, its negative impact on the inference is more salient.
As the dimension of the covariates significantly increases and becomes as large as the overall sample size, a different pattern is revealed. Table \hyperref[table2]{5.2} also reports simulation results under the same DGP except that the dimension $p$ now increases to 600, slightly smaller than the overall sample size 625. First, we compare the results obtained without cross-fitting. The simulation results demonstrate that the methods based on the POLS with no selection and those based on the existing LASSO approaches with over-selection all suffer from severe under-coverage. The proposed method, in contrast, remain among the best in terms of bias, RMSE, and coverage. When cross-fitting is performed, there is again a significant improvement across all approaches in terms of the sample coverage but it is also in a cost of a slight efficiency loss measured by the increase in SD.
Finally, we examine the performance of the two-way cluster-robust method as well as the one-way cluster-robust and non-robust methods when the true DGP features i.i.d data. The theoretical results tell that the high-dimensional methods should remain valid. The pattern in Table (ref) matches with the theories, and the non-robust high-dimensional methods performs slightly better than the cluster-robust ones. The conservativeness of the DKA variance estimator is well-known in the cluster-robust inference literature: two-way cluster-robust variance without the double-counting adjustment can ensure the semi-positive definiteness, but when the components are degenerate, it has a probability limit twice of the true asymptotic variance (see, for example, mackinnon2021wild,chen2024fixed).
In this section, I re-examine the effects of government spending on the output of an open economy following the framework of nakamura2014fiscal. It is one of the most cited empirical macro papers in the American Economic Review, and it investigates one classic quantity of interest in economics: the government spending multiplier. The question here is, can we improve the estimation and inference through more robust and flexible methods? As I will show, it is made possible by the proposed toolkit in this paper.
This framework utilizes the regional variation in military spending in the US to estimate the percentage increase in output that results from the increase of government spending by 1 percent of GDP, i.e. government spending multiplier. It is referred to as the "open economy relative multiplier" because this framework takes advantage of uniform monetary and tax policies across the regions in the US to difference-out their effects on government spending and output. The parameter of interest is a scalar, and the baseline model is identified without considering control variables, so why is the high dimensionality relevant here? As discussed below, the high dimensionality from heterogeneity and flexible modeling could be hidden.
Due to the endogeneity in the variation of the regional military procurement, nakamura2014fiscal achieves identification through an instrumental variable (IV) approach. As argued by the authors, the national military spending is largely determined by geopolitical events so it is likely exogenous to the unobserved factors of regional military spending and it affects the regional military spending disproportionally. In other words, the identifying assumption is that the buildups and drawdowns in national military spending are not due to unbalanced military development across regions. Based on this observation, a share-shift type IV is considered and the share is estimated by regressing the regional military spending on the national military spending allowing for region-specific constant slope coefficients.\footnote{All quantities, unless specifically defined, are in terms of two-year growth rate of the real per capita values. Per capita is in terms of total population. nakamura2014fiscal also presents results when per capita is calculated using the working age population as a robustness check.} To focus on the main idea, the shares are taken as given and the resulting instrument variable is treated as observable instead of generated regressors.
In this paper, to avoid the endogeneity caused by the misspecification of the function form, I extend the linear model with additive unobserved heterogeneous effects to a partial linear model with non-additive unobserved heterogeneous effects. Let $D_{it}$ be the percentage change in per capita regional military spending in state $i$ and time $t$ and $Z_{it}$ be the IV. Specifically, the baseline model from the original study and the one from this paper differ as follows:
where $\theta_0$ is the parameter of interest, i.e. the true multiplier; $X_{it}$ and $W_t$ are exogenous control variables with the latter being only time-varying; $\pi_i$ are non-random unit specific slope coefficients of $W_t$; $(c_i,d_t)$ are unobserved heterogeneous effects. In the original study, the linear model is estimated by the two-stage least square (2SLS) with two-way fixed effects. We apply the approach in Section (ref) for the extended model.
In the baseline specification of nakamura2014fiscal, $W_t$ are not included. In their alternative specifications, $W_t$ is chosen as the real interest rate or the change in national oil price. These two variables are never included together in the original study. Note that allowing the unit-specific slope coefficients for controls generates many nuisance parameters: with 51 state groups\footnote{The regions in this analysis are defined by the states. nakamura2014fiscal also presents results on regions as clusters of states.}, one control would increase 51 parameters and two controls would generate 102 parameters, without considering interactions or higher order terms. With a sample size of less than 2000, the high dimensionality in nuisance parameters could result in a noisy estimate of $\theta_0$. In this paper, I consider additional controls. As is shown in Table 3 of nakamura2014fiscal, the change in state population is likely not affected by the treatment (the regional military spending), so it is not subject to the "bad control" issue; But it could affect the treatment and the outcome, so it is included in $X_{it}$. By considering more flexible function forms and additional exogenous control variables, the excludability condition of the instruments is more plausible. On the other hand, the high-dimensionality arose from the flexible function form and the unobserved heterogeneity necessitates the use of high-dimensional methods. Moreover, state-level yearly variables of those macroeconomic characteristics are often considered to be cluster-dependent in both cross-sectional and temporal groups due to common time shocks and state unobserved effects. These concerns justify the use of the proposed robust methods in this paper.
The data is available through nakamura2014fiscal. It is a balanced (after trimming) state-level yearly panel data with 51 states from 1971-2005 years. The military spending data is collected from the electronic database of DD-350 military procurement forms of the US Department of Defense. The state output is measured by state DGP collected from the US Bureau of Economics Analysis (BEA). The state population data is from the Census Bureau. Data on oil prices is from West Texas Intermediate. The Federal Funds rate is from the FRED database of the St. Louis Federal Reserve. The state inflation measures are constructed from several sources. For more details on data construction, readers are referred to nakamura2014fiscal. \footnote{ For replication, the data and code used for this empirical study are included in the supplementary material and also publicly available at \nolinkurl{http://kaichengchen.github.io/twlasso_paneldml_replication.zip}.}
Table \href{table1_emp}{6.1} provides benchmark results for the original model with different choices of control variables. All estimates (columns 6) of are given by 2SLS with two-way fixed effects and the standard errors (s.e.) are calculated using CHS and DKA formulas given in Section (ref). The estimates of the multiplier replicate those given in nakamura2014fiscal with significant differences in the standard errors. It is because the variance estimates here account for the potential two-way dependence while the variance estimator used in nakamura2014fiscal assumes cross-sectional independence.
The main comparisons are done in Tables \href{table2_emp}{6.2} and \href{table3_emp}{6.3}. In Table \href{table2_emp}{6.2}, no cross-fitting is performed in the first stage. The number of parameters associated with regressors generated by the polynomials transformations are reported in column (4). Overall, with more controls and the polynomial transformation of the observables, the standard errors are generally larger than those in \href{table1_emp}{6.1}. With no transformations of the original regressors, the estimates obtained by four different methods are similar and are consistent with the baseline results. As the flexibility increases with the higher-order polynomial transformations, the number of selected regressors increases across all methods. While the standard errors of most approaches become larger and the estimates deviate from the baseline results, the proposed approach remains less noisy. This indicates that many higher-order polynomials included in the extended model for robustness in the function form may not matter much. While the existing approaches tend to over-select those terms under potential two-way dependence, the proposed LASSO method is robust against over-selection and the panel DML estimator remains accurate. Table \href{table3_emp}{6.3} demonstrates the comparison between various sparse methods with the clustered-panel cross-fitting \footnote{Due to a smaller sample used in the first-step estimation and multicollinearity among the polynomial terms, methods based on the POLS first-step is too noisy and so they are omitted for comparison here.}. It reveals a similar pattern as in Table \href{table2_emp}{6.2}: The variability of different methods increases as the model approximated by higher-order polynomial series, except for the two-way cluster-robust approach.
To conclude, the empirical study of the government spending multiplier using a flexible model and sparse methods illustrates the issue of hidden dimensionality. In the current example, the estimates obtained through the high-dimensional methods do not deviate much from the baseline results, so it implies the nonlinear effects omitted from the original model may not be very relevant. While the two-way cluster-LASSO and the inference procedure remain relatively accurate and provide results as a robustness check, other sparse methods tend to over-select and the DML estimators become too noisy to interpret.
The inferential theory for high-dimensional models is particularly relevant in panel data settings where the modeling of unobserved heterogeneity commonly leads to high-dimensional nuisance parameters. This paper enriches the toolbox of researchers in dealing with high-dimensional panel models. Particularly, I propose a package of tools that deal with the estimation and inference in high-dimensional panel models that feature two-way cluster dependence and unobserved heterogeneity. I first develop a weighted LASSO approach for two-way clustered panels. As is shown in the statistical analysis of the two-way cluster-LASSO, the convergence rates are slow due to the cluster dependence, making it challenging for inference purposes. However, by utilizing a cross-fitting method designed for a two-way clustered panel, the rate requirement for the first step can be substantially relaxed, making the proposed two-way cluster-LASSO a feasible first-step estimator for the panel-DML inference procedure in a high-dimensional semiparametric model. Individually, both the two-way cluster-LASSO and the clustered-panel cross-fitting can be of independent interest; Together, they extend the DML approach to panel data settings. Two concerns are also revealed regarding cross-fitting for DML in panel data setting, and alternative approaches are discussed and left for future research. In the partial linear panel model with high dimensionality as a special case, I further consider the unobserved heterogeneity and inferential theory using the full sample. The validity of the full-sample estimation and inference is established under a slightly stronger sparsity condition compared to the general case with cross-fitting.
The estimation and inferential theory are empirically relevant. I illustrate the proposed approaches in an empirical example and exemplify that high-dimensionality could be hidden in questions not traditionally considered high-dimensional. In practice, when the question is naturally high-dimensional and answered by panel data, then the proposed approaches are natural solutions. When the questions are originally not high-dimensional, it is reasonable to start with a simple model as a baseline and then extend it to a more general and flexible model for a robustness check.
The author confirms that the data supporting the findings of this study are available within the article and its supplementary materials.
For helpful comments and discussions, I thank Tim Vogelsang, Hugo Freeman, Jeff Wooldridge, Antonio Galvao, Kyoo il Kim, Shlomo Levental, Chiang Harold, Rahul Singh, Le Wang, Louise Laage, and Saera Oh. I wish to thank Whitney Newey for his enlightening lectures in 2023 Asian Summer School in Econometrics and Statistics and helpful advice. I also appreciate all helpful comments from seminar participants at Michigan State and conference participants at AMES2024, CES2024, ESIF-AIML2024, MEG2024, CES2025, ESWC2025.