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.
146,428 characters · 26 sections · 92 citation commands
Inference after discretizing time-varying unobserved heterogeneity
{\it Keywords:} Unobserved heterogeneity, k-means clustering, panel data, double machine learning, inference
\setstretch{1.5}
Accounting for unobserved heterogeneity is often critical for credible identification in both reduced-form and structural economic analyses. Among the various available dimension reduction devices, economists have increasingly relied on clustering techniques. This strategy has proven particularly effective in panel data models with time-varying unobserved heterogeneity. Early contributions developed valid inference procedures under the assumption of well-separated groups bonhomme2015grouped. Recent work considers more realistic settings with continuous heterogeneity, using clustering only to approximate the unobserved structure through discretization bonhomme2022discretizing, freeman2023linear. This approach has recently gained traction in empirical studies bonhomme2019distributional,jolivet2024structural,mahler2024lifestyle. In these settings, however, while consistent two-step estimators and their rates of convergence are available, valid inference procedures are still lacking, so that practitioners are unable to assess the statistical significance of their results.\footnote{bonhomme2022discretizing does have inference results when the unobserved heterogeneity is time-invariant, but lacks them in the case where it is time-varying, which we study here.}
This paper takes a first step toward filling this gap by establishing the asymptotic normality and unbiasedness, at the parametric rate, of a novel estimator for the slope coefficient in a linear semiparametric version of bonhomme2022discretizing's model. Specifically, for units $i=1,\dots,N$ and dates $t=1,\dots,T$, we assume that there exist unobserved fixed effects $\alpha_i \in \mathbb{R}^{K_\alpha}$ and $\gamma_t \in \mathbb{R}^{K_\gamma}$ such that
where $y_{it} \in \mathbb{R}$ is an observed outcome variable, $x_{it} = (x_{it1}, \ldots, x_{itK})^\top \in \mathbb{R}^K$ is a vector of observed covariates, ${\mathbb{E}}[x_{it}v_{it}]=0$, and $f(\alpha_i,\gamma_t)={\mathbb{E}}[y_{it}-x_{it}^\top \beta|\alpha_i,\gamma_t]$ and $h_k(\alpha_i,\gamma_t)= {\mathbb{E}}[x_{itk}|\alpha_i,\gamma_t],\ k \in \{1, \dots, K\}$, are unknown deterministic mappings from $\mathbb{R}^{K_\alpha} \times \mathbb{R}^{K_\gamma}$ to $\mathbb{R}$. The main features of the model are that the mappings $f$ and $(h_k)_{k\in \{1,\dots,K\}}$ are smooth, the fixed effects are low-dimensional, and the error terms $v_{it}\in{\mathbb{R}}$ and $u_{it} = (u_{it1}, \ldots, u_{itK})^\top\in{\mathbb{R}}^K$ are uncorrelated and sufficiently weakly dependent across $i$ and $t$.\footnote{The existence of such fixed effects can be motivated by the literature on exchangeable arrays. If $(y_{it}, x_{it}^\top)^\top,\ i = 1, \ldots, N,\ t = 1, \ldots, T$, are exchangeable arrays, then the Aldous--Hoover--Kallenberg representation theorems guarantee that $\alpha_i$ and $\gamma_t$ exist and the residuals $v_{it}$ and $u_{it}$ are i.i.d. However, the low-dimensionality of the fixed effects, the smoothness of the functions $f$ and $(h_k)_{k \in \{1, \ldots, K\}}$, and the absence of correlation between $u_{it}$ and $v_{it}$ are not implied by the Aldous--Hoover--Kallenbeng representation.} Our interest lies in the unknown regression parameter $\beta\in{\mathbb{R}}^K$. Importantly, some fixed effects may appear in one equation only.\footnote{For example, suppose that $K=1$, $y_{it} = x_{it1} \beta_1 + f^y(\alpha_i^y, \gamma_t^y) + v_{it} $, and $x_{it1} = h^x(\alpha_i^x, \gamma_t^x) + u_{it1}.$ This can be written as (ref) with $\alpha_i=((\alpha_i^y)^\top,(\alpha_i^x)^\top)^\top$, $\gamma_t=((\gamma_t^y)^\top, (\gamma_t^x)^\top)^\top$, $f(\alpha_i,\gamma_t)= f^y(\alpha_i^y, \gamma_t^y)$, and $h_1(\alpha_i,\gamma_t)= h^x(\alpha_i^x, \gamma_t^x).$ } Moreoever, the fixed effects' contribution to the outcome, $f(\alpha_i, \gamma_t)$, can be flexibly correlated with the covariates through $(h_k(\alpha_i, \gamma_t))_{k \in \{1, \dots, K\}}$. In applications, $\alpha_i$ typically represents consumers' preferences, workers' abilities, or states' political structures, while $\gamma_t$ captures macroeconomic shocks and business cycles. We consider asymptotic regimes such that $N$ and $T$ grow to infinity while $K$ is fixed.
We propose a novel two-step estimation procedure that combines well-established methods from the existing literature. The first step constructs a discrete approximation of the unobserved heterogeneity using k-means clustering of unit-specific and time-specific informative moments (e.g., cross-section or time-series average of the data). The second step is a linear regression with additively separable two-way grouped fixed effects specific to each cluster, estimated by ordinary least squares (OLS). The main contribution of this paper is to formally establish that this simple combination can enjoy parametric-rate asymptotic normality and unbiasedness because it leverages bias-reducing Neyman-orthogonal moments. Such moments are standard in the double machine learning literature chernozhukov2018double and, in the present context, mitigate the influence of estimation errors originating from the clustering step.\footnote{Robust moments, though used in the interactive fixed effects literature pesaran2006estimation,westerlund2015cross,beyhum2023factor,freeman2022linear, are novel in the present non-interactive fixed effects setting.} As a by-product, standard OLS inference routines can be applied to the second-step regression.
We establish the asymptotic normality of a cross-fitted version of the proposed estimator under the condition $\max(N, T) = o(\min(N, T)^3)$, up to logarithmic factors. Cross-fitting is a resampling technique from the double machine learning literature that simplifies the theoretical analysis. We evaluate the finite-sample performance of the proposed estimator and its cross-fitted version by means of Monte Carlo simulations. The baseline estimator exhibits excellent finite sample properties and slightly outperforms its cross-fitted variant, which itself significantly improves upon benchmark estimators. Notably, the confidence intervals from both estimators achieve nearly nominal coverage even when $T$ is much smaller than $N$, a common feature in many microeconomic datasets. We apply the methodology to assess the fiscal response of U.S. states to resource revenues and find results closely aligned with predictions from economic theory. The proposed estimators are implemented in the R package \href{https://github.com/Wei-M-Wei/pcluster}{pcluster}.
Inference procedures have been proposed for models with a fixed effects structure corresponding to special cases of ours. A first related strand of literature is that on panel data models with interactive fixed effects. pesaran2006estimation, greenaway2012asymptotic, and westerlund2015cross assume that both $f$ and $(h_k)_{k \in \{1, \ldots, K\}}$ are interactive and therefore known. bai2009panel imposes that $f$ is known and interactive but does not model covariates. A second related strand of literature is that of the special case of grouped fixed effects models bonhomme2015grouped,chetverikov2022spectralpostspectralestimatorsgrouped,mugnier2024simple, in which an exact group structure is assumed. Following bonhomme2022discretizing, we do not assume that the fixed effects follow a group pattern but instead use clustering as an approximation device.
Several recent papers consider models with nonseparable fixed effects (this structure is also called a “nonlinear factor model" in the literature). Most closely related are freeman2023linear and bonhomme2022discretizing. freeman2023linear study the same outcome model as us, leaving the relationship between the covariates and the fixed effects unrestricted. bonhomme2022discretizing consider nonlinear versions of both the outcome and covariate models with a parametric likelihood specification of the distribution of $y_{it}$ given $x_{it},\alpha_i,\gamma_t$. Both freeman2023linear and bonhomme2022discretizing derive convergence rates but do not establish asymptotic normality, thus falling short of providing the inference tools we develop here. Our novel estimation procedure shares similarities with these two papers. The first-step clustering is the same as that of bonhomme2022discretizing, but unlike the latter paper, we rely on additively separable two-way grouped fixed effects in the second step. freeman2023linear employ the same second step as we do, but the first steps differ between the two papers.
Next, feng2020causal, deaner2025inferring, and athey2025identification consider estimation of the average treatment effect on the treated of a binary treatment when the potential outcome in the absence of treatment follows a nonlinear factor model. In contrast to our proposal, these approaches allow for heterogeneous treatment effects but cannot be applied when the covariate of interest is continuous and do not rely on clustering.\footnote{Relatedly, it has been shown that synthetic control methods are valid under a linear factor model with a growing number of factors, which can be seen as the approximation of a nonlinear factor model, see arkhangelsky2021synthetic,arkhangelsky2023large.} Finally, zeleneev2020identification proposes estimators for linear and nonlinear network models with nonseparable fixed effects and obtains rates of convergence. Unlike in our work, regressors with an exact two-way structure are allowed in the latter paper.
The rest of the paper is organized as follows. Section (ref) introduces the two-step estimation procedure and discusses its link with double machine learning. We then describe the cross-fitted variant of the estimator and provide our main large sample result in Section (ref). Section (ref) displays the results of the Monte Carlo simulations. The application to fiscal policy is developed in Section (ref). Section (ref) concludes. Several key lemmas, all proofs, and results from additional Monte Carlo simulations are presented in the Appendix.
We begin by providing some intuition for the estimation strategy. By plugging the covariate model into the outcome model, we have
where $h_{K+1}(\alpha_i,\gamma_t)= f(\alpha_i,\gamma_t) + \sum_{k=1}^K \beta_k h_k(\alpha_i,\gamma_t)$ and $e_{it}=u_{it}^\top \beta+ v_{it}$. Next, note that ${\mathbb{E}}[v_{it}|\alpha_i,\gamma_t]= {\mathbb{E}}[y_{it}-x_{it}^\top \beta-f(\alpha_i,\gamma_t)|\alpha_i,\gamma_t]=0,$ because $f(\alpha_i,\gamma_t)={\mathbb{E}}[y_{it}-x_{it}^\top \beta|\alpha_i,\gamma_t]$. This, the fact that $x_{it}=h(\alpha_i,\gamma_t)+u_{it}$, and the contemporaneous exogeneity assumption ${\mathbb{E}}[x_{it}v_{it}]=0$, then imply $$ {\mathbb{E}}[u_{itk}v_{it}]=0, \quad k=1,\ldots,K, $$ so that $\beta$ is the slope coefficient in the linear regression of $e_{it}$ on $u_{it}$. Since $e_{it}$ and $u_{it}$ are both unobserved, this linear regression is infeasible but it suggests the following two-step estimation procedure for $\beta$: (i) estimate $e_{it}$ and $u_{it}$, and (ii) linearly regress the estimates of $e_{it}$ on those of $u_{it}$.
To estimate $e_{it}$ and $u_{it}$, we start by constructing a discrete approximation of unobserved heterogeneity across units and dates.\footnote{An alternative approach would be to discretize solely across one dimension (either units or dates). However, as noted in bonhomme2022discretizing and freeman2023linear, this leads to slower rates of convergence. See also beyhum2023factor for a similar argument in panel data models with interactive fixed effects. Simulations in Section (ref) confirm that discretizing along a single dimension yields much worse performance.} Following bonhomme2022discretizing, we focus on the popular k-means clustering algorithm applied to cross-sectional and time-series averages of the data. This approach can be expected to perform well if such averages are informative about the underlying unobserved heterogeneity in a way that can be exploited by the discretization method (see Section (ref) below). Other algorithms are discussed below. Estimates of $(e_{it},u_{it}^\top)^\top$ are then obtained from the residuals of the linear projection of $y_{it}$ and $x_{it}$ on cluster-specific additively separable two-way fixed effects. The two main steps of the proposed estimation procedure are formally described below. Let $z_{it}:=(x_{it}^\top,y_{it})^\top$. \paragraph{Step 1 (Two-way clustering).} Let $G$ and $C$ denote the number of unit and time clusters, respectively (a rule to select them is outlined below). Let $\|\cdot\|$ denote the Euclidean norm.
Clustering algorithm for units. Let $a_i:=\frac{1}{T} \sum_{t=1}^T z_{it}, \ i\in\{1,\dots,N\}$. Compute $$\left(\widehat{a}(1),\dots,\widehat{a}(G),g_1,\dots,g_N\right)\in\operatorname*{arg\,min}\limits_{
}\sum_{i=1}^N \left\| a_i-a(\tilde g_i)\right\|^2.$$
Clustering algorithm for dates. Let $b_t:=\frac{1}{N} \sum_{i=1}^N z_{it},\ t\in\{1,\dots,T\}$. Compute $$\left(\widehat{b}(1),\dots,\widehat{b}(C),c_1,\dots,c_T\right)\in\operatorname*{arg\,min}\limits_{
}\sum_{t=1}^T \left\| b_t-b(\tilde c_t)\right\|^2. $$ The procedures deliver unit and time cluster labels $g_1,\dots,g_N$ and $c_1,\dots,c_T$, respectively. Since some fixed effects can enter the outcome model but not the covariate model, or vice versa, it is crucial to include both $y_{it}$ and $x_{it}$ as inputs of each clustering algorithm. Else, such fixed effects would not be accounted for. Fast computational routines exist to find exact solutions to both k-means clustering problems for data sets of moderate sizes \citep[e.g.,][]{Merle1997AnIP,Aloise2009AnIC}, and local minima for others (e.g., Hartigan--Wong’s algorithm).\footnote{If the quality of the local minima raises suspicion, we recommend using hierarchical clustering approaches as outlined in Appendix~\ref{subsec.hiech_app} as a sensitivity analysis, though we leave the verification of their approximation properties for further research.} \paragraph{Step 2 (Two-way grouped fixed effect estimator).} The estimators of $e_{it}$ and $u_{it}$ are
where, for any variable $w_{it}$, we define
with $N_{g_i}:= \sum_{j=1}^N \mathbf{1}\{g_j=g_i\}$ and $T_{c_t} :=\sum_{s=1}^T \mathbf{1}\{c_s=c_t\}$. These estimators correspond to within-group transformations applied to $y_{it}$ and $x_{it}$ in a similar fashion to the standard within transformations in standard linear panel data models with two-way fixed effects. The final estimator of $\beta$ is the ordinary least squares estimator of $\widehat e_{it}$ on $\widehat u_{it}$,
which is numerically equivalent to the two-way grouped fixed effects regression coefficient $$\operatorname*{arg\,min}_{\beta\in{\mathbb{R}}^K}\min_{\delta\in{\mathbb{R}}^{N\times C}} \min_{\nu\in{\mathbb{R}}^{G\times T}} \sum_{i=1}^N\sum_{t=1}^T \left(y_{it}-x_{it}^\top\beta -\delta_{i,c_t}-\nu_{g_i,t}\right)^2.$$
In contrast, bonhomme2022discretizing considers estimators with either only unit cluster fixed effects of the form $\delta_{i,c_t}$, or interacted unit and time clusters of the form $ \xi_{g_i,c_t}$. The additively separable grouped fixed effects structure that we use here delivers better rates of convergence, see freeman2023linear for a heuristic discussion. The grouped fixed effects estimator in freeman2023linear uses a similar second step but a different first-step clustering procedure; essentially, their proposal uses clusters that approximate $f(\alpha_i,\gamma_t)$ but not $h(\alpha_i,\gamma_t)$.\footnote{The comparison between our estimator and that of freeman2023linear is similar to the relation between, respectively, the Double Lasso and the Post-Lasso estimators, see chernozhukov2024applied for a discussion. The post-Lasso estimator is not asymptotically normal since it only approximates the best linear predictor in the outcome equation. Approximating the best linear predictor in both the outcome and covariate equations, as done by the double Lasso, is key for inference.} In Section (ref), we compare our approach with these alternative estimators in simulations.
Since the two-way grouped fixed effects estimator relies on linear regression, usual standard errors (with a degree of freedom correction) can be used. We note that extending the approach to accommodate a model with unit- or time-heterogeneous slopes ($\beta_i$ or $\beta_t$) is relatively straightforward.
\paragraph{Choice of the number of clusters.} To choose the number of clusters $G$ and $C$, we use the data-driven selection procedure developed by bonhomme2022discretizing. Let $Q_g(G):=\frac{1}{N}\sum_{i=1}^N \left\| a_i-a(g_i)\right\|^2$ and $Q_c(C):=\frac{1}{T}\sum_{t=1}^T \left\| b_t-b( c_t)\right\|^2$ denote the k-means objective functions evaluated at their maxima. The quantities $Q_g(G)$ and $Q_c(C)$ measure the approximation errors made through the clustering. Let $ \widehat{V}_g:=\frac{1}{NT^2}\sum_{i=1}^N \sum_{t=1}^T\left\|z_{it}- a_i\right\|^2$ and $ \widehat{V}_c:=\frac{1}{N^2T}\sum_{t=1}^T\sum_{i=1}^N \left\|z_{it}- b_t\right\|^2$ denote empirical dispersions, which measure the fundamental noise level in the inputs of the clustering procedures. The data-driven choice of the number of clusters is $\widehat{G}:=\min_{G\ge 1}\{G:\ Q_g(G)\le \widehat{V}_g\}$ and $\widehat{C}:=\min_{C\ge 1}\{C:\ Q_c(C)\le \widehat{V}_c\}$. It aims at balancing the approximation error and the input noise. We provide some theoretical guarantees in Section (ref).
\paragraph{On the clustering algorithm.} As in bonhomme2022discretizing, the baseline approach clusters on cross-section and time-series averages using a k-means algorithm. Intuitively, this procedure requires that these averages be informative about the fixed effects. This leads to an “injectivity” condition, formalized in Assumption (ref) below, which imposes that the limit of the averages is injective in the fixed effects. Such an assumption can be relaxed or avoided. One solution is to use moments beyond averages, leading to weaker restrictions. Another approach, studied in Appendix (ref), uses hierarchical clustering on the pseudo-distance of zhang2017estimating, avoiding averaging the data before clustering. We focus on k-means clustering of averages in the main text because of its simplicity and excellent performance in simulations.
Two-step estimation procedures whose second steps are based on Neyman-orthogonal moments lie at the heart of the double machine learning literature chernozhukov2018double. Such moments are bias-reducing because they limit the influence of the errors in estimating the nuisance parameters in the first step and, therefore, make inference possible. This robustness property arises because the difference between the empirical counterpart of the Neyman-orthogonal moment and the infeasible empirical moment based on the true values of the nuisance parameters decomposes into sums of either products of estimation errors or products of an estimation error and an error term; see in particular the discussion in Section 1 of chernozhukov2018double. It turns out that the moment on which our second-step estimator is based exhibits the same type of robustness properties. To see this, note that the second-step estimator solves the empirical moment equation
Moment (ref) approximates the empirical moment equation
solved by an infeasible “oracle” OLS estimator knowing $u_{it}$ and $e_{it}$. Notice that
where
Hence, the difference between the moments (ref) and (ref) is the sum of a term $a^*$, corresponding to the sum of the products of two estimation errors, and two terms $b^*$ and $c^*$ which are sums of products of an estimation error and an error term. All of these terms are, therefore, sums of products of “small terms" and will thus be asymptotically negligible. This explains why the proposed estimator can be asymptotically normal.
In this section, we provide theoretical guarantees for a cross-fitted variant of the estimator. In Section (ref), we motivate and discuss the use of cross-fitting. Section (ref) introduces the cross-fitted version of the two-step estimator. Section (ref) provides sufficient conditions for its asymptotic normality. Section (ref) formally presents the large sample result. Section (ref) contains some results regarding the data-driven choice of the number of clusters.
Deriving the limiting distribution of the least-squares estimator $\widehat{\beta}$ is challenging, as it requires controlling the dependence between the clusters estimated in the first step and the error terms of the data used in the second step. This difficulty is a common feature of many two-step estimators based on highly nonlinear black-box first-step estimators.\footnote{In particular, without a control of the dependence between the two steps, one cannot use concentration arguments on $u_{it}$ and $v_{it}$ to bound the terms $b^*$ and $c^*$ introduced in Section (ref).}
This type of issue has also been encountered in the literature on double machine learning chernozhukov2018double. The solution taken in this research area is to use cross-fitting. The data is split into different folds, and the first-step and second-step estimations are performed on different folds. The role of the folds is then reversed, and the second-step estimators over the different folds are averaged to improve efficiency. Under independent observations, this mechanically eliminates the dependence between the first-step estimator and the data used in the second step, therefore solving the aforementioned problem.
In this section, we follow this strategy to establish the asymptotic normality at the parametric $\sqrt{NT}$-rate of a cross-fitted version of the estimator that learns clusters and estimates the slope coefficient from separate batches of the data.
We emphasize that cross-fitting is merely a proof device, and we recommend using $\widehat{\beta}$ in practice. Indeed, Monte Carlo simulations in Section (ref) demonstrate that the original estimator $\widehat{\beta}$ outperforms its cross-fitted version, which already performs very well. Cross-fitting has been shown not to improve estimator performance in simulations across various settings dukes2021inference, chen2022debiased, vansteelandt2024assumption, wang2024doubly, shi2024off. Moreover, it has been demonstrated that cross-fitting is not always essential for achieving asymptotic results in double machine learning when the learners adhere to a natural leave-one-out stability property chen2022debiased or the lasso is used chernozhukov2015post. These findings suggest that in certain contexts, cross-fitting is not only unnecessary but may even be counterproductive. Our simulation results indicate that k-means clustering is one such learner where cross-fitting can be omitted without compromising performance.\footnote{An intuition of why cross-fitting is not needed in practice is as follows. In our case, the clusters are fully determined by cross-section and time series averages of the variables. A particular observation should be only very weakly dependent on these averages (unless the variables have heavy tails) so that the clustering step does not overfit. In contrast, if one replaces k-means on averages by hierarchical clustering on a pseudo-distance as studied in Appendix (ref), the clusters do not depend only on the averages, and we see in simulations that the estimator without cross-fitting does not perform well. }
To describe the alternative estimator based on cross-fitting, let us consider a simple cross-fitting scheme with only four folds:
We also use the notation $N_d:=|\mathcal{N}_d|$ and $T_d:=|\mathcal{T}_d|$ and note that $\mathcal N_1=\mathcal N_2$, $\mathcal N_3=\mathcal N_4$, $\mathcal T_1=\mathcal T_2$, and $\mathcal T_3=\mathcal T_4$. This type of division in four folds is appropriate for panel data and also appears in freeman2023linear.\footnote{In unreported simulations, we have not found any improvement resulting from increasing the number of folds.}
We briefly outline the construction of the cross-fitted estimator, denoted $\widehat{\beta}^{\rm CF}$. For a detailed presentation, we refer to Appendix (ref). For an observation $(i,t) \in \mathcal{O}_d$, we estimate $u_{it}$ and $e_{it}$ as follows:
The final estimator $\widehat{\beta}^{\rm CF}$ is the pooled OLS estimator from regressing $\widehat{e}_{it}^d$ on $\widehat{u}_{it}^d$. This procedure determines cluster memberships using data distinct from that used for within-transformations, thereby simplifying the theoretical analysis while maintaining efficiency across the entire dataset.
Consider the following assumptions.
Assumption (ref) is a mild regularity condition on $(h_k)_{k\in\{1,\dots,K+1\}}$. Assumption (ref) is similar to Assumption 2 in bonhomme2022discretizing. It is best understood in the case of pointwise limits, where $\text{plim}_{T\to \infty}a_i^d=\varphi_d^\alpha(\alpha_i)$ and $\text{plim}_{N\to \infty}b_t^d=\varphi_d^\gamma(\gamma_t)$, which can be justified by laws of large numbers. Assumption (ref) then requires that the probability limits are injective and imposes some rate of convergence of the sample averages to these limits.
Let us first discuss the injectivity property. It requires that units (resp. time periods) with similar values of time-series (resp. cross-sectional) averages of $z_{it}$ have similar values of unit-specific (resp. time-specific) fixed effects and vice versa, with equality in the limit. Intuitively, such an injectivity property suggests that matching on observed panel data averages is sufficient to control for unobserved heterogeneity (i.e., matching on the fixed effects). It is also useful to analyze the injectivity assumption in an example. Consider the case where $K=d_\alpha=d_\gamma=1$, $\beta=0$ and $f(\alpha_i,\gamma_t)=h_1(\alpha_i,\gamma_t)=\alpha_i\gamma_t$. Then, under weak regularity conditions, $\varphi_d^\alpha(\alpha_i) =(\alpha_i,\alpha_i)^\top {\mathbb{E}}[\gamma_t]$ and injectivity fails to hold only if ${\mathbb{E}}[\gamma_t]=0$, showing that failure of injectivity is the exception rather than the norm in this setting.\footnote{Also note that this assumption is related to the full rank condition in the common correlated effects literature pesaran2006estimation, which guarantees that cross-section averages allow for the recovery of the factors.} As noted earlier, the injectivity property can be avoided by using different clustering approaches such as hierarchical clustering applied on a pseudo-distance matrix, as we study in Appendix (ref).
Next, in Assumption (ref), the rate of convergence of $a_i^d$ and $b_i^d$ to their probability limits in sup-norm is controlled by the sequences $r_\alpha$ and $r_\gamma$. Concentration inequalities boucheron2013concentration can be used to show that the bounds hold for particular values of $r_\alpha$ and $r_\gamma$ under different dependence settings and conditions on the tails of the distribution of $z_{it}$. For instance, we show in Lemma (ref) in Appendix (ref) that if, conditional on $\alpha_i$, $(z_{it})_{t\in\mathcal{T}_d}$ are independent sub-Gaussian random variables with with common mean ${\mathbb{E}}[z_{it}|\alpha_i]$ and sub-Gaussian norm bounded uniformly in $t$ and the value of $\alpha_i$, then the bound on $\max_{i\in\mathcal N_d}\left\|a_i^d-\varphi_d^\alpha(\alpha_i)\right\|^2$ in Assumption (ref)(ref) holds with $r_\alpha=\log(N)$. Under analogous conditions, the bound in Assumption (ref)(ref) holds with $r_\gamma=\log(T)$. As a result, $r_\alpha$ and $r_\gamma$ will typically be negligible with respect to $N$ and $T$, respectively.
The following assumption collects standard dependence, moment, and non-collinearity conditions that prove helpful in establishing the limiting distribution of the estimator. Let $\mathcal F_{NT}$ denote the sigma-algebra generated by $\{\alpha_i,\gamma_t:(i,t)\in\left\{1\ldots,N\right\}\times\left\{1,\ldots,T\right\}\}$.
Assumption (ref)(ref) rules out conditional cross-section or time-series dependence in the error terms, and requires errors to have zero conditional mean, i.e., that they are are mean-independent of the fixed effects, a standard assumption in the panel data literature. It implies that the data from the different folds are independent conditional on the fixed effects. Though it may be arguably strong, relaxing it would require obtaining a precise control of the dependence between the clustering algorithm's outcome and the error terms, which, as noted earlier, is particularly challenging with black-box methods such as k-means. In the simulations reported in Section (ref), we find that the estimator still performs very well under time series correlation. Note that the assumption of i.i.d. errors is commonly made in papers studying sophisticated panel data models; see, for instance, moon2015linear, chen2021quantile, bonhomme2022discretizing, and freeman2023linear. Similar to our work, these papers derive their main theoretical results under this assumption but provide simulation evidence suggesting that the restriction may not be necessary.
Assumption (ref)(ref) requires the idiosyncratic component of each equation to admit slightly more than an uniformly bounded conditional second moment across units, time periods, and regressors. This is useful to verify a Lindeberg--Feller condition and apply a central limit theorem to the dominant term in the estimator.
The first part of Assumption (ref)(ref) is a standard asymptotic non-collinearity condition on the covariates in the second-step regression. Together with the second part of Assumption (ref)(ref), it ensures that the estimator possesses a non-degenerate limiting distribution.
The following assumption specifies the relative rates at which $N$, $T$, and the numbers of clusters $G_d$ and $C_d$ can grow.
As noted above, under the standard conditions of Lemma (ref), $r_\alpha=\log(N)$ and $r_\gamma =\log(T)$ and Assumption (ref)(ref) becomes $\max(N,T)=o(\min(N,T)^3)$ up to logarithmic terms. The latter is weaker than the rate conditions on $N$ and $T$ typically found in the literature on panel data models with interactive fixed effects. For instance, bai2009panel imposes $\max(N,T)=o(\min(N,T)^2)$ to derive asymptotic normality, while the estimators in westerlund2015cross are biased as $T/N$ goes to a constant. This improvement is substantial, as it is obtained while relaxing the modeling assumption that $g$ and $h_k$ are interactive. We relax the condition in bai2009panel thanks to the use of the orthogonal moment stemming from the covariate equations, while we improve on westerlund2015cross by estimating both unit and time-specific fixed effects in the first step, while westerlund2015cross only estimate the factors (corresponding to the time-specific fixed effects in an interactive fixed effects model); see also Footnote (ref) for a related discussion. In contrast, the rate condition (ref) is stronger than that for grouped fixed effects models such as in bonhomme2015grouped, where $T$ can grow at an arbitrary polynomial rate with respect to $N$. This is because we do not assume that the data has an exact group structure, and instead use clustering as an approximation device.
Assumption (ref)(ref) stipulates that both the number of unit clusters and time clusters must be negligible with respect to $N$ and $T$, respectively. Intuitively, this is necessary because, otherwise, the within transformations applied to the data to estimate $e_{it}$ and $u_{it}$ would create non-negligible time series and cross-section dependence in the generated regressors of the second step, precluding the estimator from being $\sqrt{NT}$-consistent.
The last assumption concerns the approximation error of an infeasible “oracle” approximation procedure that would directly cluster the unobserved unit and time fixed effects. We follow bonhomme2022discretizing and define such approximation errors as, for all $d\in\{1,\ldots,4\}$,
and
Lemma (ref) in Section (ref) below suggests that, due to the injectivity condition (Assumption (ref)), the k-means clustering algorithm used in the first step achieves an approximation error close to the infeasible oracle k-means algorithm (that is $ B_\alpha^d(G_d)$, $B_\gamma^d(C_d)$). Next, we require this approximation error of the clustering algorithm to be small enough for the estimator to be asymptotically normal. This is subsumed in the next assumption below.
Assumption (ref) requires the oracle approximation error resulting from discretizing the unobserved heterogeneity to decrease sufficiently fast as the sample size increases. Intuitively, this condition requires the number of clusters to increase at a rate governed by the difficulty of the approximation problem, which itself depends on the dimensions of the fixed effects $K_\alpha$ and $K_\gamma$. As discussed in freeman2023linear and bonhomme2022discretizing, a precise dependence of the approximation error on $K_\alpha$ and $K_\gamma$ can be obtained under further regularity conditions on the distribution of $\alpha_i$ and $\gamma_t$.
Lemma (ref) shows that the approximation error decreases at a rate inversely proportional to the dimension of the underlying fixed effects. The assumption that $\alpha_i$ and $\gamma_t$ are i.i.d with compact support is only a sufficient condition that may not be necessary. While it may be restrictive for some applications and the result might hold under departures from this assumption, proving the validity of such an extension is beyond the scope of this paper. In the Monte Carlo study, the estimator continues to perform well when the time-specific fixed effects exhibit autocorrelation and have an unbounded support. We note that the assumption of i.i.d. fixed effects with compact support is invoked in Assumption S2(i) in bonhomme2022discretizing. Using Lemma (ref), we obtain the following corollary, which gives sufficient conditions for Assumption (ref).
Note that, when $N$ and $T$ grow at the same rate, the rate conditions of Corollary (ref) and Assumption (ref)(ref) can only hold together if $K_\alpha\le 3$ and $K_\gamma\le 3,$ so that we are imposing a restriction on the dimensions of the fixed effect spaces.
Our first asymptotic result is Lemma (ref) below. It states that the clustering algorithm groups together units (resp. time periods) with similar unit (resp. time) fixed effects, up to the oracle approximation error. A similar type of result is Lemma 1 in bonhomme2022discretizing.
Lemma (ref) suggests that injectivity ensures that if the approximation errors resulting from discretizing the unobserved heterogeneity based on the unobserved heterogeneity itself, $B_\alpha(G_d)$ and $B_\gamma(C_d)$, are small, then the approximation errors resulting from discretizing the unobserved heterogeneity based on discretizing time-series or cross-sectional averages of the data are small as well, as $N,T$ tend to infinity.
Next, we state the main result of the paper, that is, the asymptotic normality of the cross-fitted version of the two-step estimator.
Theorem (ref) justifies inference on $\beta$ based on Gaussian approximations of the asymptotic distribution. This contrasts with the properties of grouped fixed effects estimators in nonlinear likelihood models bonhomme2022discretizing. Indeed, classification noise affects the properties of second-step estimators in general through an incidental parameter bias. Theorem (ref) shows that under a linear structure and using a Neyman-orthogonal moment, one can construct an estimator that is free of such bias and thus allows the researcher to avoid using potentially computationally difficult and not proven valid bias reduction or bootstrap techniques for inference.
We now turn to discussing some theory for the data-driven selection rules for the number of clusters. Similarly to Section (ref), define $Q_g^d(G):=\frac{1}{N}\sum_{i\in\mathcal{N}_d} \left\| a_i^d-\widehat{a}^d(g_i^d)\right\|^2$ and $Q_c^d(C):=\frac{1}{T_d}\sum_{t\in\mathcal{T}_d} \left\| b_t^d-b^d( c_t^d)\right\|^2$, where the cluster centers $a^d(\cdot)$ and $b^d(\cdot)$ are formally defined in Appendix (ref). Let $\widehat{V}_g^d:= \frac{1}{N_dT_d^2}\sum_{i\in\mathcal{N}_d} \sum_{t\in\mathcal{T}_d}\left\|z_{it}- a_i^d\right\|^2$ and $\widehat{V}_c^d:=\frac{1}{N_d^2T_d}\sum_{t\in\mathcal{T }_d}\sum_{i\in\mathcal{N}_d} \left\|z_{it}- b_t^d\right\|^2$ denote estimators of the variance of $a_i^d$ and $b_t^d$, respectively. The data-driven selection rules are $\widehat{G}_d:=\min_{G\ge 1}\{G:\ Q_g^d(G)\le \widehat{V}_g^d\}$ and $\widehat{C}_d:=\min_{C\ge 1}\{C:\ Q_c^d(C)\le \widehat{V}_c^d\}$. The following lemma gives conditions under which $\widehat{G}_d$ and $\widehat{C}_d$ yield an approximation error decaying at a rate satisfying Assumption (ref).
The condition $\widehat{V}_g^d=O_P(1/T)$ and $\widehat{V}_c^d=O_P(1/N)$ is natural since $\widehat{V}_g^d$ and $\widehat{V}_c^d$ are variance estimators for time series and cross-section averages, respectively.
We consider Monte Carlo simulations to evaluate the finite sample performance of the estimator $\widehat{\beta}$ and its cross-fitted version $\widehat{\beta}^{\rm CF}$. All results in this section are averages over 10,000 replications. In all simulations, we use 30 random starting values and the Hartigan-Wong algorithm to optimize the k-means objective functions.\footnote{The results are not sensitive to the implementation of k-means.}
\paragraph{DGP.} First, we describe the data-generating processes (DGPs). We consider the sample sizes $N=50$ and $T\in\{10,20,30,40,50\}$. There is a single regressor, that is, $K=1$, and we set $\beta=1$. The unit fixed effects $\alpha_i$ are i.i.d. $\text{Gamma}(1,1)$ random variables (so that $K_\alpha=1$). The time fixed effects $\gamma_t$ are one-dimensional, that is $K_\gamma=1$, and follow an AR$(1)$ process with parameter $\rho\in \{0,0.7\}$ and disturbances drawn from a Gamma distribution with shape parameter $(1-\rho)^2/(1-\rho^2)$ and scale parameter $(1-\rho)/(1-\rho^2)$.\footnote{The process is initialized with a $\text{Gamma}(1,1)$ distribution, and we discard the first $10,000$ observations as a burn-in period.} Here, $\rho$ controls the degree of serial correlation in $\gamma_t$. When $\rho=0$, $\gamma_t$ simply follows an i.i.d. $\text{Gamma}(1,1)$ distribution.
The error terms $u_{it1}$ and $v_{it}$ also follow AR$(1)$ processes. Specifically, we set $u_{i11}\sim\mathcal{N}(0,1)$ and $v_{i1}\sim \mathcal{N}(0,1)$, and, for all $i\in\{1,\dots,N\}$ and $t\in\{2,\dots,T\}$, $$ u_{it1} =\kappa u_{i(t-1)1} + \mathcal{N}(0,(1-\kappa^2))\ \text{ and }\ v_{it} = \kappa v_{i(t-1)} + \mathcal{N}(0,(1-\kappa^2)), $$ where $\kappa$ is set to either $0$ or $0.7$ and controls the level of time-series dependence in the error terms.
For the functions $f$ and $h_1$, we consider two DGPs:
DGP 1 is inspired by the constant elasticity of substitution (CES) specification for time-varying unobserved heterogeneity proposed in bonhomme2022discretizing.
\paragraph{Estimators.} We start by evaluating the baseline estimator $\widehat{\beta}$, where the number of clusters $G$ and $C$ are chosen according to the rule outlined in Section (ref).\footnote{When $T\in\{10,20\},$ for a small fraction of replications, this rule yields values of $\widehat{G}$ and $\widehat{C}$ such that the number of degrees of freedom of the estimator is 0. To circumvent this problem, when the data-driven rule implies a number of unit clusters (resp. time clusters) larger than $4N/5$ (resp. $4T/5)$, we replace it by $4N/5$ (resp. $4T/5$). This ensures that the number of degrees of freedom remains strictly positive across all replications.} For inference, we use heteroskedasticity autocorrelation consistent standard errors clustered at the level of each unit à la Arellano1987, that is, the standard error for $\widehat{\beta}$ is $$\text{se}(\widehat{\beta}):=\sqrt{\frac{NT}{NT-NC-TG}}\left(\frac{1}{NT}\sum_{i=1}^N\sum_{t=1}^T\widehat{u}_{it1}^2 \right)^{-1}\left(\frac{1}{NT}\sum_{i=1}^N\left(\sum_{t=1}^T(\widehat{e}_{it}-\widehat{\beta}\widehat{u}_{it1})\widehat{u}_{it1}\right)^2\right)^{1/2},$$ where the factor $ \sqrt{\frac{NT}{NT-NC-TG}}$ is a degrees-of-freedom correction.
Then, in the same designs, we study the cross-fitted estimator $\widehat{\beta}^{\rm CF}$. For all $d\in\{1,\dots,4\}$, we set $G_d$ and $C_d$ in each fold according to the data-driven rule described in Section (ref). The standard errors are heteroskedasticity and autocorrelation consistent standard errors and computed as $$\text{se}(\widehat{\beta}^{\rm CF}):=\sqrt{\frac{NT}{df}}\left(\frac{1}{NT}\sum_{i=1}^N\sum_{t=1}^T\left(\widehat{u}_{it1}^{d^*_{it}}\right)^2 \right)^{-1}\left(\frac{1}{NT}\sum_{i=1}^N\left(\sum_{t=1}^T(\widehat{e}_{it}^{d^*_{it}}-\widehat{\beta}^{\rm CF}\widehat{u}^{d^*_{it}}_{it1})\widehat{u}^{d^*_{it}}_{it1}\right)^2\right)^{1/2},$$ where $d^{*}_{it}= \sum_{d=1}^4 d\mathbf{1}\{(i,t)\in\mathcal{O}_d\}$ is the fold corresponding to observation $(i,t)$ and $df:=\sum_{d=1}^4 \left(N_dT_d- N_dC_d-T_dG_d\right)$ is the number of degrees of freedom.
We compare our estimators against seven alternative approaches. The first benchmark, denoted $\widehat{\beta}^{\text{Bai}}$, corresponds to the estimator proposed by bai2009panel, using $\lfloor T^{1/2} \rfloor$ factors. As shown by freeman2023linear, this estimator is consistent under our model.
We next assess the two-step grouped fixed effects estimator introduced by freeman2023linear, denoted $\widehat{\beta}^{\text{GFE}}$. Our implementation follows their methodology, clustering only the first five loadings and factors using a hierarchical clustering procedure with a minimum single linkage algorithm.
We also include the classical two-way fixed effects estimator, denoted $\widehat{\beta}^{\rm TWFE}$, as well as the factor-augmented regression estimator, $\widehat{\beta}^{\rm FA}$, proposed by greenaway2012asymptotic and westerlund2015cross, where the number of factors is selected using the eigenvalue ratio estimator of ahn2013eigenvalue. Additionally, we consider the pooled CCE estimator, $\widehat{\beta}^{\rm CCE}$, introduced by pesaran2006estimation.
Finally, we evaluate the two estimators proposed by bonhomme2022discretizing in general parametric likelihood models with nonseparable two-way fixed effects, for which no inference results are available.\footnote{Specifically, we consider here the estimators of bonhomme2022discretizing corresponding to the special case where the likelihood is that of the linear model with Gaussian homoscedastic errors.} Let the number of clusters and the actual clusters be computed as in Section (ref). The first estimator, denoted $\widehat{\beta}^{\rm 1}$, is defined as
where $ \widetilde{e}_{it}^1 := y_{it} - \bar{y}_{g_it}$ and $\widetilde{u}_{it}^1 := x_{it} - \bar{x}_{g_it}$. This estimator has only unit cluster fixed effects, i.e., it applies a within-transformation with respect to unit clusters only. The second estimator, denoted $\widehat{\beta}^{\rm 2}$, is given by
where $ \widetilde{e}_{it}^2 := y_{it} - \bar{y}_{g_ic_t}, $ $ \widetilde{u}_{it}^2 := x_{it} - \bar{x}_{g_ic_t}, $ and, for any variable $w_{it}$, $$ \bar{w}_{g_ic_t} := \frac{1}{N_{g_i}T_{c_t}}\sum_{j=1}^N\sum_{s=1}^T \mathbf{1}\{g_j=g_i\} \mathbf{1}\{c_s=c_t\}w_{js}.$$ This estimator has interacted unit and time cluster fixed effects, that is it performs the within-transformation with respect to the interaction of unit and time clusters.
For all alternative estimators, we use unit-clustered heteroskedasticity-robust standard errors with a degrees-of-freedom correction.
\paragraph{Results.} The results for all estimators are reported in Tables (ref) and (ref). The columns “Bias" and “Var" report the estimators' bias and variance. The columns “Cov" and “Wid" present the coverage and width of the 95% confidence intervals based on a Gaussian approximation based on the aforementioned standard errors. For the estimator $\widehat{\beta}$, we also report the average values of $\widehat{G}$ and $\widehat{C}$ in the columns with the same names.
We find that the estimators, $\widehat{\beta}$ and $\widehat{\beta}^{\rm CF}$, exhibit small bias across all designs. The baseline estimator, $\widehat{\beta}$, achieves coverage levels close to the nominal 95% across nearly all sample sizes. In comparison, the cross-fitted estimator has slightly lower performance, though it remains reasonably close to $\widehat{\beta}$. In contrast, all benchmark estimators have much larger bias and variance and clearly lower coverage. Concerning $\widehat{\beta}^{\rm Bai}$ and $\widehat{\beta}^{\rm GFE}$, we conjecture that this is because the latter two estimators are designed to approximate the term $f(\alpha_i,\gamma_t)$ but not $h(\alpha_i,\gamma_t)$, thereby losing the robustness property discussed in Section (ref).\footnote{It can also be noted that the performance of $\widehat{\beta}^{\rm Bai}$ and $\widehat{\beta}^{\rm GFE}$ deteriorates when $\kappa=0.7$ instead of $0$. We conjecture that this is due to the fact that $\widehat{\beta}^{\rm Bai}$ is asymptotically unbiased and principal components analysis is fixed-T consistent under i.i.d. errors but not otherwise bai2003inferential,bai2009panel.} The estimators $\widehat{\beta}^{\rm TWFE}$, $ \widehat{\beta}^{\rm FA}$ and $ \widehat{\beta}^{\rm CCE}$ are biased because they estimate either an additive two-way fixed effects or a linear factor structure, which are not enough to capture the nonlinearities of our DGPs. Interestingly, inference based on $\widehat{\beta}^{\rm 1}$ and $\widehat{\beta}^{\rm 2}$ does not seem to be correct, demonstrating the advantage of our proposal using additive two-way group fixed effects in the second step. Based on these results, we recommend that practitioners primarily use the baseline estimator, $\widehat{\beta}$.
These findings affirm that cross-fitting primarily serves as a theoretical construct to facilitate asymptotic proofs, offering limited practical benefits in finite samples. The slightly weaker performance of $\widehat{\beta}^{\rm CF}$ can be intuitively attributed to its use of only half the observations for clustering. This robustness underscores the practical value of our approach in such contexts. Moreover, our results reveal that the estimators maintain strong performance under time-series dependence, suggesting that the i.i.d. assumption on the errors is also primarily a theoretical convenience.
Appendix (ref) contains additional simulation results. First, in Appendix (ref), we study, in the same simulation designs, the finite-sample performance of the estimators when $N=500$ and $T\in\{10,20,30,40,50\}$. We also find that $\widehat{\beta}$ and $\widehat{\beta}^{\rm CF}$ have very good performance in such a large $N$ small $T$ setting common in real-world datasets and outperform the alternatives. Second, in Appendix (ref), we study the sensitivity of the baseline estimator to the number of clusters and find that it is remarkably robust.
We revisit james2015us, focusing on the impact of increases in resource-based government revenues on various fiscal outcomes across U.S. states: non-resource tax revenues, income tax revenues, total expenditures, education expenditures, and public savings. The data comprises annual government revenues and expenditures, as well as private income for all U.S. states over the period from 1958 to 2008, so that $N=50$ and $T=51$.\footnote{The full dataset is available at \url{https://www.openicpsr.org/openicpsr/project/114577/version/V1/view}.}
As argued by james2015us, following economic theory, resource-based tax revenue should have a negative effect on nonresource revenue and income tax revenue, but a positive impact on total expenditure, education expenditure, and savings. A potential confounding factor is the business cycle $\gamma_t$, which can influence both non-resource and resource-based revenues. During periods of high macroeconomic output, energy consumption and private income tend to rise, leading to higher revenues. This relationship can introduce omitted variable bias. Our estimation approach addresses this issue by allowing the effect of business cycles to vary nonlinearly across states and revenue types through $\alpha_i$, reflecting differences in tax schemes and economic structures. Arguably, the effect of most unobserved state-specific characteristics such as average population density, political preferences, wealth, unemployment, culture, and institutional quality deemed time-invariant in james2015us might actually vary over 51 years.
\paragraph{Regression results.} james2015us employs the within-estimator for a two-way fixed effects model, denoted $\widehat{\beta}^{\text{TWFE}}$, regressing the ratio of the various outcomes to private income in the state-year on the ratio of resource-based government revenues to private income. We consider seven estimators: $\widehat{\beta}$, $\widehat{\beta}^{\rm CF}$, $\widehat{\beta}^{\text{TWFE}}$, $\widehat{\beta}^{\rm Bai}$, $\widehat{\beta}^{\rm GFE}$, $\widehat{\beta}^{\text{FA}}$ and $\widehat{\beta}^{\text{CCE}}$. Note that for all estimators but $\widehat{\beta}^{\text{TWFE}}$, we first standardize the outcome variable and regressor before applying the methods. The estimated coefficients and standard errors are then rescaled to correspond to the original model. All estimators and their standard errors are computed as in the simulations of Section (ref), except that we use $10,000$ initializations for the kmeans algorithms of $\widehat{\beta}$ and $\widehat{\beta}^{\rm CF}$. Table (ref) reports the results for all outcomes and estimators. Table (ref) presents the values of $\widehat{G}$ and $\widehat{C}$ for the different outcomes.
The proposed estimator $\widehat{\beta}$ always has the sign predicted by economic theory. The results for $\widehat{\beta}$ are also significant for 4 of the 5 outcomes. In contrast, each of the alternative estimators has a sign in disagreement with the theory for at least one outcome. For the outcome “Savings", $\widehat{\beta}^{\rm GFE}$ and $\widehat{\beta}^{\rm CCE}$ yield estimates larger than $1$, which are difficult to justify from an economic perspective. Overall, the estimator's conclusions often differ from the alternatives, demonstrating its ability to provide unique insights.
\paragraph{Clusters for the outcome “Savings."} We now present the clusters obtained by our method for the outcome “Savings."\footnote{We selected this outcome because it has the fewest clusters, making it easier to represent on a map.} Figure (ref) displays the 5 unit clusters on a map, with their centers listed in Table (ref). These clusters represent states with similar average savings and resource-based government revenues over the period. While there is no reason to expect that they should correspond to geographically close states (this assumption is not imposed in the data-driven estimation), it turns out that some geographical dependence is effectively captured as geographically close states often end up in the same estimated cluster. Interestingly, cluster 5 corresponds to Alaska, and cluster 4 consists of New Mexico and Wyoming. These three states are known to be particularly rich in natural resources. The information for the time clusters is given in Figure (ref) and Table (ref). We find similar patterns, with clusters seemingly capturing the business cycle.
Overall, our results indicate that controlling for flexible patterns of time-varying unobserved heterogeneity does not refute the predictions made by economic theory. Estimated clusters confirm that unobserved heterogeneity is both spatially and temporally correlated.
This paper shows how to use Neyman-orthogonal moments to build inference tools after discretizing time-varying unobserved heterogeneity in linear panel data models. The proposed procedure is intuitive and simple, but nevertheless exhibits excellent asymptotic properties and finite-sample performance. A natural extension is to consider heterogeneous slope parameters. While adapting the proposed estimation procedure to accommodate either unit- or time-specific slope coefficients is relatively straightforward, we leave the study of more flexible unit- and time-varying structures for further research.
\setcounter{section}{0}
As for $\widehat{\beta}$, we estimate the group memberships in the first step via a clustering method and compute an OLS estimator in the second step. The main difference is that the data used in each of these two steps do not intersect but the final estimator still uses variation across the full dataset. The estimation procedure to obtain the resulting cross-fitted two-way grouped fixed effect estimator is as follows. For each fold, $d\in\{1,\dots,4\}$:
The final estimator is the linear regression of the $\widehat{e}_{it}^d$ on the $\widehat{u}_{it}^d$, $$\widehat{\beta}^{\rm CF}:=\left(\sum_{d=1}^4\sum_{(i,t)\in \mathcal{O}_d} \widehat{u}_{it}^d (\widehat{u}_{it}^d)^\top \right)^{-1}\sum_{d=1}^4\sum_{(i,t)\in \mathcal{O}_d}\widehat{u}_{it}^d \widehat{e}_{it}^d, $$ which is numerically equivalent to
In summary, to obtain the unit cluster indicators $g_i^d\in\{1,\dots,G_d\}$ (resp. the time cluster indicators $c_t^d\in\{1,\dots,C_d\}$), we use the fold that contains the same units as $\mathcal{O}_d$ but different dates (resp. the same dates as $\mathcal{O}_d$ but different units). A similar trick is used for clustering time periods. We then use these clusters to estimate $e_{it}$ and $u_{it}$, before running a linear regression on such estimates.\footnote{As for our baseline estimator, the second step (ref) of our cross-fitted estimator corresponds to the second step of the cross-fitted estimator in freeman2023linear. The clustering steps differ between the two papers.}
The clustering steps are carried out using straightforward adapted versions of the algorithm introduced in Section (ref) (indexing all relevant sample $d$-dependent variable by $d$) , which we display below for completeness.
Clustering algorithm for units. Let the empirical averages $ a_i^d:=\frac{1}{T_{\tilde{d}}} \sum_{t\in \mathcal T_{\tilde{d}}} z_{it}, \ i\in\mathcal{N}_d$ be computed on fold $\tilde{d}$. We use the algorithm $$\left(\widehat{a}^d(1),\dots,\widehat{a}^d(G_d),\{g_i^d,\ i\in\mathcal{N}_{d}\}\right)\in\operatorname*{arg\,min}\limits_{
}\sum_{i\in\mathcal{N}_{d}} \left\| a_i^d-a(g_i)\right\|^2. $$
Clustering algorithm for dates. .Let the empirical averages $ b_t^d:=\frac{1}{N_{\tilde{d}}} \sum_{i\in \mathcal{N}_{\tilde{d}}} z_{it},\ t\in\mathcal{T}_d$ be computed on fold $\tilde{d}$. We use the algorithm $$\left(\widehat{b}^d(1),\dots,\widehat{b}^d(C_d),\{c_t^d,\ t\in\mathcal{T}_{d}\}\right)\in\operatorname*{arg\,min}\limits_{
}\sum_{t\in\mathcal{T}_{d}} \left\| b_t^d-b(c_t)\right\|^2. $$ In practice, we use the data-driven rule outlined in Section~\ref{sec:est} to select the number of clusters $G_d$ and $C_d$ in the different folds $d\in\{1,\dots,4\}$.
In this section, we describe an alternative clustering algorithm for the first step based on the pseudo-distance of zhang2017estimating and hierarchical clustering as in mugnier2024simple. Let us explain how units are clustered with this approach. For $i,j\in\{1,\dots,N\}$, we define the pseudo-distance $$\widehat{d}_{\infty,1}(i,j)=\frac1T\max_{\ell \in\{1,\dots,N\}\backslash\{i,j\}}\left(\left|\sum_{t=1}^T(y_{it}-y_{jt})y_{\ell t}\right|+\sum_{k=1}^K\left|\sum_{t=1}^T(x_{itk}-x_{jtk})x_{\ell tk}\right|\right).$$ Then, to obtain the unit clusters, we apply a hierarchical clustering algorithm to the $N\times N$ matrix $\widehat{D}$ such that $\widehat{D}_{ij}=\widehat{d}_{\infty,1}(i,j).$ See mugnier2024simple for a formal presentation. Similarly to mugnier2024simple, we choose the threshold $c_{NT}$ for the maximum intragroup distance equal to $$1.35\frac{\log(T)}{K\sqrt{\min(N,T)}}\check{\sigma},$$ where
To avoid having 0 degrees of freedom, if this value of $c_{NT}$ gives more than $\lfloor 2N/5\rfloor $ clusters, we set the number of unit clusters to $\lfloor 2N/5\rfloor $. The time clusters are obtained symmetrically. The proposed approach circumvents averaging the data before clustering, which should allow avoiding the injectivity assumption; see also the discussion in athey2025identification.
We consider baseline $\widetilde{\beta}$ and cross-fitted $\widetilde{\beta}^{\rm CF}$ adaptations of the estimators that use hierarchical clustering with an average linkage function on the pseudo-distance instead of k-means to obtain the clusters. In Table (ref), we report the results of simulations with these alternative estimators, where the DGPs and standard errors are as in Section (ref). Results suggest that $\widetilde{\beta}$ has poor performance, while $\widetilde{\beta}^{\rm CF}$ has almost nominal coverage but large confidence intervals compared to the version using k-means.
We present simulation results for the estimators when $N=500$ and $T\in\{10,20,30,40,50\}$ under the data-generating processes outlined in Section (ref). The results are reported in Tables (ref) and (ref).
We also study the sensitivity of the estimator to the number of groups. To do so, we implement simulations under $N=T=50$ in DGP 1 and 2 (with $\kappa=\rho=0$) of Section (ref). We simulate 10,000 datasets and compute the value of $\widehat{\beta}$ with a number of unit and time clusters $G=C$ varying between 1 and 24. The estimator and its standard error are computed as in Section (ref). Figures (ref) and (ref) present the average bias, variance, coverage and width of 95% confidence intervals, in DGP 1 and 2, respectively. As long as the number of groups is larger than 5, the bias and coverage are insensitive to $G=C$. However, as $G=C$ increases, the variance and the width also become larger.
We prove only the first statement, as the argument for the second is analogous. The proof proceeds in two steps. \paragraph{Step 1.} In this step, we establish that
By the triangle inequality and the classical inequality $ab\le (a^2+b^2)/2$ for all $a,b\in{\mathbb{R}}$,
Under Assumption (ref)(ref), arguments analogous to those in the proof of Lemma 1 in bonhomme2022discretizing yield
Next, using that $\widehat a^d(g_i^d)= \frac{1}{N^d_{g_i^d}}\sum_{j\in\mathcal N_d} \mathbf{1}\{g_j^d= g_i^d\}a_j^d $, we obtain
where the inequality follows from the triangle inequality, and the last equality is a consequence of Assumption (ref). Combining (ref)-- (ref), we obtain (ref).
\paragraph{Step 2.} In this step, we establish the result stated in the lemma. We first note that
where the first inequality follows from the triangle inequality and the second inequality is a consequence of the Cauchy--Schwarz inequality. Next, by Assumption (ref)(ref), there exists a constant $L>0$ such that
Moreover, we have
where in the last equality we used (ref). Combining the last result with (ref)--(ref), we obtain the result of the lemma.
This section concerns the proof of Theorem (ref). It is organized as follows. Section (ref) introduces the notation used in the proof. Section (ref) contains the main body of the proof of Theorem (ref), which relies on auxiliary lemmas stated and proved in Section (ref). The proofs of these auxiliary lemmas, in turn, depend on technical lemmas stated and proved in Section (ref).
For all $(i,t,k,d)\in\{1,\dots,N\}\times\{1,\dots,T\}\times\{1,\dots,K\}\times\{1,\ldots,4\}$, we let $h_{itk}:=h_k(\alpha_i,\gamma_t)$, $f_{it}:=f(\alpha_i,\gamma_t)$, and we use the notation
We have
Since $y_{it}=x_{it}^\top\beta+f_{it}+v_{it}$, this yields
We obtain
By Lemmas (ref)-- (ref) and the continuous mapping theorem, we obtain
Under Assumption (ref), by combining H\"older's and Markov's inequalities, it is not difficult to show that the following conditional Lindeberg condition holds: for all $\varepsilon>0$, as $N,T$ tend to infinity, \[ \frac{1}{NT}\sum_{i=1}^N\sum_{t=1}^T{\mathbb{E}}\left[\Vert v_{it}u_{it}\Vert^2\mathbf 1\left\{\Vert v_{it}u_{it}\Vert\geq \varepsilon\sqrt{NT}\right\}|\mathcal F_{NT}\right]\to 0. \] An application of the multivariate Lindeberg--Feller central limit theorem yields, for all $c\in{\mathbb{R}}^K$ and $z\in{\mathbb{R}}$, almost-surely, as $N$ and $T$ tend to infinity,
By the dominated convergence theorem, the sequence of unconditional cumulative distribution functions of $c^\top\Omega^{-1/2}\left(\frac{1}{\sqrt{NT}} \sum_{i=1}^N \sum_{t=1}^Tu_{it}v_{it}\right)$ evaluated at $z$ converges to the same limit. By the Cramer--Wold device, this yields, as $N$ and $T$ tend to infinity, \[ \frac{1}{\sqrt{NT}} \sum_{i=1}^N \sum_{t=1}^Tu_{it}v_{it}\overset{d}{\to}\mathcal N\left(0,\Omega\right). \] In particular, $\frac{1}{\sqrt{NT}} \sum_{i=1}^N \sum_{t=1}^Tu_{it}v_{it}=O_P(1)$ so that
and the result follows from Slutsky's lemma.
Fix $d\in\{1,\dots,4\}$. We only show the result for $\widehat{G}_d$; the proof for $\widehat{C}_d$ is similar and, therefore, omitted. We have
Following the arguments of Step 2 of the proof of Lemma (ref), we obtain that there exists $L>0$ such that
Moreover, for all $\tilde g_i\in\{1,\dots,G_d\},\ i\in\mathcal{N}_{d} $, we have
where we used the triangle inequality and the classical inequality $ab\le (a^2+b^2)/2$. By Assumption (ref)(ref), this yields $ B_\alpha^d(G_d) \le 6LQ_g^d(G_d)+ O_P\left(\frac{r_\alpha}{T}\right). $ Since $ Q_g^d(\widehat{G}_d)\le \widehat{V}_g^d=O_P(1/T)$, we obtain, by Assumption (ref)(ref), $$B_\alpha(\widehat{G}_d)= O_P\left(\frac{r_\alpha}{T}\right)=o_P\left(\frac{1}{(NT)^{1/4}}\right).$$