The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
74,561 characters
Weak instrumental variables due to nonlinearities in panel data: A Super Learner Control Function estimator
\maketitle
\begin{abstract}
A triangular structural panel data model with additive separable
individual-specific effects is used to model the causal effect of a
covariate on an outcome variable when there are unobservable confounders
with some of them time-invariant. In this setup, a linear reduced-form
equation might be problematic when the conditional mean of the endogenous
covariate and the instrumental variables is nonlinear. The reason is that
ignoring the nonlinearity could lead to weak instruments (instruments are
weakly correlated with the endogenous covariate).
As a solution, we propose a triangular simultaneous equation model
for panel data with additive separable individual-specific fixed effects
composed of a linear structural equation with a nonlinear reduced form equation.
The parameter of interest is the structural parameter of the endogenous variable.
The identification of this parameter is obtained under the assumption of
available exclusion restrictions and using a control function approach.
Estimating the parameter of interest is done using an estimator that
we call Super Learner Control Function estimator (SLCFE).
The estimation procedure is composed of two main steps
and cross-fitting. First, we estimate the control
function using a super learner. In the following step,
we use the estimated control function to control for
endogeneity in the structural equation. Cross-fitting is done across the individual dimension.
The estimator is consistent and asymptotically normal achieving a parametric rate of convergence.
We perform a Monte Carlo simulation to test
the performance of the estimators proposed.
We conclude that the Super Learner Control Function
Estimators significantly outperform 2SLS estimators. Finally, we show that
the SLCF estimator differs from both the plug-in IV estimator and a
naive plug-in 2SLS estimator, with the former being incosistent
without cross-fitting, and
the latter being inconsitent even with cross-fitting.
\end{abstract}
\keywords{Endogeneity,
Super Learner, Instrumental Variables,
Identification,
Control Function, Panel data, Individual Fixed Effects}
\section{Introduction}
A triangular structural panel data model with additive separable
individual-specific effects is used to model the causal effect of
a covariate on an outcome variable when there are unobservable
confounders, and some of them are time-invariant. In this setup,
it is common to assume that the structural is linear and to use
a linear projection in the reduced form equation. However,
a linear reduced-form equation might be problematic when the
conditional mean of the endogenous covariate and the instrumental
variables is nonlinear. The reason is that ignoring the nonlinearity
could lead to problems due to weak instruments since the strength of the
linear relationship might not be high. This paper proposes to exploit the
nonlinearties in the relationship between the endogenous covariate and
the instrumental variables as a way to avoid problems of weak instruments
due to nonlinearities. \footnote{Instruments are weak when they are
weakly correlated to the endogenous covariate.
In a cross-sectional setup, weak IVs causes that
the median of the 2SLS estimator is biased towards the OLS estimator,
2SLS t-test presents size inflation \citep{KeaneARE2024}. In panel data,
Within-2SLS estimator is consistent when the time dimension grows
\citep{CaiER2012}.}
An example is the relationship between air pollution and children's
educational outcomes, where both variables are influenced by
unobserved confounders that can be either time-varying or time-invariant.
For instance, neighborhood characteristics and family income are often
unobserved in publicly available data. While neighborhood characteristics
typically remain constant over time, family income fluctuates.
Both factors are relevant to a child's academic achievement
and exposure to air pollution, leading to endogeneity in pollution exposure.
To address this, fixed effects can be used to control for unobserved neighborhood
characteristics, while an instrumental variable can help account for time-varying
confounders. A commonly used instrumental variable is wind direction,
which has a well-established nonlinear relationship with air pollution levels
\citep{ZabrockiN2022}.
This paper studies a triangular structural panel data model in which the reduced
form equation is nonlinear and individual specific effects are present. It proposes
the use of Super Learning \footnote{It is feasible to use other machine learning
methods instead of Super Learning. However, Super Learning automatically combines
different methods in an optimal way achieving a convergence rate that is lower or
equal to the best learner included in the library of base learners.}
for the estimation of the nonlinearities in the model.
For simplicity consider the following triangular model without covariates for
$i \in \{1,...,N\}$, $t \in \{1,...,T\}$:
\begin{equation}\label{structural_equation_toy}
y_{it}=x_{1it}\beta_{1_o}+\alpha_{i,y}+\varepsilon_{it}, \quad with \quad
\mathbb{E}[\varepsilon_{it}| z_{i}, \alpha_{i,y}]=0,
\end{equation}
\begin{equation}\label{reducedform_equation_toy}
x_{1it}=g_o(z_{it})+\alpha_{i,1x}+u_{it}, \quad with \quad \mathbb{E}[u_{it}|z_{i},\alpha_{i,1x}]=0,
\end{equation}
\begin{equation}\label{relationship_errors_toy}
\varepsilon_{it}=\rho u_{it} + \omega_{it}, \quad with \quad
\mathbb{E}[\omega_{it} |u_{i}, z_i, \alpha_{i,y} ]=0,
\end{equation}
\noindent where $y_{it}$, and $x_{it}$ are continuous and bounded, $u_i' = [u_{i1}, u_{i2}, ..., u_{iT}]$, $z_i' = [z_{i1}, z_{i2}, ..., z_{iT}]$ are vectors stacking up all the observations available per individual. $\alpha_{i,y}$, and $\alpha_{i,1x}$ represent individual specific effects. \footnote{$\alpha_{i,y}$, and $\alpha_{i,1x}$ can be equal as in the cited example but for generality we allow them to be different.}
The model comprises a structural equation (\ref{structural_equation_toy}) with an endogenous regressor $x_{1it}$ that presents a non-linear relationship with the available instrumental variable $z_{it}$ (\ref{reducedform_equation_toy}). Both equations present individual effects that represent time-invariant unobserved confounders. In this setup, conditioning on the individual effects allows us to control for time-invariant unobserved confounders. But, estimating all the individual effects leads to incidental parameter bias as the number of parameters to estimate grows with the sample size. Then, in order to eliminate the individual effects, we propose to transform the model (e.g. first-differencing, or a within transformation which is equivalent to obtain the deviations of the original variables from their individual mean). On the other hand, the presence of time-varying unobserved confounders is modeled through a linear relationship between the error terms in both equations (\ref{relationship_errors_toy}). We deal with this issue by exploiting available exclusion restrictions (instrumental variables) along with a Control Function approach on the transformed model.
To estimate the causal parameter of interest $\beta_{1_o}$,
we provide a two-step estimator.
In the first step, we use a super learner \citep{VanSAMB2007} to
learn the nonlinear relationship of the transformed endogenous
variable with the instrumental variables and the observable confounders.
After this, we obtain the residuals of the transformed endogenous variable.
In the second step, we use the estimated residuals in the first stage as a
control function for the endogenous variable in the structural equation. We call this estimator Super Learner Control
Function (SLCF) estimator. The main results of the paper are that it the Super Learner Control Function
estimator is consistent and asymptotically normal with a convergence rate equal
to $\sqrt{N_T}$ \footnote{$N_T = \sum_i^N T_i$.}. An extensive Monte Carlo experiment is performed
to test the small sample properties of the proposed estimator.
We conclude that the estimator performs well provided that we can accurately
learn the nuisance parameter in the first stage.
An R package is available, if you are interested please contact the author.
Our paper contributes to the panel data literature in several ways. To the best of our knowledge,
it is the first one to study and establish the following results for panel data. First, this paper studies the use of machine learning to estimate
panel data models with an endogenous covariate (arising from unobserved confounders
or measurement error) and additive individual fixed effects using a control
function approach and exploiting nonlinearities. Second, it provides general theoretical results that apply to different
panel data transformations \footnote{Panel data transformations
are required to eliminate the individual specific effects} while using
machine learning methods and control function. Third, it shows that, under suitable
conditions, one can identify the transformed nuisance parameter and as a result
the transformed control function. Fourth, it shows that a control function estimator
with a parametric assumption in the structural equation is a double machine learning estimator.
Finally, it shows that the control function estimator is not equivalent to a
plug-in IV estimation.
In contrast, the available papers using machine learning to estimate panel data
models focus on exogenous covariates (time-varying confounders are observed) and
the presence of individual fixed effects and focus on specific panel data transformations.
For instance, \cite{klosinWP2022} proposes a
debiased LASSO estimator to estimate heterogeneous treatment effects of a
continuous treatment, which is modeled using a structural nonparametric panel
data model with additive fixed effects and strictly exogenous covariates. They
account for the individual specific effects by taking the first-difference of
the original model. \cite{ SemenovaQE2023} focuses on a structural equation
model for panel data allowing for heterogenous treatment effects with
unconfoundedness on observables (sequentially exogenous covariates) and individual specific fixed effects.
In order to deal with the individual fixed effects, the authors use a Mundlak approach. They provide an algorithm composed of three steps.
In the first step, using cross-fitting they estimate the nuisance functions of the structural and treatment equations to obtain the residuals. In the second step, they estimate the conditional average treatment effect using lasso estimators. In the third step, they provide debiased inference.
This paper also contributes to the literature on weak instruments in panel
data by being the first to focus on a setup with \textit{short time dimension} in
the presence of nonlinearities. While available papers on
weak instrumental variables in panel data focus on the asymptotic properties of
the Within 2SLS estimator showing that it is consistent even with weak instruments
when the \textit{time dimension grows to infinity}. \cite{CaiER2012} show that the Within 2SLS is consistent when the
individual and time dimension grow to infinity. They also present asymptotical
pivotal tests (Anderson-Rubin test, Kleibergen test) for longitudinal data.
In Section \ref{S_Examples}, we present motivating examples. In Section
\ref{S_model_line_nonlin}, we present the baseline model. In Section
\ref{S_Identification}, we describe the identification strategy. In Section
\ref{S_SLCF}, we present the SLCF estimator. In Section \ref{S_Properties}, we
present the large sample properties of the SLCF estimator.
In Section \ref{S_Comparison}, we present a comparison of the proposed estimator
with IV and 2SLS plugin estimators. In Section \ref{S_Simulation}, we present
the simulation experiment and its results,
in Section \ref{EA} we present the results of
the empirical application, and in
Section \ref{S_Conclusions} we present the
conclusions of the paper.
\section{Other Related literature}
This paper is related to the literature on semiparametric panel data models. \citet{RodriguezPoo2017JES} explains that one can either use a first-differencing approach or a profiling technique to estimate a semi-parametric model with an additive separable individual specific effect and a disturbance term. \citet{Li1996EL} explains that kernel estimation of the conditional expectation of the first-differenced endogenous regressor is not feasible when the dimension of the conditioning set is larger than 5. As a solution, \citet{BaltagiASF2002} presents an additive estimation of the first-differenced non-linear function.
In addition, this paper is related to the literature on orthogonal/debiased estimation, which focuses on providing $\sqrt{N}$ consistent estimators of low-dimensional parameters in the presence of high-dimensional nuisance parameters. \citet{Chernozhukov2018EJ} defines an orthogonal score as the one that presents a vanishing Gateaux derivative with respect to the nuisance parameters when evaluated at the true finite-dimensional parameter values. The orthogonal score is closely related to the score proposed by \citet{Robinson1988E}.
Finally, our work relates to the literature of mixed effects models and machine learning. Our model is related to a mixed effects model assuming that covariates are exogenous conditioning on the cluster random effects. \cite{EmmeneggerSJS2023} uses double machine learning to estimate a partially linear mixed-effects model. \cite{Young2024JRSSB} proposes a sandwich boosting for accurately estimating partially linear models for grouped data.
\subsection{Super Learner}
The Super Learner was proposed by \cite{VanSAMB2007}, it is a heterogeneous ensemble method within
Supervised Learning. It consists of different base learners that are trained on the training data set, then
their predictions are combined using a meta-learner by minimizing the cross-validated risk. The meta-learner
can be a linear regression of the predictions of each learner on the outcome variable. \cite{VanSAMB2007} shows that
the Super Learner performs asymptotically as well as or better than learner in the library of base learners.
\section{Motivating examples}\label{S_Examples}
As mentioned in the Introduction, measuring the causal effect of air pollution on
health and educational outcomes is not straightforward as it is not randomly
assigned. More specifically, time-varying unobservables such as family income might
affect air pollution and health or educational outcomes. In addition, unobserved
time-invariant variables, such as neighborhood characteristics or genetics,
might also be present. While time-invariant unobservables could be dealt with
by including individual fixed effects in the model, the presence of time-varying
unobservable confounders requires instrumental variables. In this setup, wind
direction is a commonly used as an instrumental variable. While it is well
documented that the relationship between air pollution and wind direction is
nonlinear \citep{ZabrockiN2022}, it is common to see that the nonlinearity is
ignored or modeled with a parametric assumption.
\subsection{Air pollution and educational outcomes}
\cite{GilraineNBER2022} studies the effect of air pollution on educational outcomes. They specify the following linear regression for subject $s$, cohort $c$, district $d$, year $t$:
\begin{equation}
y_{sdct}=\beta PM2.5_{dt}+X_{sdct}'\gamma+W_{dt}'\eta+\omega_s+\alpha_{d} + \theta_{c}+\lambda_t+\epsilon_{scdt},
\end{equation}
\noindent where $y_{sdct}$ represents the average grade of subject $s$, cohort $c$,
district $d$ at year $t$. $PM2.5_{dt}$ represents the average level of $PM2.5$
at district $d$ in year $t$, $X_{sdct}$ is a vector of control variables, $W_{dt}$
is a vector of weather variables, $\omega_s$ represents subject specific effects,
$\alpha_{d}$ represents district fixed effects, $\theta_{c}$ represents cohort
fixed effects, $\lambda_t$ are time fixed effects, and $\epsilon_{scdt}$ is the disturbance term.
As mentioned before, $PM2.5_{dt}$ is not randomly assigned and subject to measurement error. Thus, the authors use year-to-year changes in coal production and a shift-share instrument equal to the interaction of fuel shares used for nearby production with national growth rates along with a linear first-stage regression.
As an alternative, following \cite{DeryuginaAER2019} we could use wind direction as an instrument of air pollution, but with a nonlinear
first-stage equation as follows:
\begin{equation} PM25_{dt}=h(WD_{dt},X_{dt})+\gamma_{d}+\varepsilon_{dt},
\end{equation}
\noindent where $WD_{dt}$ is the average wind direction in district $d$ in year $t$, and $\gamma_{d}$ represent district specific effects.
\subsection{Air pollution and health outcomes}
\cite{DeryuginaNBERWP2023} studies the effect of air pollution on mortality rates. More precisely, they specify the following linear regression for county $c$, day $d$, month $m$ and year $y$:
\begin{equation}
y_{cdmy}^k=\beta SO2_{cdmy}+X_{cdmy}'\gamma+ \alpha_{cm}+ \alpha_{my}+\epsilon_{cdmy},
\end{equation}
\noindent where $y_{cdmy}^k$ represents the cumulative mortality rate $k$ days after
day $d$, $SO2_{cdmy}$ represents the level of $SO2$ at day $d$ in county $c$,
$X_{cdmy}$ is a vector of control variables that include precipitation,
wind speed, and temperature, $\alpha_{cm}$ represents county-month fixed effects,
$\alpha_{my}$ month-by-year fixed effects,
and $\epsilon_{cdmy}$ is the disturbance term.
For identification, the authors use wind direction changes as an instrument for air pollution. The authors assume that the first-stage regression presents a parametric nonlinear specification.
\begin{equation}
SO2_{cdmy}=\sum_{g=1}^{50} f^g(\theta_{cdmy})+X_{cdmy}'\delta+ \alpha_{cm} + \alpha_{my}+\varepsilon_{cd},
\end{equation}
\noindent where:
\begin{equation}
f^g(\theta_{cdmy})=\gamma_g^1 \mathbf{1}\{G_c=g\} sin(\theta_{cdmy})+\gamma_g^2\mathbf{1}\{G_c=g\} sin(\theta_{cdmy}/2).
\end{equation}
In this equation, $\mathbf{1}\{G_c=g\}$ is an indicator function that is equal to 1
if county c is member of group $g$ and 0 otherwise. $\theta_{cdmy}$ is local wind
direction measured in radians. There are 100 instrumental variables.
The authors show that their results are robust to different parametric specifications of the first-stage regression but do not allow for a nonparametric specification. However, the concern is that the parametric specifications do not fully capture the nonlinearity in the first-stage equation. Thus, we propose to use a super learner to estimate the nonlinear first-stage regression as follows:
\begin{equation}
SO2_{cdmy}=h(\theta_{cdmy},X_{cdmy}') + \alpha_{cm}+ \alpha_{my}+\varepsilon_{cd},
\end{equation}
\noindent where $h(\cdot)$ is an unknown nonlinear function.
\section{ The Setup}\label{S_model_line_nonlin}
We study a model composed of a linear structural equation, the outcome variable ($y_{it}$) explained by an endogenous variable ($x_{1it}$), $K-1$ exogenous regressors ($\tilde{x}_{it}$), and an additive individual-specific effect ($\alpha_{i,y}$). The model presents an external instrumental variable ($z_{it}$) for the endogenous regressor. The endogenous regressor is mean dependent on the exogenous regressors, the instrumental variable, and the individual specific effects $\alpha_{i,x}$. We assume that:
\sloppy
\begin{assumption}\label{A_data}
We observe a sequence of data sets $\{O_i = (y_{i},x_{1i},\ldots,x_{Ki},z_{i})\}_{i=1}^N$,
which are $T_i\times 1$ independent copies of the $K+2$ random vectors
$(\mathbf{y},\mathbf{x}_1,\mathbf{x}_2, \ldots,\mathbf{x}_K,\mathbf{z})$
taking values in $[-M_y,M_y]^T \times [-M_{x1},M_{x1}]^T \times [-M_{xK},M_{xK}]^T \times [-M_z,M_z]^T$
with $M_y<\infty$, $M_{xk}<\infty$, $M_z<\infty$. Thus, each element $t$ of
$\mathbf{y}$ is bounded and contained in $[-M_y,M_y]$. Similarly, each element of
$x_k$ is contained in $[-M_{xk},M_{xk}]$, which means that $\mathbf{x}_{kt}$
is a continuously distributed random variable with compact support $[-M_{xk},M_{xk}]$.
In addition, each element of $\mathbf{z}$ is a continuously distributed random variable with compact support $[-M_z,M_z]$.
Finally, $T_i$ is uniformly bounded such that $T_i < T_{max}$.
\end{assumption}
Assumption \ref{A_data} states that each individual $i$ presents $T_i$ observations.
The total number of observations is denoted by $N_T = \sum_i T_i$.
\begin{assumption}\label{A_structural_eq}
The structural equation has a linear form with an additive unobserved error term and an additive individual specific effect $\alpha_{i,y}$
\begin{equation}
\begin{gathered}
y_{it}=x_{1it}\beta_{1_o}+\tilde{x}_{it}'\beta_{2_o}+\alpha_{i,y}+\varepsilon_{it}, \quad i \in \{1,...,N\}, t \in \{1,...,T_i\}, \\
\mathbb{E}[\varepsilon_{it}|\tilde{x}_{i}, z_{i}, \alpha_{i,y}]=\mathbb{E}[\varepsilon_{it}]=0.
\end{gathered}
\end{equation}
\noindent where $\tilde{x}_{it}=[x_{2it}, x_{3it}, ..., x_{Kit}]'$.
\noindent The unobserved random term $\varepsilon_{it}$ is independent and identically
distributed with a continuous cummulative distribution function. It has a compact
support with zero mean conditional on $\tilde{x}_{i}=[\tilde{x}_{i1}, \tilde{x}_{i2}, \ldots,
\tilde{x}_{iT_i}]'$, $z_{i}=[z_{i1}, z_{i2},..., z_{iT_i}]'$, $\alpha_{i,y}$, and
$Var(\varepsilon_{it}|\tilde{x}_{i}, z_{i}, \alpha_{i,y})< \infty$. $\alpha_{i,y}$
are individual specific
effects with sigma-field $\mathcal{A}_y$.
The unobserved parameter
of interest is $\beta_{1_o}$.
\end{assumption}
\noindent Under assumption \ref{A_structural_eq}, $\tilde{x}_{it}$ are strictly exogenous in
the structural equation conditional on $\alpha_{i,y}$. In addition, $\{z_{is}, \quad s=1,\dots,T \}$
are uncorrelated with the structural error term $\varepsilon_{it}$ conditional on $\alpha_{i,y}$.
\begin{assumption}\label{A_reduced_form}
The regressor $x_{1it}$ is endogenous and has an unknown nonlinear reduced form equation with an
additive separable individual specific effect $\alpha_{i,x}$:
\begin{equation}
\begin{gathered}
x_{1it}=g_o(\tilde{x}_{it},z_{it})+\alpha_{i,1x}+u_{it}, \quad i \in \{1,...,N\}, t \in \{1,...,T_i\}, \\
\mathbb{E}[u_{it}|\tilde{x}_{i},z_{i},\alpha_{i,1x}]=0,\\
\end{gathered}
\end{equation}
\noindent with the unknown function $g_o: [-M_x,M_x]^K \times [-M_z,M_z]\rightarrow \mathcal{S} \subset \mathbb{R}$, with
$g_o \in \mathcal{H}(\varsigma,C)$ where $\mathcal{H}(\varsigma, C)$ represents a Holder-class with smoothness parameter
equal to $\varsigma \geq 2$\footnote{If the unknown function is linear, the SLCF estimator is still consistent
if instruments are not weak. We argue that if the unknown function $g_o$ is linear,
the library of the base learners includes the linear regression, and the
instruments are appropriately transformed, then the SLCF estimator is equivalent to the 2SLS estimator.
A proof of this statement is out of the scope of this paper.}. $\alpha_{i,x}$ are individual specific
effects measureble with respect $\mathcal{A}_y$, therefore knowledge of $\alpha_{i,y}$ implies
knowledge of $\alpha_{i,x}$. The reduced form error term $u_{it}$ is independent and identically distributed
with a continuous cummulative distribution function. It has a compact support
with zero mean conditional on $\tilde{x}_{i}=[\tilde{x}_{i1},
\tilde{x}_{i2}, \ldots, \tilde{x}_{iT_i}]'$,
$z_{i}=[z_{i1}, z_{i2},..., z_{iT_i}]'$, and $Var(u_{it}|\tilde{x}_{i}, z_{i}, \alpha_{i,1x})< \infty$.
\end{assumption}
\noindent Under assumption \ref{A_reduced_form}, $\tilde{x}_{it}$, $z_{it}$ are strictly exogenous in the reduced form equation conditional on the individual specific effect $\alpha_{i,1x}$. In addition, this assumption states that the endogenous regressor has a nonlinear relationship with the exogenous regressors $\tilde{x}_{it}$, and the variable $z_{it}$.
As a consequence of Assumptions \ref{A_structural_eq} and \ref{A_reduced_form}, $z_{it}$ is an instrumental variable.
\noindent The decomposition of $x_{1it}$ into a conditional expectation and an additive disturbance term is appropriate under the assumption that $x_{1it}$ is continuous (Assumption \ref{A_data}). Assumption \ref{A_reduced_form} needs to be modified when $x_{1it}$ is a non-continuous endogenous explanatory variable.
\begin{assumption} \label{A_eps_u_omega}
The structural error term $\varepsilon_{it}$ has a linear relationship with $u_{it}$:
\begin{equation}
\begin{gathered}
\varepsilon_{it}=\rho u_{it} + \omega_{it} , \quad i \in \{1,...,N\}, t \in \{1,...,T_i\}, \\
\mathbb{E}[ \omega_{it} |u_{i}]=0.
\end{gathered}
\end{equation}
\noindent The error term $\omega_{it}$ is independent and identically
distributed with a continuous cummulative distribution. It has a bounded compact
support with zero mean conditional on $u_{i}=[u_{i1}, u_{i2}, \dots, u_{iT_i}]'$,
$\tilde{x}_{i}=[\tilde{x}_{i1}, \tilde{x}_{i2}, \ldots, \tilde{x}_{iT_i}]'$,
$z_{i}=[z_{i1}, z_{i2},\ldots, z_{iT_i}]'$, $\alpha_{i,y}$,
and $Var(\omega_{it}|\tilde{x}_{i}, z_{i}, \alpha_{i,y}) = \sigma^2_{\omega}< \infty$.
\end{assumption}
As a consequence of Assumption \ref{A_eps_u_omega}, $\omega_{it}$ is uncorrelated with $u_{it}$.
\begin{remark}
The assumption of a linear relationship between $\varepsilon_{it}$ and $u_{it}$ is a
crucial assumption in panel data setups.
Relaxing this assumption to allow for a nonlinear relationship with $u_{it}$ is not trivial even
when there are only additive individual fixed effects. Since one needs to eliminate fixed effects
through transforming the data, one can only identify the transformed reduced form erros
under assumptions \ref{A_structural_eq}, and \ref{A_reduced_form}.
However, identification under a nonlinear function of $u_{it}$ requires knowledge of the
untransformed reduced form errors which is not feasible (See Section \ref{S_Identification}).
For example, in the case of first-differencing, we only identify the
first-differenced reduced-form errors but we would need the current and first-lag of
the reduced form errors in the required information set. These reduced-form errors
are not separetely identifiable under assumptions \ref{A_structural_eq}, and \ref{A_reduced_form}.
\end{remark}
As consequence of Assumptions \ref{A_reduced_form} and \ref{A_eps_u_omega},
the structural error term $\varepsilon_{it}$ is correlated with $x_{1it}$. The reason
is that $\mathbb{E}[x_{1it}\varepsilon_{it}]=\mathbb{E}[(g_o(\tilde{x}_{it},z_{it})+\alpha_{i,x}+u_{it})\varepsilon_{it}]$ by Assumption \ref{A_reduced_form}. The last expression is equal to $\mathbb{E}[(g_o(\tilde{x}_{it},z_{it})+\alpha_{i,x}+u_{it})(\rho_o u_{it}+\omega_{it})]=\rho_o\mathbb{E}[u_{it}^2]$
by Assumptions \ref{A_reduced_form} and \ref{A_eps_u_omega}.
\noindent The error term $\omega_{it}$ has zero mean conditional on
$x_{1it}$, $\tilde{x}_{it}$, $\alpha_{i,y}$, and $u_{it}$
($\mathbb{E}[\omega_{it}|x_{1it},\tilde{x}_{it},u_{it},\alpha_{i,y}]=0$)
because $\mathbb{E}[ \omega_{it} |x_{1i}, \tilde{x}_{i}, z_i, \alpha_{i,y}]
=\mathbb{E}[\omega_{it}|u_{1i}, \tilde{x}_{i},z_i,\alpha_{i,y}]=\mathbb{E}[\omega_{it}\mid u_i]=0$. The first equality
is a consequence of the fact that $x_{1it}$ is a one-to-one function
with $u_{it}$ conditional on $\tilde{x}_{it}$, $z_{it}$, $\alpha_{i,1x}
\subseteq \alpha_{i,y}$ \citep{Wooldridge2010Book}. The second equality is a consequence of the mean independence of
$\varepsilon_{it}$ on the exogeneous covariates, the instrumental variables,
and the individual specific effects $\alpha_{i,y}$ (Assumption \ref{A_structural_eq}).
\section{Identification: Control Function Approach} \label{S_Identification}
The identification of the parameter of interest $\beta_{1_o}$ relies on the presence of an instrumental variable $z_{it}$ conditioning on the individual-specific effects (Assumptions \ref{A_structural_eq} and \ref{A_reduced_form}). Since we run into an incidental parameter problem, we get rid of the individual-specific effects by using an appropriate transformation of the data $\tau: \mathbb{R} \rightarrow \mathbb{R}$. The typical transformations of the data are first-differencing and within transformation. The first-difference transformation is equal to $y_{it}-y_{it-1}$, and the Within transformation is equal to $y_{it}-\sum_{t=1}^{T_i} y_{it}/T_i$.
\begin{equation} \label{dif_structural}
\tau y_{it}=\tau x_{1it}\beta_{1_o}+\sum_{k=2}^K \tau x_{kit}\beta_{2k_o}+\tau\varepsilon_{it}, i\in \{1, ..., N\}, t \in \{t_a, ..., T_i\},
\end{equation}
\begin{equation} \label{dif_reduced}
\tau x_{1it}=\tau g_o(\tilde{x}_{it}, z_{it})+\tau u_{it}, i\in \{1, ..., N\}, t \in \{t_a, ..., T_i\}.
\end{equation}
\noindent where $t_a$ is the first time period available after transforming the data. In the case of first-differencing $t_a=2$, and for a Within-transformation $t_a=1$.
Now, it is tempting to replace equation \ref{dif_reduced} into \ref{dif_structural} to obtain:
\begin{equation}\label{naive_model}
\tau y_{it}=\tau g_o(\tilde{x}_{it},z_{it})\beta_{1_o}+\sum_{k=2}^K \tau x_{kit}\beta_{2k_o}+\tau u_{it}\beta_{1_o}+\tau\varepsilon_{it}, i\in \{1, ..., N\}, t \in \{t_a, ..., T_i\}.
\end{equation}
\noindent In equation \ref{naive_model}, one can believe that the estimation of the
unknown parameters is straightforward since the error term
$\tau u_{it}\beta_{1_o}+\tau \varepsilon_{it}$ in \ref{naive_model} has zero mean
conditional, on an appropriate information set, as a consequence of assumptions
\ref{A_structural_eq}, and \ref{A_reduced_form}. Moreover, one can think that
the best estimation procedure is to estimate $\tau g_o(\tilde{x}_{it},z_{it})$
in the first stage and in the second stage, plug it in equation
\ref{naive_model} and perform simple OLS. But this is not possible because the
regularization bias of the first stage estimator of
$\tau g_o(\tilde{x}_{it},z_{it})$ contaminates the estimation of $\beta_{1_o}$
and $\beta_{2_o}$ \citep{Chernozhukov2018EJ, Robinson1988E, GuoWP2022} and
produces estimates that are not $\sqrt{N_T}$-consistent (See Section \ref{S_Comparison} for more details).
In order to deal with this problem, we can use a control function approach and augment the structural equation by controlling for the unobservable $\tau u_{it}$ as follows:
\begin{equation} \label{ControlFunction_eq}
\tau y_{it}=\tau x_{1it}\beta_{1_o}+\sum_{k=2}^K \tau x_{kit}\beta_{2k_o}+\rho_o \tau u_{it} + \tau \omega_{it}, i\in \{1, ..., N\}, t \in \{t_a, ..., T_i\}.
\end{equation}
In model \ref{ControlFunction_eq}, under assumptions \ref{A_structural_eq} - \ref{A_eps_u_omega} we can guarantee zero correlation between
$\tau x_{1it}$, and $\tau \omega_{it}$ thanks to the presence of an instrumental variable in the reduced form equation that is mean independent
of the error term $\omega_{it}$.
Stacking up the time observations and letting $M_{\tau_i}$ denote
the transformation operator mapping
$\mathbb{R}^{T_i} \rightarrow \mathbb{R}^{T_{i,a}}$, with $T_{i,a}$ equal
to $T_i$ under the within transformation and $T_i-1$ under the
first-difference transformation, we obtain the following moment conditions:
\begin{equation}\label{non_orthogonal_score}
\mathbb{E}[\phi(O_i;\theta_o,\tau g_o)]=\textbf{0},
\end{equation}
\noindent where
$\phi(O_i;\theta_o,\tau g_o)$ is a vector of score functions equal to $(M_{\tau_i} H_{i})' V_i^{-1} M_{\tau_i} \omega_{i} $, $\theta_o=[\beta_{1_o} \quad \beta_{2_o}' \quad \rho_o]'$, $M_{\tau_i} \omega_{i}=M_{\tau_i} y_{i}-M_{\tau_i} x_{1i}\beta_{1_o}-\sum_{k=2}^K M_{\tau_i} x_{ki}\beta_{2k_o}-\rho_o M_{\tau_i} u_{i}$, $\textbf{0}=[0\quad \cdots \quad 0\quad 0]'$, and $M_{\tau_i} H_{i}=[M_{\tau_i} x_{1i}\quad M_{\tau_i}\tilde{x}_{i}\quad M_{\tau_i} u_{i}]$, and $V_i = \mathbb{E}[(M_{\tau_i}\omega_i)(M_{\tau_i}\omega_i)'| H_i, z_i]=\sigma^2_{\omega}M_{\tau_i}M_{\tau_i}'=\sigma^2_{\omega} \tilde{V}_i $.
\begin{remark}
If we use the first-difference transformation, $M_{\tau_i}=D_i=\begin{bmatrix}
1 & -1 & 0& \cdots & 0 \\
0 &1 & -1 & \cdots & 0 \\
\vdots & & \ddots & \ddots & \vdots \\
0 & \cdots & 0 & 1 & -1 \\
\end{bmatrix}$, with dimensions $(T_i-1) \times T_i $. If we use the within transformation, $M_{\tau_i}= W_i= I_i-1/T_i\begin{bmatrix}
1 & 1 & 1& \cdots & 1 \\
1 &1 & 1 & \cdots & 1 \\
\vdots & \vdots & \ddots & \ddots & \vdots \\
1 & \cdots & 1 & 1 & 1 \\
\end{bmatrix}$ , with dimensions $T_i \times T_i $, $I_i$ is the identity matrix.
\end{remark}
\begin{remark}
The moment condition is valid even under misspecification of $V_i$ because we assume strict exogeneity. If we would like to relax the assumption of strict exogeneity to weak exogeneity, then it is better to use $\mathbb{E}[(M_{\tau_i}H_i)'M_{\tau_i}\omega_i]=0$ as moment conditions.
\end{remark}
\begin{remark}
If $M_{\tau_i}=W_i$ (the within transformation), the moment condition can be
simply written as $\mathbb{E}[(M_{\tau_i}H_i)'M_{\tau_i}\omega_i]=0$.
\end{remark}
Finally, we identify $\tau u_{it}$ since it has zero mean conditional on an
appropriate information set $I_{t}=\{\tilde x_{it-j}, z_{it-j}, L_i \leq j \leq U_i \}$
with $L_i \geq t_i - T_i$ and $U_i \leq t_i-1$ as a consequence of Assumption
\ref{A_reduced_form}. This identification strategy was first proposed
by \cite{AvilaMarquez2023}, and it was presented at the 28th International Panel Data Conference at the University of Amsterdam in 2023.
\begin{example_un}{ First-difference transformation:}
We could follow \cite{BaltagiASF2002} by first-differencing the structural
and reduced-form equations such that we obtain:
\begin{equation} \label{dif_structural_ex}
\Delta y_{it}=\Delta x_{1it}\beta_{1_o}+\Delta\tilde{x}_{it}'\beta_{2_o}+\Delta\varepsilon_{it}, i\in \{1, ..., N\}, t \in \{2, ..., T_i\}.
\end{equation}
\begin{equation} \label{dif_reduced_ex}
\Delta x_{1it}=\Delta g_o(\tilde{x}_{it}, z_{it})+\Delta u_{it}, i\in \{1, ..., N\}, t \in \{2, ..., T_i\}.
\end{equation}
Now, we can use the proposed control function approach and augment the structural equation by controlling for the unobservable $\Delta u_{it}$ as follows:
\begin{equation} \label{ControlFunction_eq_ex}
\Delta y_{it}=\Delta x_{1it}\beta_{1_o}+\Delta \tilde{x}_{it}'\beta_{2_o}+\rho_o \Delta u_{it} + \Delta \omega_{it}, i\in \{1, ..., N\}, t \in \{2, ..., T_i\}.
\end{equation}
Next, we set up the following moment conditions:
\begin{equation}\label{non_orthogonal_score_ex}
\mathbb{E}[(D_i H_{i})' V_i^{-1} D_i \omega_{i}]=\textbf{0},
\end{equation}
\noindent where
$D_i \omega_{i}=D_i y_{i}-D_i x_{1i}\beta_{1_o}-D_i \tilde{x}_{i}'\beta_{2_o}-\rho_o D_i u_{i}$, $\textbf{0}=[0\quad \cdots \quad 0\quad 0]'$, $D_i H_{i}=[D_i x_{1i}\quad D_i\tilde{x}_{i}'\quad D_i u_{i}]'$, and $V_i = \sigma^2_{\omega} D_iD_i'$.
Finally, we identify $\Delta u_{it}$ since it has zero mean conditional on $\tilde{x}_{it}$, $\tilde{x}_{it-1}$, $z_{it}$, $z_{it-1}$ as a consequence of Assumption \ref{A_reduced_form}.
\begin{proposition}
Under assumption \ref{A_reduced_form}, $\Delta g({\tilde{x}_{it},z_{it}})$ is identified by $\mathbb{E}[\Delta x_{1it} | \tilde{x}_{it},z_{it},\tilde{x}_{it-1},z_{it-1}]$.
\end{proposition}
\begin{proof}
\begin{equation}
\mathbb{E}[\Delta x_{1it}|\tilde{x}_{it},z_{it},\tilde{x}_{it-1},z_{it-1}] = \mathbb{E}[\Delta g(\tilde{x}_{it},z_{it})|\tilde{x}_{it},z_{it},\tilde{x}_{it-1},z_{it-1}] + \mathbb{E}[\Delta u_{it}|\tilde{x}_{it},z_{it},\tilde{x}_{it-1},z_{it-1}]
\end{equation}
By the linearity property of the expectation operator, we have that the first term is equal to:
\begin{equation}
\begin{split}
\mathbb{E}[ g(\tilde{x}_{it},z_{it})|\tilde{x}_{it},z_{it},\tilde{x}_{it-1},z_{it-1}] - \mathbb{E}[ g(\tilde{x}_{it-1},z_{it-1})|\tilde{x}_{it},z_{it},\tilde{x}_{it-1},z_{it-1}] = \\
g(\tilde{x}_{it},z_{it})- g(\tilde{x}_{it-1},z_{it-1})
\end{split}
\end{equation}
Similarly for the second term, by the Law of Iterated Expectations and Assumption \ref{A_reduced_form}:
\begin{equation}
\mathbb{E}[ u_{it}|\tilde{x}_{it},z_{it},\tilde{x}_{it-1},z_{it-1}] - \mathbb{E}[ u_{it-1}|\tilde{x}_{it},z_{it},\tilde{x}_{it-1},z_{it-1}] =0
\end{equation}
\end{proof}
\end{example_un}
\begin{remark}
A similar analysis and conclusion is obtained using a Within transformation
which is equivalent to time demeaning the data. In this case, the information
set $I_{t} = \{ \tilde{x}_{it-j}, z_{it-j}, -T_i+t \leq j \leq t-1\}$.
\end{remark}
In Section \ref{S_SLCF}, we propose an estimation method based on the Control Function approach.
\section{Super Learner Control Function estimation}\label{S_SLCF}
As described in the previous section, we can write the transformed structural equation for observation $i$ at period $t$ as follows:
\begin{equation}
\tau y_{it}=\tau x_{1it}\beta_{1_o}+\tau \tilde{x}_{it}'\beta_{2_o}+\rho_o \tau u_{it} + \tau \omega_{it} , i\in \{1, ..., N\}, t \in \{t_a, ..., T_i\}.
\end{equation}
In this transformed equation, the parameters of interest are identified because $\tau x_{1it}$, $\tau \tilde{x}_{it}$ and $\tau u_{it}$ are uncorrelated with error term $\tau\omega_{it}$ under assumptions \ref{A_structural_eq} - \ref{A_eps_u_omega}.
In addition, it is clear that if $\tau u_{it}$ is observed, we could estimate the three parameters $\beta_{1_o}$, $\beta_{2_o}$, $\rho_o$ using the sample counterparts of the population moment conditions \ref{non_orthogonal_score}. But since $\tau u_{it}$ is not observed, we propose the following estimation procedure:
\noindent \textbf{First Step}:
\noindent We partition the set $\{1, 2, ..., N\}$ in B subsets $S_1$, $S_2$, ..., $S_B$, and denote $n_{T,b} = \sum_{i \in S_b} T_i$.
This sample splitting is equal to the one proposed by \cite{EmmeneggerSJS2023} and it was developed without knowledge of it. Alternatively, if the cluster structure in the data is known (individuals exhibit correlation within clusters but not across them, and both cluster membership and the number of clusters are identified) the data should be partitioned accordingly.
\noindent Then, for each $b \in{B}$ we estimate $\tau g_o(\tilde{x}_{it},z_{it}) = \mathbb{E}[\tau x_{1it}|I_t]$ using a super learner with partition $S_b^c= \{i \in S_g , g \neq b\}$, and we call the estimation $\widehat{\tau g_o}^{S_b^c}$. Then, we obtain the residuals $\widehat{\tau u}_{it}^{S_b}= \tau x_{1it}-\widehat{\tau g_o^{S_b^c}(\tilde{x}_{it},z_{it})}$ for partition $S_b =\{i \in S_b\}$.
Finally, the estimator of $\theta_o$ for partition $S_b$ is the solution of
the sample moment conditions such that:
\begin{equation}\label{sample_moments}
\frac{1}{n_{T,b}}\sum_{i \in S_b}\phi(O_i;\hat{\theta}_b,\widehat{\tau g_o}^{S_b^c})=0,
\end{equation}
\noindent where $\phi(O_i;\hat{\theta}_b,\widehat{\tau g_o}^{S_b^c}) = (\tilde{H}_{i}^{S_b})'V_i^{-1}M_{\tau_i}(y_i - \tilde{H_i}^{S_b}\hat{\theta}_b)$, $\tilde{ H}_{i}^{S_b}=[M_{\tau_i} x_{1i} \quad M_{\tau_i} \tilde{X}_{i} \quad \widehat{\tau u}_{i}^{S_b}]$, $\widehat{\tau u}_{i}^{S_b}$ is the vector stacking the residuals obtained in the first stage for partition $S_b$, $\tilde{X}_i$ a $T_i\times (K-1)$ vector collecting all time observations of $\tilde{X}_{it}$.
This is equivalent to estimate $\theta_o$ by performing generalized least squares
regression of $\tau y_{it}$ on $\tau x_{1it}$, $\tau \tilde{x}_{it}$ and
$\widehat{\tau u}_{it}^{S_b}$ or minimizing the following quadratic loss function:
$$\hat{\theta}_b=\substack{argmin\\ \theta} Q_b(O_i;\theta ,\widehat{\tau g_o}^{S_b^c}),$$
with $Q_b(O_i;\theta ,\widehat{\tau g_o}^{S_b^c})=\frac{1}{n_{T,b}}
\sum_{i \in S_b}(M_{\tau_i} y_{i}-M_{\tau_i} x_{1i}\beta_{1}-M_{\tau_i}
\tilde{X}_{i}\beta_{2}-\widehat{\tau u_{i}}^{S_b}\rho)'V_i^{-1}(M_{\tau_i}
y_{i}-M_{\tau_i} x_{1i}\beta_{1}-M_{\tau_i} \tilde{X}_{i}\beta_{2}-
\widehat{\tau u_{i}}^{S_b}\rho)$. Importantly, we can just perform ordinary
least squares estimation and use a sandwich type variance-covariance
matrix estimator for robust inference. In addition, OLS is the most efficient if $M_{\tau_i}=D_i$,
$T = 2$, and $\omega_{it}$ is homoskedastic and uncorrelated across both $i$
and $t$.
\noindent \textbf{Second Step}:
We average $\hat{\theta}_b$ to obtain the cross-fitting estimator of $\theta_o$ as:
$$\hat{\theta}_o=\frac{1}{B}\sum_{b=1}^B \hat{\theta}_b.$$
Thus, the cross-fitting estimator of our parameter of interest $\beta_{1o}$ is:
$$\hat{\beta}_{1o}=\frac{1}{B}\sum_{b=1}^B \hat{\beta}_b.$$
In order to deal with the dependency of the estimator $\hat{\beta}_{1o}$ on
the particular sample split, we
follow \cite{Chernozhukov2018EJ} by proposing to repeate the procedure $ss$
times. The estimates are aggregated either by using the mean or the median.
In addition, a correction term is added to the variance estimator of
$\hat{\beta}_{1o}$ (See Subsection \ref{Including_sample_splitting_uncertainty}).
\begin{algorithm}
\caption{Super Learner Control Function estimation}\label{Algortithm_CF1}
\begin{algorithmic}
\STATE \textbf{Input:} Panel data sets $\{y_i,X_{i},z_{i}\}_{i\in \{1, 2, ..., N\}}$ from model satisfying assumptions \ref{A_data} to \ref{A_eps_u_omega}, a natural number B.
\STATE \textbf{Output:} An estimator of the parameter of interest $\beta_{1o}$.
\FOR{$ss \in \{ 1, \dots, SS\} $}
\STATE Split the individual index set $\{1, 2, ..., N\}$ into B sets $S_1, S_2,...,S_B$.
\FOR{$b \in \{1, 2, ..., B\}$}
\STATE Estimate the conditional expectation $\mathbb{E}[\tau x_{1it}|I_t]$ using a super learner and data corresponding to $S_b^c$. \STATE Predict $\mathbb{E}[\tau x_{1it}|I_t]$ for $i \in S_b$.
\STATE Estimate the residuals $ \tau u_{it}$ for $i \in S_b$ using the prediction of $\mathbb{E}[\tau x_{1it}|I_t]$ of previous step.
\STATE Estimate $\theta_{o} = [\beta_{1o}, \beta_{2o}, \rho_o]'$ solving the sample moment conditions \ref{sample_moments} for $i \in S_b$.
\ENDFOR
\STATE Compute $\hat{\theta}_{ss}=\frac{\sum_1^B \widehat{\theta_{b,ss}}}{B}.$
\ENDFOR
\STATE Compute $\hat{\theta}_{o}=\frac{\sum_1^{SS} \widehat{\theta_{ss}}}{SS}.$
\end{algorithmic}
\end{algorithm}
\begin{example_un}{First-difference transformation (continued):}
As described in the previous section, we can write the first-differenced structural equation for observation $i$ at period $t$ as follows:
\begin{equation}
\Delta y_{it}=\Delta x_{1it}\beta_{1_o}+\Delta \tilde{x}_{it}'\beta_{2_o}+\rho_o \Delta u_{it} + \Delta \omega_{it} , i\in \{1, ..., N\}, t \in \{2, ..., T_i\}.
\end{equation}
\sloppy
First, we partition the set $\{1, 2, ..., N\}$ in B subsets $S_1$, $S_2$, ..., $S_B$,
and for each $b \in B$ we estimate $\mathbb{E}[\Delta x_{1it}|\tilde{x}_{it},\tilde{x}_{it-1},z_{it},z_{it-1}]$ using a super learner with partition $S_b^c$. Then, we get the residuals $\widehat{\Delta u}_{it}^{S_b}=\Delta x_{1it}-\widehat{\Delta g_o^{S_b^c}(\tilde{x}_{it},z_{it})}$ for partition $S_b$.
We obtain the estimator of $\theta_o$ for partition $S_b$ as the solution of the sample moment conditions such that:
\begin{equation}\label{sample_moments_ex}
\frac{1}{n_{T,b}}\sum_{i \in S_b}\phi(O_i;\hat{\theta}_b,\widehat{\Delta g}_o^{S_b^c})=0,
\end{equation}
\noindent where $\phi(O_i;\hat{\theta}_b,\widehat{\Delta g}_o^{S_b^c}) = (\tilde{H}_i^{S_b})' V_i^{-1} (D_i y_i - \tilde{H}_i^{S_b} \hat{\theta}_b) $, $\tilde{H}_{i}^{S_b}=[D_i x_{1i} \quad D_i \tilde{x}_{i} \quad \widehat{\Delta u}_{i}^{S_b}]$, and $\widehat{\Delta u}_{i}^{S_b}$ are the residuals obtained in the first stage for partition $S_b$.
Finally, since $\hat{\beta}_{1o,b} \in \hat{\theta}_{o,b}$ we estimate our parameter of interest $\beta_{1o}$ using:
$$\hat{\beta}_{1o}=\frac{1}{B}\sum_{b=1}^B \hat{\beta}_{1o,b}.$$
\begin{algorithm}
\caption{Super Learner Control Function Estimation: First-difference transformation }\label{Algortithm_CF1_ex}
\begin{algorithmic}
\STATE \textbf{Input:} Panel data sets $\{y_i,X_{i},z_{i}\}_{i\in \{1, 2, ..., N\}}$ from model satisfying assumptions \ref{A_data} to \ref{A_eps_u_omega}, a natural number B.
\STATE \textbf{Output:} An estimator of the parameter of interest $\beta_{1o}$.
\FOR{$ss \in \{1, 2, ..., SS\}$}
\STATE Split the individual index set $\{1, 2, ..., N\}$ into B sets $S_1, S_2,...,S_B$.
\FOR{$b \in \{1, 2, ..., B\}$}
\STATE Estimate the conditional expectation $\mathbb{E}[\Delta x_{1it}|\tilde{x}_{it},\tilde{x}_{it-1},z_{it},z_{it-1}]$ using a super learner and data correspondig to $S_b^c$. \STATE Predict $\mathbb{E}[\Delta x_{1it}|\tilde{x}_{it},\tilde{x}_{it-1},z_{it},z_{it-1}]$ for $i \in S_b$.
\STATE Obtain the residuals $\Delta u_{it}$ for $i \in S_b$ using the prediction of $\mathbb{E}[\Delta x_{1it}|\tilde{x}_{it},\tilde{x}_{it-1},z_{it},z_{it-1}]$ of previous step.
\STATE Estimate $\theta_{o} = [\beta_{1o}, \beta_{2o}, \rho_o]$ solving the sample moment conditions \ref{sample_moments} for $i \in S_b$.
\ENDFOR
\STATE Compute $\hat{\beta}_{1o,ss}=\frac{\sum_1^B \widehat{\beta_{1o,b,ss}}}{B}.$
\ENDFOR
\STATE Compute $\hat{\beta}_{1o}=\frac{\sum_1^{SS} \widehat{\beta_{1o,ss}}}{SS}.$
\end{algorithmic}
\end{algorithm}
\end{example_un}
\newpage
\subsection{The SLCF control function estimator is a debiased Machine Learning (DML) estimator }\label{S_Estimation}
In this subsection, we show that the SLCF estimator is a debiased Machine Learning (DML) estimator.
The reason is that the moment condition \ref{non_orthogonal_score} has a Neyman
orthogonal equivalent. Without loss of generality, we assume that there are not covariates such that the model \ref{dif_structural} simplifies to:
\begin{equation} \label{ControlFunction_eq_simple}
\tau y_{it}=\tau x_{1it}\beta_{1_o}+\rho_o \tau u_{it} + \tau \omega_{it}, i\in \{1, ..., N\}, t \in \{t_a, ..., T_i\}.
\end{equation}
If we stack up the the observations in matrix form, we obtain:
\begin{equation} \label{ControlFunction_eq_simple_matrix}
M_{\tau_i} y_i=M_{\tau_i} x_{1i}\beta_{1_o}+\rho_o M_{\tau_i} u_i + M_{\tau_i} \omega_i, i\in \{1, ..., N\}, t \in \{t_a, ..., T_i\}.
\end{equation}
Now, we premultiply \ref{ControlFunction_eq_simple_matrix} by the projection matrix in the orthogonal space of $\tau u_i $ ($A_i = I_i - M_{\tau_i} u_i( u_i'M_{\tau_i}'M_{\tau_i} u_i)^{-1} u_i'M_{\tau_i}'$), we obtain:
\begin{equation} \label{ControlFunction_eq_simple_matrix_transformed}
A_i M_{\tau_i} \tau y_i= A_i M_{\tau_i} x_{1i}\beta_{1_o}+ A_i M_{\tau_i} \tau \omega_i, i\in \{1, ..., N\}, t \in \{t_a, ..., T_i\}.
\end{equation}
Using the transformed model \ref{ControlFunction_eq_simple_matrix_transformed}, we know that the score for identification of $\beta_{1}$ is:
\begin{equation}\label{orthogonalScore}
\psi(O_i;\beta_1, \tau g) = x_{1i}'M_{\tau_i}'A_i V_i^{-1} A_i M_{\tau_i} \omega_i.
\end{equation}
\begin{proposition}
The score \ref{orthogonalScore} is Neyman Orthogonal (see appendix for the proof \ref{Proof_OrthogonalScore}).
\end{proposition}
\section{Large sample properties of the SLFC estimator}\label{S_Properties}
\begin{theorem} {Consistency of the SLCF estimator}\label{T_Consistency}
\noindent If $\{O_i = (y_i, X_i, z_i) \}_i\in [N]$ satisfy assumptions
\ref{A_data} to \ref{A_eps_u_omega}, assumptions of subsection
\ref{Annex_assumptions} in the Annex hold, i) $\Omega_i$ is positive
semi-definite matrix and $\Omega_i\mathbb{E}[\phi(O_i; \theta, \tau g_o)]=0$
only if $\theta = \theta_o$, ii) $\theta \in \Theta$ with $\Theta$
compact, iii) $\phi(O_i; \theta, \tau g_o)$ is continuous at each
$\theta \in \Theta$, iv) the function class $\mathcal{F}=\{\phi(O_i; \theta,
\tau g_o\}$ satisfies a uniform integrability condition
$\mathbb{E}[\operatorname*{sup_{ \theta \in \Theta}} \left \| \phi(O_i;
\theta, \tau g_o)\right \|] < \infty$, then $\hat{\theta}_o
\overset{p}{\to} \theta_o$.
\end{theorem}
The proof is presented in the appendix (Subsection \ref{Consistency_proof}).
\begin{theorem}{Asymptotic normality of the SLCF estimator of $\beta_{1o}$}\label{Theorem_AN}
\noindent If $\{O_i = (y_i, X_i, z_i) \}_i\in [N]$ satisfy assumptions \ref{A_data} to \ref{A_eps_u_omega}, assumptions of subsection \ref{Annex_assumptions} in the Annex hold,
then,
\begin{equation}
\sqrt{N_T} (\hat{\beta}_{1o} - \beta_{1o}) \xrightarrow{d} \mathcal{N}(0, \sigma_o^2),
\end{equation}
with $\sigma_o^2 = \mathbb{E}_P[x_{1i}'M_{\tau_i} A_i \tilde{V}_i^{-1}
A_i M_{\tau_i} x_{1i}]^{-1} \mathbb{E}_{P} [ x_{1i}' M_{\tau_i} A_i
\tilde{V}_i^{-1} A_i M_{\tau_i} V_i M_{\tau_i} A_i \tilde{V}_i^{-1}
A_i M_{\tau_i}x_{1i}] \mathbb{E}_P[x_{1i}'M_{\tau_i} A_i \tilde{V}_i^{-1} A_i M_{\tau_i} x_{1i}] ^{-1}$.
\end{theorem}
The proof is presented in the appendix (Subsection \ref{Proof_AsymptoticNormality}).
\begin{theorem}{Asymptotic normality of the SLCF estimator of $\theta_{o}$}\label{Theorem_AN_theta}
\noindent If $\{O_i = (y_i, X_i, z_i) \}_i\in [N]$ satisfy assumptions \ref{A_data} to \ref{A_eps_u_omega}, assumptions of subsection \ref{Annex_assumptions} in the Annex hold,
then,
\begin{equation}
\sqrt{N_T} (\hat{\theta}_{1o} - \theta_{1o}) \xrightarrow{d}
\mathcal{N}(0, \Sigma_o),
\end{equation}
with $\Sigma_o = \mathbb{E}_P[H_{i}' M_{\tau_i} \tilde{V}_i^{-1}
M_{\tau_i} H_i]^{-1} \mathbb{E}_{P} [ H_{i}' M_{\tau_i} \tilde{V}_i^{-1}
\omega_i \omega_i ' \tilde{V}_i^{-1} M_{\tau_i} H_{i}] \mathbb{E}_P[H_{i}'
M_{\tau_i}
\tilde{V}_i^{-1} M_{\tau_i} H_i]^{-1}$.
\end{theorem}
The proof is not presented as it follows a similar argument as the proof of Theorem \ref{Theorem_AN}.
\subsection{Estimator of the variance-covariance matrix}
The estimator of the variance-covariance matrix of the SLCF estimator of $\theta_o$ is given by:
\begin{equation}
\hat{\Sigma}_o = \hat{J}_o^{-1}
\frac{1}{B} \sum_{b=1}^B \mathbb{E}_{n,b} [ \hat{H'}_{i}^{S_b} \tilde{V}^{-1}_i M_{\tau_i}
(M_{\tau_i} y_{i}^{S_b}
- \hat{H}_{i}^{S_b} \hat{\theta}_o )(M_{\tau_i} y_{i}^{S_b} - M_{\tau_i}
\hat{H}_{i}^{S_b} \hat{\theta}_o )'\tilde{V}^{-1}_i M_{\tau_i} \hat{H}_{i}^{S_b} ] \hat{J}_o^{-1} ,
\end{equation}
where:
\begin{equation}
\hat{J}_o^{-1} = \frac{1}{B} \sum_{b = 1}^B \mathbb{E}_{n,b} [\hat{H'}_{i}^{S_b}\tilde{V}^{-1}_i\hat{H}_{i}^{S_b}].
\end{equation}
\subsubsection{Including sample splitting uncertainty}\label{Including_sample_splitting_uncertainty}
In order to take into account the sample splitting uncertainty, we can repeat the
cross-fitting process several times as explained in Section \ref{S_Estimation}
and algorithm \ref{Algortithm_CF1}. We also need to add correction term to
the variance-covariance estimator $\hat{\Sigma}_o$ given by:
\begin{equation}
C = \frac{1}{SS}\sum_{ss = 1}^{SS}(\hat{\theta}_{ss} - \hat{\theta}_o)(\hat{\theta}_{ss} - \hat{\theta}_o)',
\end{equation}
with $\hat{\theta}_o = \frac{1}{SS} \sum_{ss = 1}^{SS}\hat{\theta}_{ss}$.
\begin{algorithm}\label{Algortithm_VarCov}
\caption{Estimator of the variance-covariance of the Super Learner Control Function estimator}\label{your_label}
\begin{algorithmic}
\STATE \textbf{Input:} Panel data sets $\{y_i,X_{i},z_{i}\}_{i\in \{1, 2, ..., N\}}$
from model satisfying assumptions \ref{A_data} to \ref{A_eps_u_omega}, a natural number B.
The selected index set $\{1, 2, ..., N\}$ for each set $S_1, S_2,...,S_B$ across all samples splits (SS) obtained
in algorithm \ref{Algortithm_CF1}. The estimated $\hat{\theta}_{ss}$ obtained in algorithm \ref{Algortithm_CF1}.
\STATE \textbf{Output:} An estimator of the variance-covariance matrix of $\theta_{1o}$.
\FOR{$ss \in \{1, 2, ..., SS \}$}
\FOR{$b \in \{1, 2, ..., B\}$}
\STATE Estimate $\hat{J}_{o,b}^{-1} = \mathbb{E}_{n,b} [\hat{H'}_{i}^{S_b}\hat{H}_{i}^{S_b}]$.
\STATE Estimate $\mathbb{E}_{n,b} [ \hat{H'}_{i}^{S_b} M_{\tau_i}
(M_{\tau_i} y_{i}^{S_b}
- \hat{H}_{i}^{S_b} \hat{\theta}_{ss} )(M_{\tau_i} y_{i}^{S_b} - M_{\tau_i}
\hat{H}_{i}^{S_b} \hat{\theta}_{ss} )' M_{\tau_i} \hat{H}_{i}^{S_b}] $.
\ENDFOR
\STATE Compute $\hat{\Sigma}_s.$
\ENDFOR
\STATE Compute $C=\frac{1}{SS}\sum_{ss = 1}^{SS}(\hat{\theta}_{ss} -
\hat{\theta}_o)(\hat{\theta}_{ss} - \hat{\theta}_o)'.$
\STATE Compute $\hat{\Sigma}_o = \frac{1}{SS} \sum_{ss=1}^{SS}
(\hat{\Sigma}_s) + C.$
\end{algorithmic}
\end{algorithm}
\section{Comparison of the SLCFE with plug-in IV and naive plug-in 2SLS estimators}\label{S_Comparison}
\subsection{Naive plug-in 2SLS estimator is inconsistent and it is not numerically equivalent to SLCFE}
The Super Learner Control function estimator is not numerical equivalent to a 2SLS estimator that
plugs in the estimated transformed nuisance parameter $\tau g_o(\tilde{x}_{it}, z_{it})$.
Without loss of generality, we assume that exogenous covariates are not present (model
\ref{structural_equation_toy} to \ref{relationship_errors_toy}), and $V_{i} = I_{i}$.
As described in Section \ref{S_Identification}, one could be tempted to estimate the parameter of
insterest using the naive 2SLS estimator:
\begin{equation}
\hat{\beta}_{2SLS, 1_o} = \sum_b^B \frac{1}{B}(\sum_{i \in S_b} \widehat{\tau g_o}^{S_b^c}( z_i)' \widehat{\tau g_o}^{S_b^c}( z_i))^{-1} (\sum_{i \in S_b}\widehat{\tau g_o}^{S_b^c}( z_i)' M_{\tau_i} y_i).
\end{equation}
\begin{theorem}
If $\{O_i = (y_i, X_i, z_i) \}_i\in [N]$ satisfy assumptions
\ref{A_data} to \ref{A_eps_u_omega}, assumptions \ref{A_ThetaCompact} - \ref{A_several} of subsection
\ref{Annex_assumptions} in the Annex hold, Naive plug-in 2SLS estimation is not
numerically equivalent to SLCFE when $g_o(.)$ is nonlinear
\end{theorem}
\begin{proof}
Assuming that $B = 1$, we can write the estimators as:
\begin{equation}
\hat{\beta}_{2SLS, 1_o} = \Big(\sum_{i} \widehat{\tau g_o}( z_i)'\widehat{\tau g_o}^{S_b^c}( z_i)\Big)^{-1}
\Big(\sum_i\widehat{\tau g_o}( z_i)' M_{\tau_i} y_i\Big),
\end{equation}
\begin{equation}
\hat{\beta}_{SLCFE, 1_o} = \Big(\sum_i (M_{\tau_i}x_{1i})' \hat{A}_{i} M_{\tau_i} x_{1i} \Big)^{-1}
\Big(\sum_i (M_{\tau_i} x_{1i}) ' \hat{A}_i M_{\tau_i} y_i\Big),
\end{equation}
\noindent where $\hat{A}_i = I_i - \widehat{\tau u_i} ( \widehat{\tau u_i}'\widehat{\tau u_i})^{-1} \widehat{\tau u_i}'$.
\noindent Because $ \hat{A}_i M_{\tau_i} x_{1i} \neq \widehat{\tau g_o}( z_i)$ , the two estimators are not equivalent.
They are numerically equivalent if and only if $\tau g_o( z) \in \text{col}(\tau z)$.
\end{proof}
\begin{theorem}
If $\{O_i = (y_i, X_i, z_i) \}_i\in [N]$ satisfy assumptions
\ref{A_data} to \ref{A_eps_u_omega}, assumptions \ref{A_ThetaCompact} - \ref{A_several} of subsection
\ref{Annex_assumptions} in the Annex hold in particular that
$\left\| \tau g - \tau g_o \right\|_{P,2} \leq \delta_N$
with $\delta_N \sqrt{N_T} \to \infty$,
$\{\delta_N\}$ a sequence of positive numbers that converge
to 0 as $N \rightarrow \infty$
, the Naive plug-in 2SLS estimator is inconsistent.
\end{theorem}
\begin{proof}
The naive plug-in 2SLS estimator with cross-fitting is given by:
\begin{equation}
\hat{\beta}_{2SLS,1_o} = \frac{1}{B}\sum_b^B\Big(\sum_{i \in S^b} \widehat{\tau g_o}^{S_b^c}( z_i)'\widehat{\tau g_o}^{S_b^c}( z_i)\Big)^{-1}
\Big(\sum_{i \in S^b}\widehat{\tau g_o}^{S_b^c}( z_i)' M_{\tau_i} y_i\Big).
\end{equation}
Replacing $M_{\tau_i} y_i$ by $\tau g_o(z_i) \beta_{1_o} + M_{\tau_i}u_i \beta_{1_o} + M_{\tau_i}
\varepsilon_i$, adding and substracting $\widehat {\tau g_o}^{S_b^c}(z_i)\beta_{1_o}$, it can be written as:
\begin{equation}
\hat{\beta}_{2SLS,1_o} = \beta_{1_o} + \frac{1}{B}\sum_b^B\Big(\sum_{i \in S^b} \widehat{\tau g_o}^{S_b^c}( z_i)'\widehat{\tau g_o}^{S_b^c}( z_i)\Big)^{-1}
\Big(\sum_{i \in S^b}\widehat{\tau g_o}^{S_b^c}( z_i)'( \xi_i^{S_b^c} \beta_{1,o} + M_{\tau_i}u_i + M_{\tau_i} \varepsilon_i)\Big),
\end{equation}
with $\xi_i^{S_b^c} = \tau g_o(z_i) - \widehat{ \tau g_o} (z_i)^{S_b^c}$.
\noindent Now, adding and substracting $\tau g_o(z_i)$ to $\widehat{ \tau g_o}^{S_b^c} (z_i)$ we obtain that:
\begin{equation}
\hat{\beta}_{2SLS,1_o} = \beta_{1_o} + \frac{1}{B}\sum_b^B\Big(\sum_{i \in S^b} \widehat{\tau g_o}^{S_b^c}( z_i)'\widehat{\tau g_o}^{S_b^c}( z_i)\Big)^{-1}
\Big(\sum_{i \in S^b}(\xi^{S_b^c}_i + \tau g_o( z_i))'( \xi^{S_b^c}_i \beta_{1,o} + M_{\tau_i}u_i + M_{\tau_i} \varepsilon_i)\Big),
\end{equation}
The second term does not converge to zero at rate faster than $\sqrt{N_T}$ because $\tau g_o(z_i)$
is not centered at 0.
Thus, the naive plug-in estimator presents a first order bias term that does not vanish
using sample-splitting. Thus, the naive plug-in 2SLS estimator
is inconsistent.
\end{proof}
\subsection{Plug-in IV estimation is not numerically equivalent to SLCFE}
The IV estimator is given by:
\begin{equation}
\hat{\beta}_{IV, 1_o} = \sum_b^B \frac{1}{B}(\sum_{i \in S_b} \widehat{\tau g_o}^{S_b^c}( z_i)' M_{\tau_i} x_{1i})^{-1} (\sum_{i \in S_b}\widehat{\tau g_o}^{S_b^c}( z_i)' M_{\tau_i} y_i).
\end{equation}
\begin{theorem}
If $\{O_i = (y_i, X_i, z_i) \}_i\in [N]$ satisfy assumptions
\ref{A_data} to \ref{A_eps_u_omega}, assumptions of subsection
\ref{Annex_assumptions} in the Annex hold, plug-in IV estimation is not
numerically equivalent to SLCFE when $g_o(.)$ is nonlinear
\end{theorem}
\begin{proof}
Assuming that $B = 1$, we can write the estimators as:
\begin{equation}
\hat{\beta}_{IV, 1_o} = \Big(\sum_{i} \widehat{\tau g_o}( z_i)' M_{\tau_i} x_{1i}\Big)^{-1}
\Big(\sum_i\widehat{\tau g_o}( z_i)' M_{\tau_i} y_i\Big),
\end{equation}
\begin{equation}
\hat{\beta}_{SLCFE, 1_o} = \Big(\sum_i (M_{\tau_i}x_{1i})' \hat{A}_{i} M_{\tau_i} x_{1i} \Big)^{-1}
\Big(\sum_i (M_{\tau_i} x_{1i}) ' \hat{A}_i M_{\tau_i} y_i\Big),
\end{equation}
\noindent where $\hat{A}_i = I_i - \widehat{\tau u_i} ( \widehat{\tau u_i}'\widehat{\tau u_i})^{-1} \widehat{\tau u_i}'$.
\noindent Because $ \hat{A}_i M_{\tau_i} x_{1i} \neq \widehat{\tau g_o}( z_i)$ , the two estimators are not equivalent.
They are numerically equivalent if and only if $\tau g_o( z) \in \text{col}(\tau z)$.
\end{proof}
\begin{theorem}
If $\{O_i = (y_i, X_i, z_i) \}_i\in [N]$ satisfy assumptions
\ref{A_data} to \ref{A_eps_u_omega}, assumptions of subsection
\ref{Annex_assumptions} in the Annex hold in particular that
$\left\| \tau g - \tau g_o \right\|_{P,2} \leq \delta_N$
with $\delta_N \sqrt{N_T} \to \infty$,
$\{\delta_N\}$ a sequence of positive numbers that converge
to 0 as $N \rightarrow \infty$
, the plug-in IV estimator without the use of cross-fitting is inconsistent.
\end{theorem}
\begin{proof}
The plug-in IV estimator without cross-fitting is given by:
\begin{equation}
\hat{\beta}_{IV,1_o} = \Big(\sum_{i} \widehat{\tau g_o}( z_i)'M_{\tau_i} x_{1i}\Big)^{-1}
\Big(\sum_{i }\widehat{\tau g_o}( z_i)' M_{\tau_i} y_i\Big).
\end{equation}
Replacing $M_{\tau_i} y_i$ by $M_{\tau_i} x_{1i} \beta_{1_o} + M_{\tau_i}
\varepsilon_i$, it can be written as:
\begin{equation}
\hat{\beta}_{IV,1_o} = \beta_{1_o} + \Big(\sum_{i } \widehat{\tau g_o}( z_i)'M_{\tau_i} x_{1i}\Big)^{-1}
\Big(\sum_{i}\widehat{\tau g_o}( z_i)'( M_{\tau_i} \varepsilon_i)\Big),
\end{equation}
\noindent Now, adding and substracting $\tau g_o(z_i)$ to $\widehat{ \tau g_o} (z_i)$ we obtain that:
\begin{equation}
\hat{\beta}_{IV,1_o} = \beta_{1_o} + \Big(\sum_{i } \widehat{\tau g_o}( z_i)'M_{\tau_i} x_{1i}\Big)^{-1}
\Big(\sum_{i }(\xi_i + \tau g_o( z_i))'( M_{\tau_i} \varepsilon_i)\Big),
\end{equation}
The second term does not converge to zero at rate faster than $\sqrt{N_T}$ because $\tau g_o(z_i)$
is not learned using cross-fitting.
Thus, the plug-in IV estimator
without using sample-splitting presents a first order bias term that does not vanish. Thus, the plug-in IV estimator
is inconsistent when estimation of the transformed nuisance parameter is done without cross-fitting.
\end{proof}
\section{Monte Carlo simulation }\label{S_Simulation}
We test the performance of the proposed estimation methods using a Monte Carlo simulation experiment. As competing methods, we consider Within OLS, and Within 2SLS estimators including polynomial transformations of the covariates and instruments of degree 1 and 5.
\subsection{Settings}
\subsubsection{Linear structural equation and non-linear secondary equation }
We generate 100 samples from a triangular simultaneous model with a linear structural equation and a non-linear reduced-form equation with $N=1000$ and $T=2$. Both equations present additive disturbance terms, and individual-specific effects. The structural equation presents one endogenous regressor $x_{1it}$ and one exogenous covariate $x_{2it}$.
\begin{equation}\label{DGP1_structural_eq}
\begin{split}
y_{it}=\beta_1 x_{1it}+\beta_2 x_{2it}+\alpha_i+\epsilon_{it},
\end{split}
\end{equation}
\noindent where:
\begin{equation}\label{DGP1_exog}
x_{2it}=\alpha_i+\zeta_{it}, \quad with \quad \zeta_{it}\sim U(-2,2),
\end{equation}
\begin{equation}\label{DGP1_ivs}
z_{it}=\alpha_i+\nu_{it}, \quad with \quad \nu_{it} \sim U(-2,2),
\end{equation}
\begin{equation}\label{DGP1_reduced_form}
x_{1it}=g(x_{2it},z_{it})+\alpha_i+ u_{it},
\end{equation}
\begin{equation}
u_{it}\sim U(-1,1),
\end{equation}
\begin{equation}
\epsilon_{it}=0.9*u_{it}+\tilde{\zeta_{it}}, \quad with \quad \tilde{\zeta}_{it}\sim U(-1,1).
\end{equation}
\begin{equation}
\alpha_{i}\sim U(-1,1).
\end{equation}
\noindent The unknown function $g(x_{2it},z_{it})$ is similar to the ones proposed by \citet{GuoWP2022}, and depends on the parameter $a$ that controls the nonlinearity in the functional form and the strength of the linear correlation between $x_{1it}$ and $z_{it}$. As in \citet{GuoWP2022}, a larger $a$ causes higher nonlinearity and lower linear correlation between $x_{1it}$ and $z_{it}$. However, we cannot compare our results to theirs because they consider a model for cross-sectional data.
\begin{table}[h]
\addtolength\tabcolsep{2pt}
\caption{ Monte Carlo Experiment}
\addtolength\tabcolsep{2pt}
\resizebox{\textwidth}{!}{
\begin{tabular}{@{}l@{\hspace{25pt}}ccccc@{}}
\hline\noalign{\smallskip}
Scenario& \centering $g(\cdot)$ & $\beta_1$ &$dim(x_2)$&$\beta_2$ \\
\noalign{\smallskip}\hline\noalign{\smallskip}
1&$-a\cdot|z_{it}|-2tanh(x_{2it})+1/a*z_{it}$ & 1 & 1 & 1 \\\\
\noalign{\smallskip}\hline\noalign{\smallskip}
\end{tabular} }
\label{TableScenariosCF}
\end{table}
\begin{figure}[H]
\centering
\subfloat[\centering ]{\includegraphics[width=7cm]{FDxend_z_DGP_Buhlmann_1_a_10_.png}}
\qquad
\subfloat[\centering ] {{ \includegraphics[width=7cm]{FDxend_z_DGP_Buhlmann_1_a_1_.png}}}
\qquad
\caption{ a) Relationship between $\Delta x_{1it}$ and $z_{it}$ with $a=10$, b) Relationship between $\Delta x_{1it}$ and $z_{it}$ with $a=1$}
{\small }
\label{FunctionalFormDGP1}
\end{figure}
\subsection{ The results}
In this section, we present the results of the simulation experiment.
In figure \ref{ResultsAllaDGP1}, we present the average estimated parameter
of interest $\beta_{1o}$ for different values of $a$. The methods are Within
OLS (WOLS), Within 2SLS (W2SLS) using polynomial of degree 1 of the covariate and
the instrumental variables, Within 2SLS using polynomial of degree 5 of the
covariate and the instrumental variables (W2SLS_polynomial), the Super Learner
Control Function estimator on the First-Differenced data (FDSLCF), and the Super
Learner Control Function estimator on the Within-transformed data (WSLCF).
We also plot the p-value of the estimated control
functions using the first-differenced and the within transformed data. As $a$ is
larger, the nonlinearity in the reduced form equation is more important.
This translates in a lower correlation between the endogenous regressor and
the instrumental variables. As a result, the Within-2SLS estimators present
a decay in their performance. In contrast, the Super Learner Control Function
estimators remain stable around the true value of the parameter of interest
$\beta_{1o}=1$.
In panel (a) of figure \ref{GraphResultsDGP1}, we present the boxplot of
the estimated parameter of interest $\beta_{1o}$ for a highly nonlinear
reduced form equation (DGP 1 with $a=5$), and in panel (b) for a linear
reduced form equation (DGP 1 with $a=1$). The estimated parameters with
the Super Learner Control Function estimators are tightly concentrated
around the true value of the parameter of interest. In contrast, the
Within estimates present wide variation around the true value.
For the estimation of the nuisance parameters, we used a Super
Learner. The base learners used for the Super Learner are a linear model, a neural network, and the mean.
The architecture of the neural network is one hidden layer with two neurons, the activation function
is the sigmoid one, and the output function is logistic, the regularization is 0, and the maximum number of
iterations is 100.
\begin{remark}
Implementation: in practice, a richer library of base learners requires
a larger number of observations. In addition, it is better to avoid overfitting
by using regularization.
\end{remark}
\begin{figure}[H]
\centering
\includegraphics[width=15cm]{lineplot_a_DGP_Buhlmann_1_N_1000_T_2.png}
\caption{Average estimates of $\beta_{1o}$ for 100 samples simulated with different values of $a$, $N=1000$, $T=2$}
\label{ResultsAllaDGP1}
\end{figure}
\begin{figure}[H]
\centering
\subfloat[\centering ]{\includegraphics[width=7cm]{boxplot_a5_DGP_Buhlmann_1_N_1000_T_2_}}
\qquad
\subfloat[\centering ] {{ \includegraphics[width=7cm]{boxplot_a1_DGP_Buhlmann_1_N_1000_T_2_}}}
\caption{a) Estimates of $\beta_{1o}$ for a DGP highly nonlinear (a=5), b) Estimates of $\beta_{1o}$ for a DGP highly linear (a=1), $N=1000$, $T=2$}
{\small }
\label{GraphResultsDGP1}
\end{figure}
\begin{table}[h]
\addtolength\tabcolsep{2pt}
\caption{ Coverage results: Nominal coverage 95\%}
\addtolength\tabcolsep{2pt}
\resizebox{\textwidth}{!}{
\begin{tabular}{@{}l@{\hspace{25pt}}cccccccccccccc@{}}
\hline\noalign{\smallskip}
N & T & a & SS & K & $\beta_{1o}$& WOLS & W2SLS & W2SLS polynomial & FDCF
& WCF \\% & WCFOS & WCFOS2 \\
\hline\noalign{\smallskip}
1000 & 2 & 1 &10 & 5 & 1 & 100 & 100 & 100 & 71
& 76 \\% & 51 & 76\\
1000 & 2 & 2 &10 & 5 & 1 & 100 & 100 & 100 & 94
& 90 \\% & 51 & 90\\
1000 & 2 & 3 &10 & 5 & 1 & 100 & 100 & 100 & 99
& 99 \\%& 71 & 98 \\
1000 & 2 & 4 &10 & 5 & 1 & 100 & 100 & 100 &98
& 100 \\% &77 & 99\\
1000 & 2 & 5&10 & 5 & 1 & 100 & 100 & 100 & 97
&100 \\%&81 &99 \\
1000 & 2 & 6 &10 & 5 & 1 & 100 & 100 & 100 & 97
& 100 \\% &89 &99 \\
1000 & 2 & 7 &10 & 5 & 1 & 100 & 100 & 100 & 95
& 100 \\% & 91 & 98 \\
1000 & 2 & 8 &10 & 5 & 1 & 100 & 100 & 100 &97
&100\\% &92 &98 \\
1000 & 2 & 9 &10 & 5 & 1 & 100 & 100 & 100 & 96
& 100 \\%& 95 & 97\\
1000 & 2 & 10 &10 & 5 & 1 & 100 & 100 & 100 & 95
& 100\\% &95 &99 \\
1000 & 2 & 1 &100 & 5 & 1 & 100 & 100 & 100 & 79
&81 \\% &57 & 83 \\
1000 & 2 & 2 &100 & 5 & 1 & 100 & 100 & 100 & 95
& 96 \\%& 57 & 96 \\
1000 & 2 & 3 &100 & 5 & 1 & 100 & 100 & 100 & 95
& 96\\%& 54 & 96\\
1000 & 2 & 4 &100 & 5 & 1 & 100 & 100 & 100 &96
& 97 \\%& 58 &97 \\
1000 & 2 & 5 &100 & 5 & 1 & 100 & 100 & 100 & 99
&100 \\% &100 &100 \\
1000 & 2 & 6 &100 & 5 & 1 & 100 & 100 & 100 & 96
& 100 \\%& 92 & 100\\
1000 & 2 & 7 &100 & 5 & 1 & 100 & 100 & 100 & 97
&100 \\%&94 &100 \\
1000 & 2 & 8 &100 & 5 & 1 & 100 & 100 & 100 & 99
& 100 \\% & 94 &99 \\
1000 & 2 & 9 &100 & 5 & 1 & 100 & 100 & 100 &100
& 100 \\%& 92 &100 \\
1000 & 2 & 10 &100 & 5 & 1 & 100 & 100 & 100 &99
&100 \\%&94 &100 \\
\noalign{\smallskip}\hline\noalign{\smallskip}
\end{tabular} }
\label{TableScenariosCF}
\end{table}
\section{Empirical Application: Air Pollution and Educational Outcomes in the U.S. (Revisited)} \label{EA}
As an illustration, we estimate the causal effect of air pollution on student
performance in the United States. The data on air pollution are obtained from
\cite{DeryuginaAER2019}. This dataset presents daily data on air pollutants:
PM 2.5, ozone ($O_3$), carbon monoxide (CO), sulfur dioxide ($SO_2$),
nitrogen dioxide ($NO_2$). It also contains data on wind direction,
wind speed, and weather conditions at county level.
We aggregate the data to yearly level by averaging across days.
The data on student performance are obtained from Stanford Educational
Data Archive (SEDA, \citealt{FahleSEDA2024}) for county-yearly level.
We merge both datasets by county-year indicator. The final dataset is
a panel data set with 787 counties for the period 2009 to 2013. We transform the
data using the first-difference and within transformations.
We specify the following linear regression for county $d$, subject $s$,
cohort $c$ at year $t$:
\begin{equation}
y_{scdt}=\beta PM2.5_{dt}+X_{dt}'\gamma+\psi_s+\phi_c+\alpha_{d} +\tilde{\alpha}_g+\epsilon_{dst},
\end{equation}
\noindent where $y_{dst}$ represents the average grade of subject $s$,
cohort $c$, county $d$, subject $s$ at year $t$. $PM2.5_{dt}$ represents the
yearly average of $PM2.5$ concentration at county $d$ in year $t$, $X_{dt}$
is a vector of control variables at county-year level including average school
characteristics, social characteristics, and weather variables,
$\phi_c$ corresponds to cohort specific effects, $\tilde{\alpha}_g$ grade specific effects, and
$\psi_s$ corresponds to subject specific effects. $\alpha_{d}$
represents county fixed effects.
As mentioned before, $PM2.5_{dt}$ is endogenous. Thus, following
\cite{DeryuginaAER2019} we use wind direction as an instrument of
air pollution along with a nonlinear first-stage equation as follows:
\begin{equation} PM25_{dt}=h(WD_{dt},X_{dt},\psi_s, \phi_c, \tilde{\alpha}_g)
+\eta_{d}+\varepsilon_{dt},
\end{equation}
\noindent with $WD_{dt}$ is the yearly average wind direction in county
$d$ in year $t$, $\eta_{d}$ are county specific effects.
Using a First-difference transformation, we estimate that an increase in the
concentration of PM 2.5 causes a reduction in test scores of 0.976
standard deviations that is statistically significant at a 5\% significance level.
The base learners used for the Super Learner are the mean, a linear model, a neural network, and a random forest.
The architecture of the neural network is one hidden layer with two neurons, the activation function
is the sigmoid one, and the output function is linear, the regularization is 0, and the maximum number of
iterations is 100. The random forest is composed of 100 trees with minimum leaf size of 5, and two
variables to select at each split. The number of sample splits is equal to 10,
and the number of folds is 5. We use $\tilde{V}_i = I_i$ which is the most efficient in the case of the
Within transformation. Inference is still valid as we used a sandwich type estimator
of the variance-covariance matrix of the coefficients of the structural model.
\begin{table}[H]
\centering
\caption{Estimated effect of air pollution (PM 2.5 concentration) on student performance}
\begin{tabular}{cccccc}
\hline
WOLS & W2SLS & FD SLCF & W SLCF \\
\hline
-0.072 & -0.844 & -0.979 & -0.859 \\
(0.011) & (0.762) & (0.101) & (0.085) \\
\hline
\end{tabular}
\caption*{Note: US counties, (N = 787). }\label{tab:my_label}
\end{table}
\section{Conclusions}\label{S_Conclusions}
This paper proposes a triangular simultaneous equation model for panel data
with additive separable individual-specific fixed effects composed of a
linear structural equation with a nonlinear reduced form equation to deal with
the issue of weak instrumental variables due to nonlinearities in panel data
settings. The parameter of interest is the structural parameter of the endogenous
variable. The identification of this parameter is obtained under the assumption
of available exclusion restrictions and using a control function approach.
Estimating the parameter of interest is done using an estimator that we call
Super Learner Control Function estimator (SLCFE). The estimation procedure
is composed of two main steps and cross-fitting.
We estimate the control function with a super learner using sample
splitting. In the following step, we use the estimated control function
to control for endogeneity in the structural equation.
Cross-fitting is done across the individual dimension.
The estimator is consistent and asymptotically normal achieving a
parametric rate of convergence.
We perform a Monte Carlo simulation to test the performance of
the estimators proposed. We conclude that the estimator perform
well, provided that we can learn accurately the nuisance parameter
in the first stage.