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.
107,390 characters · 19 sections · 94 citation commands
Inference on Linear Regressions with Two-Way Unobserved Heterogeneity
We are interested in inference on $\beta $ in the model,
where $X_{it}\in \mathbb{R}^{K}$ are a set of observed covariates, $\alpha _{i}$ and $\gamma _{t}$ are unobserved finite--dimensional fixed effects that potentially correlate with $X_{it}$, and $\varepsilon _{it}$ is an idiosyncratic shock that is mean--independent of $X_{i,t}$, $\alpha _{i}$ and $\gamma _{t}$. The unknown function $g(\cdot ,\cdot )$ adheres to certain smoothness conditions but is otherwise left unrestricted and treated as a nonparametric object.
Our framework includes as special cases well-known parametric models such as the two-way additive fixed-effects model, $g(\alpha _{i},\gamma _{t})=\alpha _{i}+\gamma _{t}$ and the interactive fixed-effects model, $g(\alpha _{i},\gamma _{t})=\sum_{r=1}^{R}u_{r}\left( \alpha _{i}\right) v_{r}\left( \gamma _{t}\right) =\sum_{r=1}^{R}\lambda _{ir}f_{tr}$, where $\lambda _{ir}=\sigma _{r}u_{r}(\alpha _{i})$ and $f_{tr}=v_{r}(\gamma _{t})$. In both cases, the (transformed) individual effects and time effects can be controlled for with, e.g., the least square estimator Bai2009,MoonWeidner2015, or correlated common effects estimator Pesaran2006 and the resulting estimator of $\beta $ is $\sqrt{NT}$ --asymptotically normally distributed.
However, the nonparametric case presents a challenging scenario that does not trivially lead to $\sqrt{NT}$-consistent estimates for $\beta $. For example, freeman2023linear develop least-squares estimator of $\beta $ where a sieve--estimator of $g$ is employed, but find that the resulting estimator of $\beta $ does not enjoy $\sqrt{NT}$-asymptotic normality; see also fernandez2021low.
The contribution of this paper is two--fold: First, we propose a general estimation procedure for $\beta $ that takes as input any first-step estimators of $g$ and $(\alpha _{i},\gamma _{t})$. Under weak conditions on the first--step estimators, we show that the resulting estimator of $\beta $ is $\sqrt{NT}$-asymptotically normally distributed, thereby allowing for standard inference tools to be employed. Second, we propose novel estimators of $g$ and $(\alpha _{i},\gamma _{t})$ and verify that these satisfy the regularity conditions for the theory of our general estimation procedure to hold.
The general estimation procedure relies on two main ingredients: First, we follow freeman2022multidimensional and develop moment conditions for the estimation of $\beta $ that are Neyman orthogonal w.r.t. to $g$; see, e.g., chernozhukov2022locally for an overview of this method in the context of semiparametric estimation. Unsurprisingly, the Neyman orthogonal moment conditions take the same form as the ones used for estimation of the partially linear model, where $\left( \alpha _{i},\gamma _{t}\right) $ are treated as observed co--variates; see, e.g., robinson1988root. In order to operationalise these moment conditions, the researcher will have to plug in first-step estimators of the nonparametric functions $g_{X}\left( \alpha _{i},\gamma _{t}\right) :=\mathbb{E}\left[ X_{it}|\alpha _{i},\gamma _{t}\right] $ and $g_{Y}\left( \alpha _{i},\gamma _{t}\right) :=\mathbb{E} \left[ Y_{it}|\alpha _{i},\gamma _{t}\right] $ together with estimators of the unknown fixed effects $\left( \alpha _{i},\gamma _{t}\right) $. In the ideal scenario where $\left( \alpha _{i},\gamma _{t}\right) $ are observed, the Neyman orthogonality of the moment conditions w.r.t. $g_{Z}:=\left( g_{X},g_{Y}\right) $ ensure that the nonparametric estimation of $g_{Z}$ has no first-order effect on the estimator of $\beta $. In particular, it ensures that the estimated fixed effects used in the estimation of $g_{Z}$ do not contribute to the asymptotic variance the estimator of $\beta $.
However, the resulting estimator of $\beta $ will generally suffer from well-known incidental parameter biases due to estimated fixed effects. The leading term of these biases could be adjusted for using standard techniques, e.g., analytical bias adjustment or the Jackknife. The second ingredient of our procedure avoids any such post--estimation bias adjustment by combining sample splitting and a generalised version of the two--way nonparametric regression estimator of freeman2023linear that removes the incidental parameter biases. Under weak regularity conditions on the first-step nonparametric estimator, we show that the resulting estimator of $ \beta $ will be $\sqrt{NT}$-asymptotically normally distributed and with the same distribution as the oracle estimator where $g$ is known and $\left( \alpha _{i},\gamma _{t}\right) $ are observed.
Next, we develop a specific estimator of $\left( \alpha _{i},\gamma _{t}\right) $ and $g_{Z}$ that satisfy the regularity conditions for our general estimation theory to hold. These novel estimators impose weak conditions on $g_{Z}$ and the fixed effects. Under weak regularity conditions, FreemanKristensen2026 show that a finite number of leading eigenfunctions of $g$ evaluated at $\left( \alpha _{i},\gamma _{t}\right) $ can be used as proxies for $\left( \alpha _{i},\gamma _{t}\right) $. Our proposed estimator is then obtained in two steps: In the first step, we employ the approximate factor model estimator of freeman2023linear to obtain estimates of the leading eigenfunctions of $g$ . In the second step, we run a nonparametric regression of $Y_{it}$ and $ X_{it}$, respectively, onto the estimated leading eigenfunctions from the first step to obtain our final estimator of $g_{Z}$. The nonparametric regression approach is related to the nonparametric smoothing technique employed in freeman2024multidimensional.
The nonparametric regression in the second step comes in two versions: The first version constructs a multivariate index of the eigenfunctions that are then used as covariates in a nonparametric regression. The second one employs an nonparametric additive regression procedure with the attractive feature of not suffering from any curse--of--dimensionality. The two estimators rely on different assumptions on $g_{Z}$ and come with different convergence rates: The first one imposes weaker conditions and comes with smaller biases but bigger variances compared to the second one. We demonstrate that both of the two estimators exhibit sufficiently fast rate of convergence so that they can be combined with the Neyman orthogonal moment conditions to obtain the desired result.
Alternative estimators of $\left( \alpha _{i},\gamma _{t}\right) $ and $ g_{Z} $ are proposed in deaner2025inferring and beyhum2025inference with the latter using their proposal to estimate $ \beta $ based on the Neyman orthogonal moments described above. The regression estimator of deaner2025inferring is closest in spirit to ours but employ a different proxy for the fixed effects, namely an estimated pseudo--distance, while beyhum2025inference use sample moments of $ Z_{it}$ to estimate the fixed effects and then $k$-means methods to estimate $g_{Z}$. We expect our estimators of $g_{Z}$ to come with smaller biases compared to the one of beyhum2025inference if $g_{Z}$ is a smooth function. In addition, the formal assumptions under which the fixed effects estimators of deaner2025inferring and beyhum2025inference are valid are different from ours. As such, the three papers complement each other.
One particularly attractive feature of our proposed estimator of $\beta $ over the one of beyhum2025inference is that it will in general be more efficient:\ Suppose that
where $(\tilde{\alpha}_{i},\tilde{\gamma}_{t})$ are additional latent individual and time effects that are not present in ((ref)). The proposed algorithm of beyhum2025inference will control for both $ g_{X}(\alpha _{i},\gamma _{t})$ and $\tilde{g}_{X}(\tilde{\alpha}_{i},\tilde{ \gamma}_{t})$ in its first step and so in effect use an estimator of $e_{it}$ as regressors in their proposed estimator of $\beta $. In contast, our procedure uses an estimator of $X_{it}-g_{X,1}(\alpha _{i},\gamma _{t})- \mathbb{E}\left[ g_{X,2}(\tilde{\alpha}_{i},\tilde{\gamma}_{t})|\alpha _{i},\gamma _{t}\right] $ as regressors in estimation of $\beta $. Importantly, the variance of the latter will be larger than the one of $ e_{it}$ and so leads to lower asymptotic variance of our estimator of $\beta $ compared to the one of beyhum2025inference. Moreover, our procedure suffers from a lower curse--of--dimensionality since it only involves learning about/estimating $(\alpha _{i},\gamma _{t})$; in contrast, the one of beyhum2025inference requires learning about $(\tilde{\alpha}_{i}, \tilde{\gamma}_{t})$ in addition. Thus, unless $\tilde{g}_{X}(\tilde{\alpha} _{i},\tilde{\gamma}_{t})=0$, our procedure should dominate theirs both in small and large samples.
We carry out an extensive simulation study that support our theoretical results: The estimator of $\beta $ that relies on our novel eigenfunction--based estimators of the fixed effects suffers from only small finite-sample biases and dominates the estimators of $\beta $ that takes as input fixed effects estimators based on aforementioned pseudo-distance and sample moments across a range of different DGP's. In particular, in scenarios where $(\alpha _{i},\gamma _{t})$ are multivariate, our proposal significantly outperforms these alternative fixed effects estimators.
The remainder of the paper is organised as follows: In Section (ref) , we present the general theory for two--step estimation of $\beta $ that allows for a broad range of first--step estimators of $g_{Z}$ and fixed effects. We present our novel estimators of $g_{Z}$ and fixed effects in Section (ref). The asymptotic properties of this estimator are analysed in Section (ref). The results of our simulation study are presented in Section (ref) and we conclude in Section (ref). All proofs have been relegated to the Appendix.
We here first present our general estimation approach, the Neyman Orthogonalised estimator of $\beta $, and show that it can achieve $\sqrt{NT} $-asymptotic normality under weak conditions on the first-step estimation of the nonparametric component $g$. However, the estimator will generally suffer from incidental parameter biases. We show how a Jackknife-type procedure can be used to remove such biases when combined with sample splitting.
Recall the following definitions,
and write $\Gamma _{X,it}=g_{X}\left( \alpha _{i},\gamma _{t}\right) $ and $ \Gamma _{Y,it}=g_{Y}\left( \alpha _{i},\gamma _{t}\right) $ for brevity. Note that $\Gamma _{Y,it}=\Gamma _{X,it}^{\prime }\beta +\Gamma _{it}$, where $\Gamma _{it}:=g\left( \alpha _{i},\gamma _{t}\right) $. Moreover,
holds by definition of $g_{X}\left( \alpha _{i},\gamma _{t}\right) $. For ease of notation, we let $Z_{it}=\left( Y_{it},X_{it}\right) $, $ g_{Z}=\left( g_{Y},g_{X}\right) $, $\Gamma _{Z,it}=\left( \Gamma _{Y,it},\Gamma _{X,it}^{\prime }\right) ^{\prime }$ and $\varepsilon _{Z,it}=\left( \varepsilon _{it},\eta _{it}\right) $ so that $Z_{it}=\Gamma _{Z,it}+\varepsilon _{Z,it}$, where $\mathbb{E}[\varepsilon _{Z,it}|\alpha _{i},\gamma _{t}]=0$.
Our estimator of $\beta $ will be based on the following moment function,
where we note that
Note that this is the same moment function used in the estimation of the partially lineear model with cross-sectional data, c.f. robinson1988root. It is easily checked that the data--generating parameter, denoted by $\beta _{0}$, is identified as the unique solution to $ \mathbb{E}\left[ m\left( Z_{it},\beta ,\Gamma _{Z,it}\right) \right] =0$ w.r.t. $\beta $,
under the following assumptions:
We here impose strict exogeneity. Identification holds under the weaker condition of contemporaneous exogeneity, $\mathbb{E}[\varepsilon _{Z,it}|X_{it},\alpha _{i},\gamma _{t}]=0$, but parts of our theoretical analysis will make use of strict exogeneity to simplify the arguments.
Importantly, $m\left( Z_{i,t},\beta ,\Gamma _{Z,it}\right) $ is Neyman-orthogonal w.r.t. $\Gamma _{Z,it}$ at $\beta =\beta _{0}$,
An important consequence of this feature is that estimation of $\beta $ based on these moment conditions is less sensitive to the first-step estimation of $\Gamma _{Z,it}$.
To develop such estimator, we here take as given any first--step estimators $ \hat{\Gamma}_{Z}=(\hat{\Gamma}_{Y},\hat{\Gamma}_{X})$ of $\Gamma _{Z}$ as chosen by the researcher. We then define our Neyman--orthogonal estimator of $\beta $ as the solution to $\sum_{it}m(Z_{it},\hat{\beta}_{NO},\hat{\Gamma} _{Z,it})=0$, where $\sum_{it}:=\sum_{i=1}^{N}\sum_{t=1}^{T}$, which takes the form
Due to the moment conditions being orthogonal w.r.t. $\hat{\Gamma}_{Z}$, this estimator of $\beta $ will be $\sqrt{NT}$-asymptotically normally distributed convergence as long as $\hat{\Gamma}_{Z}-\Gamma _{Z}=o_{P}\left( (NT)^{-1/4}\right) $. This is a well-known result from the literature on semiparametric estimators.
However, in our setting, above rate result is in general not achievable when we treat $\alpha _{i}$ and $f_{t}$ as unknown fixed effects since in this case $\hat{\Gamma}_{Z}$ will involve estimates of these. To see this, suppose that in fact $g$ in eq. ((ref)) is known to us, in which case we have fully parametric fixed effects model whose parameters could be estimated by $(\hat{\beta},\{\hat{\alpha}_{i}\}_{i=1}^{N},\{\hat{\gamma} _{t}\}_{t=1}^{T})=\arg \min_{\beta ,\{\alpha _{i}\}_{i=1}^{N},\{\gamma _{t}\}_{t=1}^{T}}\sum_{it}(Y_{it}-X_{it}^{\prime }\beta -g(\alpha _{i},\gamma _{t}))^{2}$. This in turn yields the estimator $\hat{\Gamma} _{it}=g(\hat{\alpha}_{i},\hat{\gamma}_{t})$. Except in a few special cases, such as when $g(\alpha _{i},\gamma _{t})$ is linear or multiplicative in its arguments, it is well-known that $\hat{\beta}$ will generally suffer from incidental parameter biases, $\mathbb{E[}\hat{\beta}]\simeq \beta _{0}+B_{\alpha }/T+B_{\gamma }/N$, where $B_{\alpha }$ and $B_{\gamma }$ are due to the incidental parameter biases caused by $\{\hat{\alpha} _{i}\}_{i=1}^{N}$ and $\{\hat{\gamma}_{t}\}_{t=1}^{T}$, respectively; see, e.g., Theorem 4.1 of FernandezValWeidner2016. Obviously, in our setting with $g$ unknown, we cannot hope to do better than in the parametric submodel with $g$ known, unless we are willing to impose further restrictions on the unknown fixed effects $\alpha _{i}$ and $\gamma _{t}$ or the mapping $g$.
In light of above, we develop a general asymptotic result for $\hat{\beta} _{NO}$ when $\hat{\Gamma}_{Z,it}=\hat{g}_{Z}(\hat{\alpha}_{i},\hat{\gamma} _{t})$ contains non--negiglible biases due to $\left( \hat{\alpha}_{i},\hat{ \gamma}_{t}\right) $ being used in place of $(\alpha _{i},\gamma _{t})$. We first develop an expansion of $\hat{\beta}_{NO}$ w.r.t. $\hat{\Gamma}_{Z}$: Applying the mean--value theorem to $\sum_{it}m(Z_{it},\hat{\beta}_{NO},\hat{ \Gamma}_{Z,it})=0$,
where
A second--order Taylor expansion of the right--hand side of above equation w.r.t. $\hat{\Gamma}_{X,it}$ at $\Gamma _{X,it}$ yields
where, assuming $\sum_{it}\left\Vert \eta _{it}\right\Vert ^{2}/\left( NT\right) \rightarrow ^{p}\mathbb{E}\left[ \left\Vert \eta _{it}\right\Vert ^{2}\right] $ and $\sum_{it}\left\Vert \hat{\Gamma}_{it}-\Gamma _{X,it}\right\Vert ^{2}/\left( NT\right) =o_{P}\left( 1\right) $,
and similar for the other terms. Thus, $\frac{1}{NT}\sum_{it}(X_{it}-\hat{ \Gamma}_{X,it})(X_{it}-\hat{\Gamma}_{X,it})^{\prime }=\frac{1}{NT} \sum_{it}\eta _{it}\eta _{it}^{\prime }+o_{P}\left( 1\right) $.
Next, another second--order Taylor expansion combined with $Y_{i,t}-\Gamma _{Y,it}=\beta _{0}^{\prime }\eta _{it}+\varepsilon _{i,t}$ gives us
where the last three terms contain the first and second--order effects of $ \hat{\Gamma}_{Z,it}$. The following result provides conditions under which $ \hat{\beta}_{NO}$ is asymptotically normally distributed but suffers from incidental parameter biases due to the presence of these terms:
The result is quite general and allows for a broad class of estimators ${ \hat{\Gamma}}_{Z,it}$. The main restriction is the implicit assumption that the leading bias components of this estimator are of order $O\left( 1/T\right) +O\left( 1/N\right) $. The $o_{P}(1/\sqrt{NT})$ terms in equations ((ref))--((ref)) capture both the full error from the nonparametric estimation of $g$ together with the variance components of the fixed effects estimators. As we shall see, this rate result will hold under weak restrictions on $g$, $\alpha _{i}$ and $ \gamma _{t}$ and their estimators. The asymptotic biases $B_{\alpha }/T$ and $B_{\gamma }/N$ will be present due to first--order biases of the estimators of $\alpha _{i}$ and $\gamma _{t}$ together with potentially covariation between their estimation errors and $(\eta _{it},\varepsilon _{it})$, c.f. Theorem 4.1 of FernandezValWeidner2016.
Below, we first give more primitive conditions under which ((ref))--((ref)) will hold. Next, we show how sample splitting combined with bias adjustment can be employed to remove the incidental parameter biases.
We here focus on a particular class of first-step estimators of $\Gamma _{Z}$ that combine fixed-effects estimation with nonparametric regression techniques and for these provide a partial characterisation of the incidental parameter biases introduced in ((ref))--((ref)). The class of estimators and accompanying theory includes as a special case the particular estimators that we develop in Section (ref). The estimators of beyhum2025inference also fit into our general framework. However, parts of our analysis requires the estimator of $ g_{Z}$ to be sufficiently regular (smooth) while their procedure involves discretisation and so these parts do not apply to their estimator.
We first note that, without further normalisations and restrictions on the model, we cannot separately identify $g_{Z}$, $\alpha _{i}$ and $\gamma _{t}$ . We will here assume that suitably identifying restrictions/normalisations have been developed so that the following conditions are satisfied:
First, we will require that there exists a normalised version of $g_{Z}$, denoted $g_{0,Z}$, and associated normalised version of the fixed effects $ \left( \alpha _{i},\gamma _{t}\right) $, denoted $\left( \lambda _{i},f_{t}\right) $, so that
Second, we will require that fixed-effects estimators $(\hat{\lambda}_{i}, \hat{f}_{t})$ of $\left( \lambda _{i},f_{t}\right) $ are available to us. The precise forms of the normalised versions and their estimators depend on the assumptions the researcher is willing to impose on the model. One example of such is provided in Section (ref); see FreemanKristensen2026 for further details. We will focus on the ideal setting where the estimators are regular in the sense that they satisfy
where $\mathbb{E}\left[ \psi _{\lambda }\left( Z_{it}\right) \right] = \mathbb{E}\left[ \psi _{f}\left( Z_{it}\right) \right] =0$ and $b_{\lambda ,i}$ and $b_{f,t}$ capture the leading bias terms of the two estimators; see FernandezValWeidner2016 for sufficient conditions for above to hold in a parametric setting. The fixed-effects estimators of beyhum2025inference also satisfy above with $b_{\lambda ,i}=b_{f,t}=0$.- Under above assumption, we will now provide a characterisation of the biases appearing in ((ref))--((ref)). We expect all qualitative conclusions reached in the following to also apply to irregular estimators but the precise forms of the biases contained in the first and second-order terms in ((ref))--((ref)) will in this case become more complicated, including their rates which no longer will be of the parametric kind.
Given this set-up, we will then consider estimators of $\Gamma _{Z,it}$ that runs a nonparametric regression of $Z_{it}$ onto the estimated fixed effects $(\hat{\lambda}_{i},\hat{f}_{t})$ to obtain $\hat{g}_{0,Z}$. This estimator could take many forms; one example would be kernel regression, in which case
where $K_{j}$, $j=1,2$, are kernel functions and $h_{\lambda },h_{f}>0$ are bandwidths, but the subsequent theory allows for other nonparametric regresion techniques. It is here important to note that the estimated fixed effects enter $\hat{g}_{0,Z}(\hat{\lambda}_{i},\hat{f}_{t})$ through two channels: First, as the values $(\hat{\lambda}_{i},\hat{f}_{t})$ of the argument $\left( \lambda ,f\right) $ in $\hat{g}_{0,Z}\left( \lambda ,f\right) $; second, as the regressors used to compute $\hat{g}_{0,Z}\left( \lambda ,f\right) $. Each channel could potentially generate incidental parameter biases. The first channel is well-known, c.f. FernandezValWeidner2016, while the second one is novel and has only been analysed in the context of nonparametric estimation of the density of one--way fixed effects; see, e.g., Okui2019 and Barras2021.
We first analyse the first--order terms appearing in eqs. ((ref))--((ref)). For a given function $g_{Z}$, define
Let $\hat{g}_{0,Z}$ denote the feasible nonparametric estimator of $g_{0,Z}$ that takes as input $\{\hat{\lambda}_{i}\}_{i=1}^{N}$ and $\{\hat{f} _{t}\}_{t=1}^{T}$, such as the kernel regression estimator in ((ref)), while $\hat{g}_{0,Z}^{\ast }$ is the infeasible version that takes as input $\{\lambda _{i}\}_{i=1}^{N}$ and $\{f_{t}\}_{t=1}^{T}$. We can then decompose the vector of first-order error terms in eqs. ((ref))--((ref)) as
The first two terms on the right--hand side of above display contain the estimation errors due to $\{\hat{\lambda}_{i}\}_{i=1}^{N}$ and $\{\hat{f} _{t}\}_{t=1}^{T}$, while the third term contains the error of the infeasible nonparametric estimator $\hat{g}_{0,Z}^{\ast }$; under weak regularity conditions, the third term will be negiglible:
The Donsker-type condition of Andrews1994 cited in above lemma to ensure stochastic equicontinuity is satisfied by many nonparametric regression estimators, including above kernel regression estimators. On the other hand, it is unclear whether the k-means clustering algorithm of beyhum2025inference satisfies a similar condition; in particular, their estimator is by construction non--smooth and so will not belong to the $ L_{2} $-Sobolev space of Andrews1994. As such, their estimator may not be covered by above lemma.
What remains is the analysis of the first two terms in eq. ((ref)). Consider the first component of the vector $\hat{\nu}(\hat{g} _{0,Z})-\hat{\nu}^{\ast }(\hat{g}_{0,Z})$, denoted $\hat{\nu}_{Y}\left( \hat{ g}_{0,X}\right) -\hat{\nu}_{Y}^{\ast }(\hat{g}_{0,X}):=\frac{1}{NT} \sum_{it}\varepsilon _{it}\{\hat{g}_{X}(\hat{\lambda}_{i},\hat{f}_{t})-\hat{g }_{0,X}(\lambda _{i},f_{t})\}$. A second--order Taylor expansion yields
where the remainder term $R_{N,T}$ in great generality will be of higher order than the other terms in the final expression. A similar expansion holds for the second component of $\hat{\nu}(\hat{g}_{0,Z})-\hat{\nu}^{\ast }(\hat{g}_{0,Z})$, $\hat{\nu}_{X}\left( \hat{g}_{0,Y}\right) -\hat{\nu} _{X}^{\ast }(\hat{g}_{0,Y}):=\frac{1}{NT}\sum_{it}\eta _{it}\{\hat{g}_{0,Y}( \hat{\lambda}_{i},\hat{f}_{t})-\hat{g}_{0,Y}(\lambda _{i},f_{t})\}$. The leading terms of these expansions will create incidental parameter biases if the first-order bias and/or variance components of $\hat{\lambda}_{i}$ or $ \hat{f}_{t}$ in ((ref)) correlate with $\varepsilon _{Z,it}$. The precise expressions of these bias components can be derived using arguments similar to the ones in FernandezValWeidner2016. At the same time, $\hat{\nu}(\hat{g}_{0,Z})-\hat{\nu}^{\ast }(\hat{g}_{0,Z})$ will generally not contribute to the variance of $\hat{\beta}$ due to the use of Neyman orthogonality of the moment conditions.
The analysis of the second term in eq. ((ref)) requires us to take a firmer stand on the nonparametric regression method being employed. For example, if kernel regression with a second--order kernel is employed, then we expect
for some functions $b_{g,1}\left( \lambda ,f\right) $ and $b_{g,2}\left( \lambda ,f\right) $. The expressions of these can be obtained by combining ( (ref)) with the arguments of Okui2019 and Barras2021. Given that $\varepsilon _{it}$ and $\eta _{it}$ do not correlate with $\lambda _{i}$ and $f_{t}$, these terms will not contribute to the incidental parameter biases in ((ref))--((ref)). We expect similar results to hold for other nonparametric regression estimators. In conclusion, the incidental parameter biases $ B_{\alpha ,k}/T+B_{\gamma ,k}/N$, $k=2,3$, in Theorem (ref) are expected to arrive from $\hat{\nu}(\hat{g}_{0,A})-\hat{\nu}^{\ast }(\hat{ g}_{0,A})$ if $\hat{\lambda}_{i}$ or $\hat{f}_{t}$ correlate with $ \varepsilon _{Z,it}$. Moreover, again due to the Neyman orthogonal moment conditions, the variance component of $\hat{\nu}(\hat{g}_{0,Z})-\hat{\nu} ^{\ast }(\hat{g}_{0,Z}^{\ast })$ will generally be asymptotically negiglible.
For the second--order term in ((ref)), we can carry out a similar decomposition. Assuming that $\hat{g}_{0,Z}^{\ast }$ is sufficiently regular, $\frac{1}{NT}\sum_{it}\left\Vert \hat{g}_{0,Z}^{\ast }(\lambda _{i},f_{t})-g_{0,Z}(\lambda _{i},f_{t})))\right\Vert ^{2}=o_{P}(1/\sqrt{NT})$ ; most known nonparametric regression estimators will satisfy this under standard regularity conditions. Moreover, we expect the second--order quadratic term will contain an incidental parameter bias terms on the form
c.f. FernandezValWeidner2016. We also expect $\hat{g}_{0,Z}\left( \lambda ,f\right) -\hat{g}_{0,Z}^{\ast }\left( \lambda ,f\right) $ to contribute; for example, for kernel--based estimators, we expect $\frac{1}{NT }\sum_{it}(\hat{g}_{0,Y}(\lambda _{i},f_{t})-\hat{g}_{0,Y}^{\ast }(\lambda _{i},f_{t}))(\hat{g}_{0,X}(\lambda _{i},f_{t})-\hat{g}_{0,X}^{\ast }(\lambda _{i},f_{t}))^{\prime }=O_{P}(1/N)+O_{P}(1/T)$ under weak conditions, c.f. Okui2019 and Barras2021.
One could now attempt to carry out a more complete characterisation of the incidental parameter biases and then use the resulting expressions of the biases to carry out bias adjustment of $\hat{\beta}_{NO}$ in Theorem (ref). We will refrain from carrying out such an analysis since this will require us to restrict ourselves to the ideal scenario of regular estimators satisfying ((ref)), and take a stand on the precise form of $\hat{g}_{0,Z}$.
As explained in the previous subsection, the incidental parameter biases will arise due to the use of $\hat{g}_{0,Z}(\hat{\lambda}_{i},\hat{f}_{t})$ in place of $\hat{g}_{0,Z}^{\ast }\left( \lambda _{i},f_{t}\right) $. We here explain how these biases can be adjusted for by employing a generalised version of the two-way kernel regression estimator of freeman2023linear,freeman2022multidimensional.
For a given nonparametric regression technique, let $\hat{g}_{0,Z}(\lambda ,f)$ be the full--sample version that regresses $Z_{js}$ onto $(\hat{\lambda} _{j},\hat{f}_{s})$ for $j=1,...,N$ and $s=1,...,T$. Next, for each $ t=1,...,T $, let $\hat{g}_{0,Z}^{\left( 1\right) }(\lambda ,f_{t})$ be the estimator obtained by regressing $Z_{jt}=g_{0,Z}(\lambda _{j},f_{t})+\varepsilon _{Z,it}$ onto $\hat{\lambda}_{j}$ for $j=1,...,N$; that is, we run $T$ nonparametric regressions, each along the cross-sectional dimension. Finally, for each $i=1,...,N$, let $\hat{g} _{0,Z}^{\left( 2\right) }(\lambda _{i},f)$ be the estimator obtained by regressing $Z_{is}$ onto $\hat{f}_{s}$ for $s=1,...,T$; that is, we run $N$ nonparametric regressions, each along the time dimension. We then combine these to obtain the following two--way nonparametric regression estimator,
When the kernel regression estimator in ((ref)) is employed, it takes the form
The two--way estimator has the following two attractive features: First, if ( (ref)) holds, then under great generality,
Thus, $\mathbb{E}[\hat{\Gamma}_{0,Z,it}^{TW}|\hat{\lambda}_{i},\hat{f} _{t}]=o\left( 1/T\right) $ and so the leading biases due to $\{\hat{\lambda} _{j},\hat{f}_{s}\}_{js}$ being used in the nonparametric regression has been removed. One can think this of this as a type of Jackknifing. Second, assuming the first order derivatives of $\hat{g}_{0,Z}$ and $\hat{\lambda} _{i},\hat{f}_{t}$ are consistent,
and similarly for the second-order derivatives. As such, $\hat{\Gamma} _{0,Z,it}^{TW}$ is Neyman orthogonal to $(\hat{\lambda}_{i},\hat{f}_{t})$. As a consequence, if we carry out the expansions leading to eqs. ((ref)) and ((ref)) with $\hat{g}_{0,Z}(\hat{\lambda}_{i},\hat{f} _{t})$ replaced by $\hat{\Gamma}_{0,Z,it}^{TW}$, all the partial derivatives in the final expressions are now $o_{P}\left( 1\right) $ and so the leading biases due to $(\hat{\lambda}_{i},\hat{f}_{t})$ being used as arguments in $ \hat{g}_{0,Z}(\lambda ,f)$ also become asymptotically negiglible. All together, we expect $\hat{\beta}^{NO}$ that takes $\hat{\Gamma} _{0,Z,it}^{TW} $ as input not to suffer from any incidental parameter biases.
The two--way estimator could also be employed in other settings. For example, deaner2025inferring feeds the following sample pseudo--metrics of zhang2017estimating,
into a kernel regression estimator. Their kernel regression estimator could then be combined with above two--way estimator to remove parts of the estimation error due to the "generated" kernel regressors ${\hat{d}} _{Y,ij}^{(1)}$ and ${\hat{d}}_{Y,st}^{(2)}$. Simulation evidence strongly suggests that this indeed leads to substantial improvements in the resulting regression estimator. This metric is, however, sensitive to the dimension of $\alpha _{i}$ and $\gamma _{t}$.\footnote{ See the results of our simulation study in Section (ref).}
In order to formalise the arguments in the past two subsections and obtain precise expressions of the incidental parameter biases, we would need to impose restrictions on the precise form of $(\hat{\lambda}_{i},\hat{f}_{t})$ and $\hat{g}_{0,Z}$. We here develop sample--splitting procedures that avoids us having to do so: First, sample splitting removes the leading biases caused by the first-order terms under very weak conditions on the chosen nonparametric regression method. Second, sample--splitting allows us to show that the two--way regression method removes the leading biases of the second--order term under weak conditions on $(\hat{\lambda}_{i},\hat{f} _{t})$ and $\hat{g}_{0,Z}$.
We develop two sample--splitting procedures. The first is standard in the literature and works for any, possibly non--linear and possibly non--regression--based estimators of $\Gamma _{Z,it}$; this version is able to remove the incidental parameter biases incurred from the two first-order terms in ((ref))--((ref)), but the biases arising from the second--order term in ((ref)) may still be present. The second version appears to be new to the panel data literature and applies to linear estimators of $g_{0,Z}$; this version will be able to remove the biases from both the linear and second-order terms when combined with the two--way estimator.
We follow freeman2023linear, amongst others, and here consider the following sample splitting procedure: Define the following sample splits in terms of the indices of the data points,
and $\mathcal{I=}\left\{ i=1,....,N,t=1,....,T\right\} $. We then compute our first--step estimators so it does not use data from $\mathcal{I} _{k_{1},k_{2}}$,
for $k_{1},k_{2}=1,2$, and use these to obtain
We can still apply the same arguments that lead to Theorem (ref) to $\hat{\beta}_{NO}^{SS}$ and so the rate conditions stated in ( (ref))--((ref)) still hold, except that in the left-hand side expressions $\sum_{it}$ and $\hat{\Gamma}_{Z,it}$ are now replaced by $\sum_{\left( i,t\right) \in \mathcal{I}_{k_{1},k_{2}}}$ and $ \hat{\Gamma}_{Z,it}^{\left( k_{1},k_{2}\right) }$, respectively. Eqs. ((ref))--((ref)) will hold under the following weak restrictions on the variances of the two error terms, where $ \varepsilon ^{\left( k_{1},k_{2}\right) }=\left\{ \varepsilon _{it}:\left( i,t\right) \in \mathcal{I}_{k_{1},k_{2}}\right\} $ and $\eta ^{\left( k_{1},k_{2}\right) }=\left\{ \eta _{it}:\left( i,t\right) \in \mathcal{I} _{k_{1},k_{2}}\right\} $:
This generalises Theorem 1 in beyhum2025inference to allow for a very broad class of first--step estimators $(\hat{\Gamma}_{X,it},\hat{\Gamma} _{Y,it})$ and to allow for possible time series and cross-sectional dependence. In particular, we expect that ((ref)) will hold under suitable weak dependence conditions as explored in lunde2019sample.
The theorem shows that sample splitting alone allows removes incidental parameter biases from the two first-order terms in great generality. However, the second--order term in ((ref)) may still generate incidental parameter biases.
We here develop an alternative sample splitting procedure that not only removes incidental parameter biases arising from the first-order terms but also from the second--order term when combined with the two--way estimator where $\hat{g}_{0,Z}^{\left( 1\right) }(\hat{\lambda}_{i},f_{t})$, $\hat{g} _{0,Z}^{\left( 2\right) }(\lambda _{i},\hat{f}_{t})$, and $\hat{g}_{0,Z}( \hat{\lambda}_{i},\hat{f}_{t})$ are obtained using a linear nonparametric regression procedure. One example is ((ref)), but other options are also possible, such as a two--way series estimator.
Take a linear smoother $\mathcal{W}_{it}^{\left( k_{1},k_{2}\right) }$ computed from data in $\mathcal{I}\backslash \mathcal{I}_{k_{1},k_{2}}$ only,
and let
denote data from $\mathcal{I}_{k_{1},k_{2}}$. We then define
where $\hat{\Gamma}_{Z}$ simply concatenates the partitioned matrices of estimates. This sample split ensures the smoother $\mathcal{W}_{it}^{\left( k_{1},k_{2}\right) }$ has no stochastic dependence with $\eta _{it}$ and $ \varepsilon _{it}$ for $\left( i,t\right) \in \mathcal{I}_{k_{1},k_{2}}$ under independence across $\left( i,t\right) $. Note that whilst the freeman2023linear sample split in (ref) is sufficient for this, with mutually independent $\eta _{it}$ and $\varepsilon _{it}$ a leave-one-out split would also be sufficient but requires $N\times T$ estimations, i.e. for each $i$ and $t$, so we do not implement this version.
The sample splitting estimator still takes the form ((ref)) and so we could again apply Theorem (ref) as we did in the previous subsection. However, when linear estimators are employed in the first-stage, a more precise expansion can be obtained that leads to sharper rate restrictions on the first--stage. For any two matrices $A,B\in \mathbb{R}^{N\times T}$, let $\langle A,B\rangle _{F}:=\sum_{it}A_{it}B_{it}/\left( NT\right) $ denote the scaled entrywise Frobenius inner product. With $\hat{\Gamma}=\mathcal{W}(\Gamma +\varepsilon )\in \mathbb{R}^{N\times T}$, where $\Gamma _{it}:=g(\alpha _{i},\gamma _{t}) $, we have $\hat{\Gamma}_{Y}=\hat{\Gamma}_{X}\beta +\hat{\Gamma}$ and so for $\hat\Omega_X := (NT)^{-1} \sum_{it} (X_{it}-\hat{\Gamma} _{X,it})(X_{it}-\hat{\Gamma}_{X,it})^\prime$,
Hence, for the estimator to be asymptotically normally distributed without any asymptotic biases, we need ((ref)) to hold together with\
Importantly, compared to the general case, $\hat{\Gamma}_{Y}-\Gamma _{Y}$ has been replaced by $\hat{\Gamma}-\Gamma $. Sample splitting will now ensure that $\langle \Gamma _{X}-\hat{\Gamma}_{X},\varepsilon \rangle _{F}=o_{p}(1/\sqrt{NT})$ and $\langle \Gamma _{X}-\hat{\Gamma} _{X},\varepsilon \rangle _{F}=o_{p}(1/\sqrt{NT})$. What remains is to show ( (ref)). We conclude:
We establish $\Vert \mathcal{W}^{\ast }\mathcal{W}\eta_k \Vert _{F}^{2} =o_p(NT)$, and $\mathbb{E}[\langle\Gamma - \mathcal{W}\Gamma,\Gamma_X - \mathcal{W}\Gamma_X\rangle_F] = o_p(NT)^{-1/2}$ for our specific estimators in Section (ref). Condition $\langle\mathcal{W}\mathbb{E} [\varepsilon|\mathcal{F}], \mathcal{W}\mathbb{E}[\eta|\mathcal{F}] \rangle _{F} = o_p(NT)^{-1/2}$ holds in many weak dependence settings, see Remark (ref). It trivially holds for iid $\varepsilon$, since then $ \mathbb{E}[\varepsilon|\mathcal{F}] = 0$.
Condition $\Vert \mathcal{W}^{\ast }\mathcal{W}\eta_k \Vert _{F}^{2} =o_p(NT) $ is satisfied for many $\mathcal{W}$. Take a cross-sectional smoother for simplicity, such that $\mathcal{W} = w\in\mathbb{R}^{N\times N}$ . When weights are set to $1/N$, we get the usual $\mathcal{W}\eta_k = N^{-1} \iota_N \cdot \iota_N^\prime \eta_k = \iota_N \iota_T^\prime \cdot O_p(1/\sqrt{N})$. When nonparametric smoothers are used with bandwidth $h$, $ \mathcal{W}\eta_k = \iota_N \iota_T^\prime \cdot O_p(1/\sqrt{Nh})$ under common regularity conditions. When additive nonparametric smoothers are used, with bandwidth $h$ and $R$ additive terms, this becomes $\mathcal{W} \eta_k = \iota_N \iota_T^\prime \cdot O_p(R/\sqrt{Nh})$. These rates get slower with more flexible smoothers, but are easily all $o(1)$.\footnote{ The condition for additive smoother in this example is $R^2/h = o(N)$. Since we only consider $R \to \infty$ very slowly as a function of $\min\{N,T\}$, this is not difficult to satisfy.}
In this section, we import a novel identification result for $g_{Z}$ and the fixed effects from the companion paper FreemanKristensen2026 and use this to develop a kernel regression-based estimators of these. We then proceed to verify that this estimator satisfies the high--level conditions in Theorem (ref).
We first introduce a singular value decomposition (SVD) of $g$: Let, for a given multi-index $\iota =\left( \iota _{1,1},\dots ,\iota _{1,d_{\alpha }},\iota _{2,1},\dots ,\iota _{2,d_{\gamma }}\right) \in \mathbb{N} _{0}^{d_{\alpha }+d_{\gamma }}$,
be the mixed partial derivative. Let the norm $\Vert g\Vert _{L_{f}^{2}(\Omega _{\alpha }\times \Omega _{\gamma })}$ denote the usual $ L_{2}$-norm over support of $\left( \alpha _{i},\gamma _{t}\right) $, denoted $\Omega _{\alpha }\times \Omega _{\gamma }$ taken with respect to joint distribution of $\left( \alpha _{i},\gamma _{t}\right) $,
where $\pi _{\alpha }$ and $\pi _{\gamma }$ denote the densities of $\alpha _{i}$ and $\gamma _{t}$. We will then assume that:
Under this assumption, we obtain the following SVD of $g$,
where $\sigma _{1}>\sigma _{2}>....$ are the ordered singular values, and $ u_{r}\in H_{f}^{p}(\Omega _{\alpha }\times \Omega _{\gamma })$ and $v_{r}\in H_{f}^{p}(\Omega _{\alpha }\times \Omega _{\gamma })$ are eigenfunctions. These constitute an orthonormal basis of $L_{f}^{2}(\Omega _{\alpha })$ and $ L_{f}^{2}(\Omega _{\gamma })$, respectively, so that
We refer to griebel2014approximation and freeman2023linear for further details. Importantly, the eigenfunctions are identified, c.f. freeman2023linear.
Next, we impose the following injectivity condition on $g$:
In discussing Assumption (ref), it is important to note that there exist many observational equivalent representations of $g\left( \alpha _{i},\gamma _{t}\right) $. Suppose, for example, that $g\left( \alpha _{i},\gamma _{t}\right) =\tilde{g}\left( A\alpha _{i},B\gamma _{t}\right) $, where $A\in \mathbb{R}^{\tilde{d}_{\alpha }\times d_{\alpha }}$, $B\in \mathbb{R}^{\tilde{d}_{\gamma }\times d_{\gamma }}$ and $\tilde{g}_{Z}: \mathbb{R}^{\tilde{d}_{\alpha }}\times \mathbb{R}^{\tilde{d}_{\gamma }}$ with $\tilde{d}_{\alpha }<d_{\alpha }$ and $\tilde{d}_{\gamma }<d_{\gamma }$ . In this case, a lower--dimensional observational equivalent representation is $\tilde{g}\left( \tilde{\alpha}_{i},\tilde{\gamma}_{t}\right) $, where $ \tilde{\alpha}_{i}=A\alpha _{i}$ and $\tilde{\gamma}_{t}=B\gamma _{t}$. More generally, if, for any given tuple $(\alpha _{1},...,\alpha _{d_{\gamma }})$ , the mapping $\gamma \mapsto \left( g\left( \alpha _{1},\gamma \right) ,....,g\left( \alpha _{d_{\gamma }},\gamma \right) \right) $ does not have full rank or, for any given tuple $(\gamma _{1},...,\gamma _{d_{\alpha }})\in R^{d_{\gamma }\times d_{\gamma }}$, the mapping \newline $\alpha \mapsto \left( g\left( \alpha ,\gamma _{1}\right) ,....,g\left( \alpha ,\gamma _{d_{\alpha }}\right) \right) $ does not have full rank, we say that $g\left( \alpha _{i},\gamma _{t}\right) $ is reducible. In either case the finite--dimensional distribution of $\left\{ Z_{it}\right\} _{1\leq i\leq d_{\gamma },1\leq t\leq d_{\alpha }}$ can be represented by some lower--dimensional mapping $\tilde{g}\left( \tilde{\alpha}_{i},\tilde{\gamma} _{t}\right) $, where $\tilde{d}_{\alpha }<d_{\alpha }$ or $\tilde{d}_{\gamma }<d_{\gamma }$.
In the light of these observations, Assumption (ref) is quite weak. If Assumption (ref) is violated, and so $g\left( \alpha _{i},\gamma _{t}\right) $ is reducible, then we can work with the lower--dimensional observational equivalent representation of $g$ that will satisfy above assumption. As such, we find that the assumptions (ref) imposes very weak restrictions on $g\left( \alpha _{i},\gamma _{t}\right) $.
The following theorem shows that under Assumption (ref) a linear combination of a finite number of the leading eigenfunctions are valid proxies. Moreover, the associated regression function inherits the smoothness properties of its mother.
Estimation of $g_{0,Z}$ is a high--dimensional nonparametric regression problem if $p_{\alpha }$ and/or $p_{\gamma }$ are large, and so may suffer from the well-known curse--of--dimensionality which leads to large errors in the nonparametric estimation. We therefore now introduce restrictions under which this curse is less of a concern. To simplify notation, we here assume that $d_{\alpha }=d_{\gamma }=d$.
Next, we develop two--step regression estimators of $g_{Z}$ based on above identification result: In the first step, we estimate the leading eigenfunctions of $g$ that in the second step are used as proxies for $ \alpha _{i}$ and $\gamma _{t}$ in a nonparametric regression procedure.
This section presents our first-step estimators of the leading eigenfunctions of the SVD\ representation of $g$ in ((ref)). Substituting ((ref)) into ((ref)) and truncating the singular value decomposition at some $ R_{1}\geq 1$ chosen by the econometrician yields
where $\lambda _{ir}=\sigma _{r}u_{r}(\alpha _{i})$, $f_{tr}=v_{r}(\gamma _{t})$, and $e_{R_{1},it}=\sum_{r=R_{1}+1}^{\infty }\lambda _{ir}f_{tr}$ is the truncation error. We follow freeman2023linear and obtain first-step estimates of $\lambda \in \mathbb{R}^{N\times R_{1}}$ and $f\in \mathbb{R}^{T\times R_{2}}$ by applying the estimator of Bai2009 to above approximate factor model,
where we impose the normalisations from Bai2009,bai2023approximate, $ N^{-1}\sum_{i}\lambda _{1:R_{1},i}\lambda _{1:R_{1},i}^{\prime }$ is diagonal and $T^{-1}\sum_{t}f_{1:R_{1},t}f_{1:R_{1},t}^{\prime }=\mathbb{I} _{R_{1}}$.
Above algorithm delivers estimators of the leading $R_{1}$ eigenfunctions, $ \hat{\lambda}_{1:R_{1}}$ and $\hat{f}_{1:R_{1}}$. Importantly, compared to the alternative estimation procedure of beyhum2025inference, the algorithm effectively reduces the dimension of the fixed effects to be controlled for since it takes into account the presence of $\beta ^{\prime }X_{it}$ in the model. If the DGP for $X_{it}$ takes the form ((ref)), then the algorithm of beyhum2025inference will generally take into account not only the fixed effects $(\alpha _{i},\gamma _{t})$ that enter the model for $Y_{it}$ but also the fixed effects $(\alpha _{i}^{\left( 2\right) },\gamma _{t}^{\left( 2\right) })$ that are specific to $X_{it}$. In contrast, above algorithm does not suffer from such shortcomings since it controls for $\beta ^{\prime }X_{it}$ and so directly targets $(\alpha _{i},\gamma _{t})$.
The algorithm also delivers estimators $\hat{\beta}_{LS}$ and $\hat{\Gamma} _{it}=\sum_{r=1}^{R_{1}}\hat{\lambda}_{ir}\hat{f}_{tr}$ of $\beta $ and $ \Gamma _{it}$, respectively. However, these estimators suffer from large errors due to the truncation error $e_{R,it}$ and so $\hat{\beta}_{LS}$ will not enjoy $\sqrt{NT}$-asymptotic normality: The factor model approach to Neyman Orthogonal estimator uses $\Gamma -\hat{\Gamma}=(\mathbb{I}-P_{\hat{ \lambda}})g(\alpha ,\gamma )(\mathbb{I}-P_{\hat{f}})$, which leads to,
c.f. bai2023approximate and Section (ref). The first term $O_{p}(R_{1}^{1-\rho })$ is due to the truncation error and shrinks as $R_{1}$ grows, and is also decreasing in smoothness $\rho $, i.e. smoother functions lead to smaller bias. Variance term $O_{p}(R_{1}^{2+2\rho }\min \{N,T\}^{-1})$, however, is increasing in model complexity $R_{1}$, and also increasing in smoothness $\rho $. Setting $R_{1}=c\cdot \min \{N,T\}^{\frac{1}{1+3\rho }}$ to balance bias and variance leads to $\frac{1 }{\sqrt{NT}}\Vert \Gamma -\hat{\Gamma}\Vert =\min \{N,T\}^{\frac{1-\rho }{ 1+3\rho }}$. Hence, we can at best obtain $||\hat{\beta}_{LS}-\beta ||=O_{P}\left( \min \{N,T\}^{\frac{2-2\rho }{1+3\rho }}\right) $ which is too slow for $\sqrt{NT}$-inference.
Figure (ref) shows the two limiting components to the first-step estimator. Either singular values decay too slowly, and this leaves a large bias from $e_{it}$, or singular values decay too fast and variation from $\sum_{r=1}^{R_{1}}\lambda _{ir}f_{tr}$ is indiscernible from $\varepsilon _{it}$, leading to higher variance. Figure (ref) shows a signal to noise comparison for the distribution of singular values generated from a function and from Gaussian noise. When signal drops below noise, factors from the function are no longer estimable, or estimated with noise. We see in the left panel that for non-smooth functions with slowly decaying singular values, we can potentially estimate and control for many factors, however, there is a large error that persists in the tail of the approximation. This leads to large bias. In the right panel, whilst the approximation error in the tail is small, variation from the function quickly becomes indiscernible from noise, hence estimates are noisy. As as a consequence, in either scenario, the over--all estimation error of $\hat{g}$ is too large and do not vanish at the rate required in Theorem (ref).
We here develop two nonparametric regression algorithms that both take as input the subset of the first $R_{2}$ of the $\ R_{1}$ estimated leading eigenfunctions in the first step, where again $R_{2}$ is chosen by the econometrician. With some abuse of notation, we let $\hat{\lambda}_{i}=(\hat{ \lambda}_{1,i},.....,\hat{\lambda}_{R_{2},i})^{\prime }$ and $\hat{f}_{t}=( \hat{f}_{1,t},.....,\hat{f}_{R_{2},t})^{\prime }$ denote these final $R_{2}$ estimated eigenfunctions. We will require $R_{2}\geq R_{0}$ so that we can apply Theorems (ref) and (ref) and obtain consistent estimators based on the representation results in equations ((ref)) and ((ref)), respectively.
The following kernel regression estimator is a consistent estimator of $ g_{0,Z,k}\left( \lambda ,f\right) $ defined in (ref),
where $K_{\hat{A}_{k},h_{1}}(\hat{\lambda}_{i}-\lambda ):=K_{1}(\hat{A}_{k}\{ \hat{\lambda}_{i}-\hat{\lambda}_{i_{0}}\}/h_{1})/h_{1}^{d_{\lambda }}$, $K_{ \hat{B}_{k},h_{2}}(\hat{f}_{t}-f)=K_{2}(\hat{B}_{k}\{\hat{f} _{t}-f\}/h_{2})/h_{2}^{d_{f}}$, $h_{1}$ and $h_{2}$ are bandwidths, $K_{1}$ and $K_{2}$ are kernels, and
Above estimator is a so--called multi--index regression estimator with generated regressors $\hat{\lambda}_{i}$ and $\hat{f}_{t}$. When the regressors are observed without errors, this has estimator has been analyzed in, among others, Ma2012 and Ma2013. This estimator will suffer from a curse-of--dimensionality of order $\max \left\{ p_{\alpha },p_{\gamma }\right\} $.
Theorem (ref) allows us to use additive nonparametric estimation techniques to remove the curse--of--dimensionality that above multi--index estimator suffers from. We here propose to employ the kernel--based backfitting algorithm for nonparametric additive models to remove this curse of dimensionality. We refer to opsomer1999root for details on this algorithm in a cross--sectional setting. We here extend it to a panel data setting.
Figure (ref) shows a comparison of the linear projection method from a factor model versus a nonparametric difference, where true $ g(\alpha_i,\gamma_t)$ is observed. We see that whilst error for the first three directly estimated terms is zero for linear projection methods, and positive for nonparametric methods, there is a huge benefit in the tail of residual singular values. By definition, factor components are exactly orthogonal to tail eigenfunctions, however the nonparametric estimator can still difference these out as a result of the additive model. In contrast, nonparametric methods allow for infinite linear dimension, albeit with some additional smoothness conditions.
Here we combine the first-step estimators with a backfitting algorithm that iterates over weighted-within transformations from freeman2022multidimensional to obtain our final estimates. The weighted-within transformation performs dimension specific differencing of the fixed-effects. Take the set of estimates $\{\hat{\lambda}_{ir},\hat{f} _{tr}\}_{r=1}^{R_{1}}$. Smoother weights are formed for the $i$, respectively $t$ direction as,
Dependent and independent variables are sequentially residualised with respect to weighted-differences according to the sequence of weights for $ r=1,\dots ,R_{2}$. For our asymptotic theory, we require a backfitting update. With $\check{Y}^{(0)}=Y-\bar{Y}$, step $\ell $ in the backfitting iteration can be written,
The algorithm works as follows. Take a generic smoothing function $s:\mathbb{ R}^{n}\times \mathbb{R}^{n}\rightarrow \mathbb{R}^{n\times n}$, and $\hat{ \lambda}_{ir}$ and $\hat{f}_{tr}$, $r=1,...,R_2$, We use backfitting iteration from Algorithm (ref):
If convergence occurs at $L$ iterations, the procedure can be written,
We define the estimator under these weighted differences as $\hat{\beta}_{W}$ :
In this section we formally analyse asymptotic rates for the two--step estimators of $g_{Z}$.
freeman2023linear derive convergence rates of the least squares estimator in Section (ref). These, however, can be improved by using refined estimates from griebel2019singular, which gives faster decay in the singular values for the tail term, $\sum_{r=R_{1}+1}^{\infty }\lambda _{ir}f_{tr}$. Let $\rho :=p/\min \{d_{\alpha },d_{\gamma }\}$, where $p$ is introduced in (ref), and $\{d_{\alpha },d_{\gamma }\}$ are dimensions of $\alpha $, and $\gamma $. This comprises Lemma (ref).
Assumption (ref) refines Lemma 1 in freeman2023linear by applying Corollary 3.4 and Proposition 3.5 from griebel2019singular. Here we state regularity conditions from freeman2023linear,MoonWeidner2015 on covariates and $\varepsilon$.
Condition $\| \mathbb{E}[ee^{\prime }|X]\|_2 = O_p(1)$ bounds cross and serial correlation. For example, i.i.d. $\varepsilon_{it}$ implies $\mathbb{E }[ee^{\prime }|X] = \mathbb{E}\varepsilon_{it}^2 \mathbb{I}_{NT}$ such that $ \| \mathbb{E}[ee^{\prime }|X]\|_2 = \mathbb{E}\varepsilon_{it}^2 = O(1)$. With weak correlation in $i$ and $t$ such that $\sum_{j\neq i}\sum_{s\neq t}| \mathbb{E}[\varepsilon_{it}\varepsilon_{js}|X]| = O(1)$, then Ostrowski theorem implies $\| \mathbb{E}[ee^{\prime }|X]\|_2 = O_p(1)$.
We first analyse the first-step estimators in ((ref)):
For $N\sim T$, this is at best $O_{p}(NT)^{-1/4}$ when $\rho \rightarrow \infty $; this is too slow for standard inference tools to be valid.
In Appendix (ref) we state regularity conditions on $\{\eta _{it},\varepsilon _{it}\}$ and their correlation with $\{\lambda _{i},f_{t}\} $ to apply results from Bai2009 and bai2023approximate for factor model estimates. Results in choi2025high may also apply to this setting. Define,
where $H$ are rotation matrices defined in the appendix, see bai2023approximate.\footnote{ These are not important to our analysis so we refer to discussion found in existing literature.}
In Corollary (ref) we show $R_2^{2}\cdot \xi _{f}^{2}\xi _{\lambda }^{2}=o_{p}(\min \{N,T\}^{-1})$ and $\frac{R_2^{6}}{NT}\cdot \xi _{f}^{2}\xi _{\lambda }^{2}=o_{p}(\min \{N,T\}^{-2})$, which is sufficient for our result stated in Corollary (ref) below. Since $\hat{\beta}$ is arbitrary in the statement of Lemma (ref), we can set a different number of factors for $\beta $ estimate and factor component estimates.
It may be possible to attain similar, or potentially faster rates, by iterating between factor estimation and our final $\hat\beta_{NO}$ estimator. However, the above rate is sufficient, and easy to verify, so we do not confirm this.
We here only analyze the additive version of the two proposed kernel regression estimators of $g_{0,Z}$. Since both the additive and the multi-index version are kernel regressions, our analysis of the additive version will carry over the multi-index estimator with obvious modifications. To handle the estimation of the multi-index parameters, we can apply the techniques developed in Ma2012 and Ma2013.
We add some regularity to the kernel functions and distribution of eigenfunctions.
Assumption (ref) are standard restrictions. We regularise the distribution of proxies.
Assumptions (ref), (ref) are standard nonparametric estimation assumptions, see Assumptions 1 and $2^{\prime }$-$4^{\prime }$ from opsomer1999root. Assumptions (ref) can be justified with estimated eigenfunctions of functions that adhere to Assumption (ref), which pointwise converge to the true eigenfunctions, which are in a compact set.
We state convergence results for $\widehat\beta_W^{SS}$, the estimator from Section (ref) using the sample splitting described in Section (ref) to estimate smoothing weight matrices $S^{(1)}$ and $S^{(2)}$. Lemma (ref) applies results from opsomer1999root for the semiparametric additive model. Define $\xi_{\Gamma X}$ as:
We verify $RMSE(\xi_{\Gamma X}) = o_p(NT)^{-1/2}$, which is sufficient for Theorem (ref). Assume without loss that $X_{it}$ is scalar, and expand ${\Gamma }_{X}-\widehat{\Gamma }_{X} = \widetilde{\Gamma }_{X} - \hat\eta$, where
Define $\widetilde{\Gamma } $ and $\hat\varepsilon$ similarly. In turn the definition of $\xi_{\Gamma X}$ implies,
To show any variance added to the estimates of $\hat\beta_{NO}^{SS}$ from $ \hat\eta$ and $\hat\varepsilon$ is sufficiently small order, i.e. $ o_p(NT)^{-1/2}$, it is convenient to analyse bias and variance of each term in this expression directly.
In the following Lemma we show the results hold under mutually independent proxies, which leads to much simpler theoretical properties. We present Lemma (ref) as it may be of practical interest.
Next we state a bound on the variance for the multiplicative error:
Finally, we use the rate results derived in the previous subsection to verify the general conditions of Theorem (ref) when our novel estimator of $g_Z$ is employed:
In Corollary (ref), we show for $R_2 = O(\min\{N,T\}^{\frac{1}{ 6\rho}})$ that $R_2^2\cdot\xi_\lambda^2\xi_f^2 = o_p(\min\{N,T\}^{-1})$, and $R_2^6 (NT)^{-1}\xi_\lambda^2\xi_f^2 = o_p(\min\{N,T\}^{-2})$ for $N\sim T$. For the range of $h_\lambda$, and $h_f$ stated in Corollary (ref) all other terms in bias are $o_p(\min\{N,T\}^{-1})$, and in variance are $o_p(\min\{N,T\}^{-2})$.
Recall in (ref) that $\Omega_X := \text{plim} (NT)^{-1} \sum_{it} \eta_{it}\eta_{it}^{\prime }$, and $\Sigma$ from $(NT)^{-1/2} \sum_{it} \eta_{it}\varepsilon_{it} \xrightarrow[]{d}N(0,\Sigma)$.
Estimation of the preliminary model, along with final estimation of the multi step debiasing estimators involved the following hyperparameters: $ \{R_1, R_2, h_\lambda, h_f\}$.
We motivate that $R_2$, i.e. the number of eigenfunctions used in the final estimation step, should be relatively small with respect to sample size. In practice, we advocate these being moderate, but fixed, in large enough samples. For example, if $\alpha,\gamma \in \mathbb{R}^d$ then $R_2 \geq d$ should in most cases control all latent fixed-effect variation. Hence, as long as eigenfunctions estimated in the first stage $R_1 \geq R_2 \geq d$, then our model should identify all latent variation in $g(\alpha, \gamma)$.
Our first-step estimation of eigenfunction proxies, however, showed a clear bias variance trade-off in $R_1$, regardless of the ambient dimension of $ \alpha,\gamma \in \mathbb{R}^d$. This is because variation in the tail terms of $g(\alpha, \gamma) = \sum_{r =1}^\infty \sigma_r u_r(\alpha) v_r(\gamma)$ may still influence estimation of leading terms, even as $\sigma_r \to 0$. We conjecture here, without formal proof, that tests of a matrix rank observed with noise should estimate the $R_1$ that optimally trades-off bias and variance in this setting. That is, since in finite samples we would optimally set $R_1$ equal to the highest $r$ such that $\sigma_r > \left\| \varepsilon\right\|_2$, information criterion tests from e.g. BaiNg2002 or eigenvalue ratio test in ahn2013eigenvalue and Section 5 of ke2024robust for linear regression models specifically, should work well. These tests are naturally designed to choose the point at which the distribution of eigenvalues are related to noise, hence should do well to establish that optimal $R_1$ for estimation.
Lastly, for bandwidths $\{h_\lambda, h_f\}$ we again conjecture without proof that usual split sample cross-validation arguments should pick the optimal bandwidths. Since bandwidths from Section (ref) required to perform inference are the bandwidths that optimise mean squared error, the objective function of these cross-validation hyperparameter optimisers are aligned with our purposes. We leave formalising this for future work.
Numerical implementations in simulation Section (ref) and Appendix (ref) simply use the asymptotically optimal rates for $ \{R_1, R_2, h_\lambda, h_f\}$ supposed by the theory. We make no strong claim to estimating optimal hyperparameters, just that there exist some that conform to these rates which produce good numerical results.
Lastly, in finite samples, there are degrees of freedom corrections to consider. For the factor model, we use freeman2023linear adjustments to rescale variance estimators by,
Likewise, for the nonparametric estimators,
where $tr\{W^{(1)}\}$ and $tr\{W^{(2)}\}$ are effective degrees of freedom from nonparametric estimation, see hastie2009introduction Section 7.5.2, and $W^{(1)}, W^{(2)}$ are the nonparametric kernel weights derived in Section (ref) and (ref). We also implement the correction akin to rueda2013degrees:
This correction produces wider confidence intervals, with better coverage under dependence structures in the noise term. We make no formal claim to the validity of these corrections, but note in simulations that calculated standard errors approximate simulated standard deviations well.
We implement the weighted-within by using the multi-index weights in Section (ref) to initialise estimates before using the backfitting in Section (ref). Estimates under multi-index weights in Section (ref) perform well with nominal coverage, but are dominated by the backfitting in Section (ref), so we report only those results. For the higher dimensional simulations in Appendix (ref) we also find it helps to initialise the backfitting with the standard factor model.
Data is generated according to
where $\beta = 2$, $\alpha_i, \gamma_t \sim U(-1,1)$. For $\mu_{it}\sim N(0,1), \mu_{it}^* \sim N(0,4)$, $\nu_{it} \sim N(0,\eta_{it}^2)$, noise terms $\varepsilon _{it}$, and $\eta_{it}$ are generated,
then both normalised to have variance 1. Term $\mu_{it}^*$ normalises $ \eta_{it}$ to not have too much correlation. In this way they admit conditional heteroskedasticity and correlation over $i$, and $t$. To simulate that the cross-sectional ordering is unknown to the econometrician, we randomise the order of $i$, but maintain the order of $t$. Finally,
Terms $g(\alpha _{i},\gamma _{t})$ and $g_{X}(\alpha _{i},\gamma _{t})$ are normalised to variance four.\footnote{ Rescaling variance of these terms does not impact asymptotic results, but elucidates differences across estimators with much smaller sample sizes in simulations.} Figure (ref) displays the results graphically, Table (ref) tables the bias, root mean squared error, and coverage for 95% nominal test.
Standard errors across all estimators use the partial sum estimator from BaiNg2006 in conjunction with the newey1987hac kernel estimator in the time dimension.
Tuning parameters are chosen as follows. For the factor model estimators, denote $R_1=\min \{N,T\}^{1/3}$. For the nonparametric weighted-within estimator these are set to $R_2= 4$, i.e. constant, since the optimal rate $ R_2 = \min\{N,T\}^{1/6\rho}$ implies such slow rate of growth for reasonably smooth functions, that in practice optimal $R_2$ would not change over the sample sizes we consider. Notice here that $R_2 < R_1$, which follows our theoretical arguments: we find for first-stage preliminary estimators the bias and variance is better traded off with higher number of estimated components than in the second stage. Bandwidths for the nonparametric weighted-within estimator are set to $h=(25/\min \{N,T\})^{1/2}$ and for the Oracle estimator are set to $h=(1/\min \{N,T\})^{1/2}$. Both weighted-within and Oracle estimators use the Gaussian kernel.
In addition to the estimators compared in Table (ref) and Figure (ref) we also implement our estimator with the zhang2017estimating pseudo-metric from (ref), and also with cross-sectional/time-serial first moments, e.g. those used in bonhomme2021discretizing. We report these results in Appendix (ref). As conjectured, smoothers using the zhang2017estimating pseudo-metric performs well when $\alpha _{i}$ and $ \gamma _{t}$ are scalars, but scales poorly as their dimension increases. \footnote{ The zhang2017estimating pseudo-metric is computationally very costly, so we restrict Monte Carlo rounds to 1,000 and limit sample size to at most $ N=T=300$.} Our eigenfunction based method performs well for higher-dimensional $\alpha _{i}$ and $\gamma _{t}$, but does require smoother functions. This is predicted in our theory, where error decreases in $\rho =p/\min \{d_{\alpha },d_{\gamma }\}$. The moment based estimator performs poorly regardless of dimension for $\alpha _{i}$ and $\gamma _{t}$, since no cross-sectional or time-serial moments are injective for radial type functions proposed in (ref).\footnote{ In Appendix (ref) we implement DGPs that produce nominal coverage for the moment based estimator, and compare performance in those settings.}
In this paper we present novel theory for the linear panel setting with a general function specification of unobserved heterogeneity. Using the double-debias approach proposed in freeman2022multidimensional, deriving from chernozhukov2022locally methods, along with preliminary estimators from freeman2023linear,freeman2022multidimensional, we show that under low level regularity conditions on the function of unobserved heterogeneity, statistical inference on estimates follows. In particular, we show properties of eigenfunctions that offer an incidental debias in the tail terms in the functional singular value decomposition of functions in a Hilbert space, when the weighted-within transformation from freeman2022multidimensional is implemented in this setting. The second orthogonalisation over covariates in the Neyman orthogonal estimator allows this incidental debias to be asymptotically weaker than just applying the weighted-within transformation to the equation for dependent variables. The Neyman orthogonal estimator in general allows for weaker asymptotic convergence for estimates of fixed-effects in either the dependent or independent variables.