Fixed Effects as Generated Regressors

Fixed effects as generated regressors

Introduction

In this paper, we will consider a model in which the researcher is interested in a parameter $\mu_0$ that can be characterized by a set of moment conditions:

where $W_i$ will be observables that are specific to an individual $i$, while $\alpha_i$ is unobservable to the researcher.

If the researcher has no information about $\alpha_i$, one can derive observable implications from ((ref)) which typically results in partial identification of $\mu_0$ (see entropic,ot,li2021). However, oftentimes the researcher has a panel data model that specifies $\alpha_i$ as fixed effects. Formally, for each individual $i$, the researcher observes $(Y_{it},X_{it}')_{t=1}^T$ that are related to $\alpha_i$ through the linear panel data model:

with $u_{it}$ satisfying sequential exogeneity to accommodate dynamic panel data models.\footnote{Models with additional time fixed effects can be reduced to the setup by first taking the within transformation for each time period (subtracting the mean of each variable within each time period. See Section (ref) for an example.} We do not impose any independence assumption between $(Y_{it},X_{it}')_{t=1}^T$ and $W_i$.

Many interesting economic problems share the above structure. Typical examples include value-added for teachers chetty2014a and the match quality between firm-product pairs kim2024. More concretely, chetty2014 and jackson study whether students who are matched with teachers with higher value-added in primary school might also have better long-run outcomes in college attendance and earnings. The latent variable $\alpha_i$ represents the value-added on students' test scores for teacher $i$, and $\mu_0$ is the slope coefficient of a linear regression of students' average (for each teacher $i$) long-run outcomes on $\alpha_i$'s. The linear panel data model would be specified using $Y_{it}$ as the test score for student $t$ who has teacher $i$, and $X_{it}$'s are covariates that include lagged test scores and socio-economic status of the student. In the context of firm dynamics, kim2024 show that the firm-product match quality is an important predictor for firm's product adding and dropping behaviors. In their specification, $\alpha_i$ is the firm-product match quality where $i$ represents a firm-product pair and $t$ represents time. The panel data model would use product-firm specific value of shipment as $Y_{it}$, product tenure and firm size as the covariates $X_{it}$. The parameter of interest $\mu_0$ is the average marginal effect of a shift in product-firm match quality on the probability that the firm adds or drops that product.

Using the panel data model, one can easily recover an estimate $\hat{\alpha}_i$ of $\alpha_i$, but directly plugging the estimated $\hat{\alpha}_i$ into the moment condition ((ref)) and ignoring the potential measurement error in $\hat{\alpha}_i$ might render the inference results invalid. This paper aims to provide an inference procedure for $\mu_0$ that takes into account the error in estimating $\alpha_i$'s properly. Before describing the approach, we first consider a special case of ((ref)) which has drawn a lot of attention in the literature CGK2025,xie2025,battaglia2024,deeb2021.

The teacher value-added example in chetty2014 follows the above structure where $W_i$ is the average long-run outcome of the students who have teacher $i$ while $\alpha_i$ is teacher $i$'s value-added.

The need to estimate $\alpha_i$'s and to use the estimates in the moment conditions ((ref)) suggests having small average mean squared error (MSE) $\frac{1}{N} \sum_{i=1}^N \mathbb{E}(\hat{\alpha}_i-\alpha_i)^2$ could be important for estimating $\mu_0$. Recognizing this fact, Empirical Bayes (EB) corrections of $\hat{\alpha}_i$'s are suggested to alleviate the potential impact of the measurement error angrist2017,CGK2025. Results using EB corrected versions of $\alpha_i$'s typically abstract away from the panel data structure that defines $\alpha_i$ and start by assuming that an estimator $\hat{\alpha}_i$ of $\alpha_i$ is available which has a normal distribution with mean $\alpha_i$ and some variance $\sigma_i^2$ xie2025,CGK2025. The EB corrections are then applied to $\hat{\alpha}_i$'s. xie2025 establishes a central limit theorem where estimates of $\mu_0$ are obtained by replacing $\alpha_i$'s in the moments by the EB corrected estimates in Example (ref). However, the above result requires i) independence between $u_{it}$ and $v_i$; ii) $\alpha_i$ being independent of $\mbox{Var}_i(\sum_{t=1}^T u_{it})$. CGK2025 recognizes that the latter independence assumption is restrictive. Instead, they advocate combining the classical correction of error-in-variables regression in deaton1985 with bootstrapped confidence interval, which they showed is a valid inference procedure that relaxes the independence between $\alpha_i$ and $\mbox{Var}_i(\sum_{t=1}^T u_{it})$. However, they still need to assume orthogonality condition between $u_{it}$ and $v_i$, and their approach is hard to extend to nonlinear moments.

The required independence assumption $u_{it}$ and $v_i$ may not hold in some applications of interest. For example, in the teacher value-added example, if $W_i$ represents average college attendance rate of the students, it can be correlated with $u_{it}$ for the same set of students. In the firm-product match quality example, an idiosyncratic shock in demand could affect the product value of shipment and also induce the firm to drop the product.

To address this issue, we instead exploit the panel data structure that defines $\alpha_i$'s. Although it can be more restrictive in some settings where estimators of $\alpha_i$ does not come from a panel data model, it enables us to relax the exogeneity assumption between $u_{it}$ and $v_i$. We borrow insights from weidner2024 and consider asymptotics in which both the dimension of individuals $N$ and time periods of panel data $T$ go to infinity. This allows us to use debiased machine learning (DML) techniques in cher2022 to construct orthogonal moment conditions that eliminate any first order effect of the measurement error in $\hat{\alpha}_i$ as long as $T$ goes to infinity sufficiently fast compared with $N$. More specifically, we require $\lim_{N\rightarrow \infty} \frac{N^{\frac{1}{2}}}{T}=0$ and $\lim_{N\rightarrow \infty} \frac{N}{T}=\infty$. As long as the above limits hold, the asymptotic distribution of the estimator $\hat{\mu}$ using the orthogonalized moments will not depend on the specific rate at which $T$ goes to infinity relative to $N$. Constructing the orthogonal moment conditions requires the sequential exogeneity assumption to derive the adjustment term to the moment function as shown in ichimura. This shows that taking the panel data structure into account can offer potential improvements on the estimator $\hat{\mu}$. Furthermore, shrinkage methods like EB and Stein's Unbiased Risk Estimation (SURE) can be readily applied to estimates of $\alpha_i$, because with orthogonal moments how $\alpha_i$'s are estimated does not affect the asymptotic distribution.\footnote{SURE has been studied in xie2012, and for applications in econometrics, see kwon.} We further derive a central limit theorem that characterizes the limit distribution of $\hat{\mu}$ based on orthogonal moments. The result is the first that unites the cross-sectional moments and linear panel data model without imposing parametric assumptions on the panel data part of the model. \footnote{For the case of the panel data model being parametric, see weidner2024. }

One stark difference of our approach from the standard DML literature is in the procedure of cross-fitting. DML requires the nuisance parameter $\alpha_i$ to be learned using samples from other folds and then use it to construct the sample moments $m(W_i,\hat{\alpha}_i,\mu)$ for the current fold. However, this is not possible in our context, because only data about individual $i$ contains information about $\alpha_i$. Hence any reasonable estimator of $\alpha_i$ must depend on the data for individual $i$. We adapt the cross-fitting procedure by first estimating $\beta_0$ in the panel data model using samples from other folds (denote the estimator as $\tilde{\beta}_l$), and then an estimator $\tilde{\alpha}_{il}$ of $\alpha_i$ is constructed as the difference between averages of $Y_{it}$ over time and averages of $X_{it}'\tilde{\beta}_l$ over time. The difference $\tilde{\alpha}_{il}-\alpha_i$ depends on $i$-th sample only through the panel data residuals $u_{it}$. We showed that with suitable controls on the correlation between $u_{it}$ and $\frac{\partial m(W_i,\alpha_i,\mu_0)}{\partial \alpha_i}$, using $i$-th sample to estimate both $\alpha_i$ and $\mu_0$ will not affect the limiting distribution of the orthogonal moments (see condition (ii) of Assumption (ref) for a precise statement).

In a simulation study with linear dynamic panel and linear regression model for ((ref)), we showed that when $u_{it}$ and $W_i$ are correlated, the approaches used in CGK2025 and xie2025 that ignore the potential endogeneity issue have non-negligible bias and consequently larger root mean squared error (RMSE). In addition, standard test statistics based on their estimators do not properly control the null rejection probabilities. In contrast, estimators that employ orthogonal moment conditions have null rejection probabilities closer to the nominal size, and employing an EB (or SURE) corrections on the $\hat{\alpha}_i$'s will perform slightly better in terms of null rejection probabilities but may incur larger RMSE. However, using orthogonal moment conditions results in larger variance in the estimator in the simulation, and can lose power compared with CGK2025 if endogeneity concerns are not present.\footnote{The estimator used in CGK2025 achieves the semi-parametric efficiency bound under their assumptions.}

To demonstrate the power of the method in dealing with non-linear moment conditions, we further consider estimating coefficients of a logit model as an empirical illustration. The moment conditions in this context set the expectation of the score function to zero. Specifically, we focus on the experimental site selection of the agricultural catastrophe insurance (ACI) program in China. policy have documented consistent patterns of positive site selection in the universe of policy experiments conducted by Chinese central government, but they did not commit to a particular explanation of the positive selections. As implementation details are relegated to local government and good practices discovered during the experimentation phase are encouraged to be rolled out across the whole country, the central government has incentive to pick experimental sites that are more specialized in the specific policy domain and could serve as “model sites”. We test the hypothesis in the context of ACI, where the key variable of interest ($\alpha_i$) is the level of specialization in grain production for each county. As the specialization level is unobserved to the researcher, it is captured as fixed effects in a linear panel data model with $Y_{it}$ being the grain production per rural employment. We then focus on the parameter associated with $\alpha_i$ that enters a linear logit model where the response is whether the county is selected as the experimental site. Results show that the coefficient associated with the level of specialization is positive and significant, but using orthogonal moments reduce the magnitude of the coefficient from around 0.56 (using simple plug-in methods) to around 0.4. EB and SURE corrections have more limited impacts on the estimates. Therefore, the evidence favors the theory that governments are choosing sites that are more experienced in agricultural production, and calls for usage of orthogonal moments when endogeneity between $u_{it}$ and $W_i$ is of concern.

Literature

This paper is related to literature on correcting measurement error for estimated latent variables. In addition to the classical approach like error-in-variables regression as in deaton1985 and the generated regressor problem in pagan, more recent contributions include CGK2025, xie2025, battaglia2024 and deeb2021. battaglia2024 uses an asymptotic regime so that the measurement error is similar in magnitude to the sampling error of the model ((ref)) and derives formulas for bias correction of the estimator. They focused on the case where $\alpha_i$'s are estimated using Machine Learning Methods. deeb2021 focuses on the teacher value-added context and derives moment conditions that can help estimate $\mu_0$ in Example (ref). The results of these papers rely on the linear regression structure that defines $\mu_0$, and generalizations to more general models are not straightforward.

The long panel data asymptotics with both $N$ and $T$ going to infinity has been studied in hahn2002. See also hahn2004 for the non-linear panel data context. The main objective is to derive bias correction methods for $\beta_0$ to combat the incidental parameter problem of estimating infinitely many $\alpha_i$'s neyman. Other approaches to circumvent the incidental parameter problem include reparametrization and integrating out the $\alpha_i$'s in a Bayesian framework lancaster.

This paper also uses techniques from the DML literature as in dml,cher2022. The paper that is most closely related is weidner2024, which also features a panel data first stage to recover latent variables that is going to be used in a model like ((ref)). However, their focus is on non-linear panel data models that are fully parametric and enables the researcher to exploit the likelihood function and construct higher-order orthogonal moments. This paper extends their approach to a semi-parametric context without making any distributional assumptions on $u_{it}$. The construction of the orthogonal moment functions and estimation procedure follows from cher2022 and ichimura.

The paper is also related to the literature on EB corrections for a parallel decision problem. Examples that use Empirical Bayes to correct for potential measurement error include angrist2017. Standard EB corrections may not work well when the error distribution is far from normal. Recent literature suggests choosing the shrinkage parameter based on the idea of Stein's Unbiased Risk Estimation (SURE) kwon,xie2012,brown. In our framework, any corrections using EB or SURE on the estimates of $\alpha_i$ can be incorporated, will not affect the asymptotic distribution of $\hat{\mu}$, and could be beneficial to the performance of the estimator in finite sample.

The plan of the paper is as follows. Section 2 formally introduces the model and the concept of Neyman orthogonality. Section 3 describes the orthogonal moment conditions and the cross-fitting procedure, while Section 4 gives the assumptions that establishes a central limit theorem for the estimator proposed in Section 3. Section 5 conducts a simulation study of the estimator based on orthogonal moments compared with other estimators in the literature. The empirical illustration of experimental site selection is detailed in section 6. All proofs are relegated to the appendices.

Individual fixed effects in cross-sectional moments

We will consider a long-panel data setting in which the cross-sectional dimension is indexed by $i$ and the time dimension is indexed by $t$. Suppose we have a balanced panel with $N$ number of individuals observed across $T$ time periods. For each individual $i$ at time $t$, the researcher observes a vector of observables $Z_{it}=(Y_{it},X_{it}')$ and for each individual $i$, the researcher observes a vector $W_i$ that only varies at the cross-sectional level. Assume that in the first-stage, the following fixed effect model is estimated to recover the set of fixed effects $\alpha_i$ for each individual $i$.\footnote{We will not restrict to a particular estimator for the first-stage model, but provide sufficient conditions on the estimates $\beta$ and $\alpha_i$ that can eliminate the estimation bias in the second stage.}

where $Y_{it}\in \mathbb{R}$ is the scalar outcome variable, and $X_{it}\in \mathbb{R}^p$ is a vector of observables.

One approach in dealing with the moment conditions is to use the plug-in principle by replacing $\alpha_i$'s with their corresponding estimates $\hat{\alpha}_i$, which gives us the sample analogue estimating equation for $\hat{\mu}$:

\[\sum_{i=1}^N m(W_i,\hat{\alpha}_i,\hat{\mu})=0\] under Assumption (ref).

Assuming the moment conditions $m(\cdot)$ are sufficiently smooth, a standard Taylor expansion of the above moment condition around $(\hat{\alpha}_i,\mu_0)$ yields: \[-\frac{1}{N} \sum_{i=1}^N m(W_{i},\hat{\alpha}_i,\hat{\mu})=\mathbb{E}\left( \frac{\partial m(W_i,\hat{\alpha}_i,\mu^*)}{\partial \mu})+o_p(1)\right) (\hat{\mu}-\mu_0), \] for some intermediate value $\mu^*$. The RHS of the above equation is standard in M-estimation, and the main difficulty lies in establishing the asymptotic distribution of the LHS.

Following insights from weidner2024, we can do a further Taylor Expansion for the LHS of the above display (for simplicity drop the minus sign):

under Assumption (ref) and also assuming the Hessian of $m(W_i,\alpha_i,\mu_0)$ with respect to $\alpha_i$'s is uniformly bounded almost surely.

The first term is the sum of moment conditions evaluated at the true parameter values, and a standard central limit theorem applies when we rescale the whole expression by $\sqrt{N}$. In standard DML setting, the term $II$ in ((ref)) will be $o_p(1)$ by using sample splitting and cross-fitting method. For example, by estimating $\hat{\alpha}_i$ on half of the sample and then estimating $\mu_0$ using the other half of the sample. However, in the current context, $\hat{\alpha}_i$ has to depend on the data for individual $i$ and hence could be correlated with $\frac{\partial m(W_i,\alpha_i,\mu_0)}{\partial \alpha_i }-\mathbb{E}(\frac{\partial m(W_i,\alpha_i,\mu_0)}{\partial \alpha_i })$. Therefore, standard cross-fitting techniques do not apply in the current context and they need to be adapted to make $II$ negligible asymptotically. We will show how to do this in the next section. The term $III$ will be zero if $\mathbb{E}(\frac{\partial m(W_i,\alpha_i,\mu_0)}{\partial \alpha_i })=0$. This condition is called Neyman Orthogonality and is widely used in the literature newey1994,dml,weidner2024. Finally, the term $O_p(\sup_i|\hat{\alpha}_i-\alpha_i|^2)=O_p(\frac{1}{T})$, and when $T$ increases sufficiently fast with $N$, it will not affect the asymptotic distribution of $\hat{\mu}$. We now formally introduce the definition of Neyman Orthogonality:

Inference procedure

Constructing Orthogonal Moment Conditions

As shown in ((ref)), Neyman orthogonality is useful in controlling $\frac{1}{\sqrt{N}} \sum_{i=1}^N m(W_i,\hat{\alpha}_i,\mu_0)$ when the true $\alpha_i$'s are replaced by their corresponding estimates. Here we will address the question of how to modify moment conditions $\mathbb{E}(m(W_i,\alpha_i,\mu_0))=0$ to achieve Neyman Orthogonality, but still preserve identification of $\mu_0$.

Before entering into the details, first note that ((ref)) implies the following moment restrictions on $\alpha_i$:

for any function $b(X_i,\alpha_i)$ with $X_i=(X_{i1}',\cdots,X_{iT}')'$.\footnote{One can also exploit sequential exogeneity by choosing a different period $t'$ than the final period $T$ and restricting the set of $X_i'$ to be the $X_{it}'$ up to period $t'-1$. This corresponds to the endogenous orthogonality conditions in ichimura which requires additional regularity conditions. For simplicity, we do not pursue further in that direction.}

One way to construct orthogonal moment conditions is through modifying the original moment conditions by adding the influence function of the original moment condition with respect to $\alpha_i$ as an adjustment term like cher2022, ichimura and weidner2024. More specifically, we can choose a function $a(X_i,\alpha_i,\mu)$ and form a modified moment condition $m^*(W_i,\alpha_i,Z_i,\beta,\mu,a)$ as:\footnote{Here $Z_i$ denotes $(Z_{i1}',\cdots,Z_{iT}')'$, with $(Y_{it},X_{it}')$ being a component of $Z_{it}$. } \[m^*(W_i,\alpha_i,Z_i,\beta_0,\mu,a)=m(W_i,\alpha_i,\mu)+a(X_i,\alpha_i,\mu) (Y_{iT}-X_{iT}'\beta_0-\alpha_i).\]

By the sequential exogeneity assumption, we have for any $\mu$,

hence the adjustment term will always have mean zero. Therefore, $\mu$ satisfies \\ $\mathbb{E}[m^*(W_i,\alpha_i,Z_i,\beta_0,\mu,a)]=0$ if and only if $\mathbb{E}[m(W_i,\alpha_i,\mu]=0$. As long as $\mu_0$ can be identified by the original moment conditions $\mathbb{E}[m(W_i,Z_i,\mu)]=0$, then it can also be identified by the modified ones $\mathbb{E}[m^*(W_i,\alpha_i,Z_i,\beta_0,\mu,a)]=0$. This is because $\mathbb{E}[m^*(W_i,\alpha_i,Z_i,\beta_0,\mu,a)]=\mathbb{E}[m(W_i,\alpha_i,\mu]$ for any $\mu$ and $a(X_i,\alpha_i,\mu)$.

Differentiating with respect to $\alpha_i$, and then taking conditional expectation with respect to $(X_i,\alpha_i)$, we have

The last step follows by sequential exogeneity assumption in ((ref)). The Neyman Orthogonality Condition can be achieved by choosing $a(X_i,\alpha_i,\mu_0)$ so that the above display equals zero:\footnote{The fact that $\mathbb{E}(m^*(W_i,\alpha_i,Z_i,\beta_0,\mu_0,a))=0$ can be easily checked via law of iterated expectation by conditioning on $(X_i,\alpha_i)$ first.}

\[a(X_i,\alpha_i,\mu_0)= \mathbb{E}\left( \frac{\partial m(W_i,\alpha_i,\mu_0)}{\partial \alpha_i}|X_i,\alpha_i \right)\]

In addition, \[\mathbb{E}\left( \frac{\partial m^*(W_i,\alpha_i,Z_i,\beta_0,\mu_0,a)}{\partial a} \right) =-\mathbb{E}(Y_{iT}-X_{iT}'\beta-\alpha_i)=0\] so the new moments $m^*(W_i,\alpha_i,Z_i,\beta_0,\mu_0,a)$ is also orthogonal to the new nuisance function $a(X_i,\alpha_i,\mu_0)$. With some abuse of notation, from now on, $a(X_i,\alpha_i,\mu_0)$ is used exclusively to denote $\mathbb{E}\left( \frac{\partial m(W_i,\alpha_i,\mu_0)}{\partial \alpha_i}|X_i,\alpha_i \right)$.

This is the strategy proposed and implemented for nuisance parameters that are typically non-parametric functions of the DGP (See e.g. dml,cher2022,ichimura) rather than fixed effect estimates. The component $X_i$ in $a(\cdot)$ concatenates $X_{it}$'s which will be a high-dimensional object under the data generating process of $N,T\rightarrow \infty$.\footnote{weidner2024 use the above construction but in their applications, $W_i$ only involves $X_i$. Consequently, $\mathbb{E}\left( \frac{\partial m(W_i,\alpha_i,\mu_0)}{\partial \alpha_i}|X_i,\alpha_i \right)=\frac{\partial m(W_i,\alpha_i,\mu_0)}{\partial \alpha_i}$, which is a significant simplification. } Some sparsity restriction is needed for the conditional expectation $\mathbb{E}\left( \frac{\partial m(W_i,\alpha_i,\mu_0)}{\partial \alpha_i}|X_i,\alpha_i\right)$ to be estimated with reasonable precision by employing standard high-dimensional regression tools like Elastic Net.

Example 1 (continued)

From the moment conditions in ((ref)), let $m_1(W_i,\alpha_i,\mu_0)=W_i -\mu_{01}-\mu_{02} \alpha_i$, and $m_2(W_i,\alpha_i,\mu)=\alpha_i(W_i-\mu_{01}-\mu_{02} \alpha_i)$ be the two moment functions. By taking derivatives, we see that the $a_1(X_i,\alpha_i,\mu_0)= -\mu_{02}$ and $a_2(X_i,\alpha_i,\mu_0) = \mathbb{E}(v_i - \mu_{02} \alpha_i|X_i,\alpha_i)$.

Implementation

Constructing the orthogonal moment $m^*(W_i,\alpha_i,Z_i,\beta_0,\mu_0,a)$ requires estimates of $\alpha_i$, $\beta_0$, and $a(X_i,\alpha_i,\mu_0)$ which in turn needs a preliminary estimator of $\mu_0$. In order to avoid the recursiveness in estimating $\mu_0$ and control the bias due to using the same data to estimate both the nuisance parameters and the parameter of interest $\mu_0$, cross-fitting procedures are widely adopted in the literature dml,cher2022. Although standard cross-fitting cannot deal with the estimation error in $\alpha_i$, it is still useful to control the error in the estimation of $a(X_i,\alpha_i,\mu_0)$ (see Assumption (ref) below). In particular, sample splitting and cross-fitting enable us to form an estimate of the function $a(X_i,\alpha_i,\mu_0)$ without using observation $i$. Analysis of the plug-in orthogonal moments $m^*(W_i,\hat{\alpha}_i,Z_i,\hat{\beta},\mu,\hat{a})$ can then be performed by conditioning on the observations used to estimate $a(X_i,\alpha_i,\mu_0)$, and hence one can treat $\hat{a}(X_i,\alpha_i,\mu_0)$ as a fixed function and does not need to account for the effect of individual $i$ on estimating $a(X_i,\alpha_i,\mu_0)$. Adapting the standard cross-fitting procedure to our context gives the following version:

Before presenting the conditions that guarantee the validity of algorithm (ref), the following remarks are helpful for implementation.

Asymptotic Theory

In this section, we give some high-level conditions that guarantee that the estimator $\hat{\mu}$ proposed in Section (ref) will be $\sqrt{N}$-asymptotically normal under the asymptotics that both $N,T \rightarrow \infty$, and $\hat{V}$ will be a consistent estimator of the asymptotic variance. Importantly, as in cher2022, I will first give conditions such that the orthogonal moment $\hat{m}^*(\mu_0)$ defined in ((ref)) satisfies

Condition (i) imposes a rate of convergence for the parameter $\beta$ in the first-stage panel data regression. The rate of convergence $O_p(\frac{1}{\sqrt{NT}})$ is standard for the panel data setting. Condition (ii) is a high-level convergence rate for prediction error of the function $a(X_i,\alpha_i,\mu_0)$. The restriction $\zeta\geq 1-g$ is needed to control $(\tilde{\alpha}_{il}-\alpha_i)(\tilde{a}_{il}-a(X_i,\alpha_i,\mu_0)$. As $\tilde{\alpha}_{il}-\alpha_i = O_p(\frac{1}{\sqrt{T}})=O_p(N^{-\frac{g}{2}})$, the condition $\zeta\geq 1-g$ ensures that $(\tilde{\alpha}_{il}-\alpha_i)(\tilde{a}_{il}-a(X_i,\alpha_i,\mu_0)=o_p(\frac{1}{\sqrt{N}})$. More primitive conditions can be specified depending on different types of learning algorithms that are used to learn the function $a(X_i,\alpha_i,\mu_0)$. Condition (iii) focuses attention on functions that are square-integrable which is standard in non-parametric regressions. Condition (iv) is implied by condition (ii) by observing that conditional on $\mathcal{W}_l^C$, we have

\[\frac{1}{|I_l|} \sum_{i\in I_l} [(\tilde{a}_{il}-a(X_i,\alpha_i,\mu_0))^2-\mathbb{E}((\tilde{a}_{il}-a(X_i,\alpha_i,\mu_0))^2|\mathcal{W}_l^C)]=O_p(N^{-1/2})=o(N^{-\zeta}),\] where the first equality observes that it is an average of iid zero mean random variables and uses lemma 6.1 of dml to convert conditional convergence to unconditional convergence. The second one uses $\zeta<1/2$ in Assumption (ref).

Condition (i) of Assumption (ref) imposes a rate of convergence of the loss for the parallel estimation problem for all $\alpha_i$'s with $i\in I_l$. Although asymptotically using the original $\tilde{\alpha}_{il}$ based on averages over time leads to the same asymptotic results, for finite sample performance, using the Empirical Bayes variants and SURE can lead to null rejection probabilities that are closer to the nominal size. In addition, taking square of condition (i) on both sides, we see that \[\frac{1}{N} \sum_{i,j\in I_l} (\tilde{\alpha}_{il}^*-\alpha_i)^2 (\tilde{\alpha}_{jl}^*-\alpha_i)^2=O_p(\frac{N}{T^2}),\] which implies \[\frac{1}{N} \sum_{i\in I_l} (\tilde{\alpha}_{il}^*-\alpha_i)^4=O_p(\frac{N}{T^2})=o_p(1).\] When $T \propto N^g$ for $g>1/2$, then $O_p(\frac{N}{T^2})=o_p(1)$.

Condition (ii) is needed to ensure that $\frac{1}{\sqrt{N}} \sum_{i\in I_l} (\frac{\partial m(W_i,\alpha_i,\mu_0) }{\partial \alpha_i}- a(X_i,\alpha_i,\mu_0))(\tilde{\alpha}_{il}^*-\alpha_i)=o_p(1)$. Even if $\tilde{\alpha}_{il}$ is constructed using data on individual $i$, the condition makes the correlation of the estimation error with the mean zero term $\frac{\partial m(W_i,\alpha_i,\mu_0) }{\partial \alpha_i}- a(X_i,\alpha_i,\mu_0)$ negligible. The moment conditions are not restrictive as it involves the time averages $\bar{\bar{u}}_i$. Typically, under standard assumptions on weakly dependent data $\sqrt{T} \bar{\bar{u}}_i=O_p(1)$, hence one would expect $T \mathbb{E}((\frac{\partial m(W_i,\alpha_i,\mu_0) }{\partial \alpha_i}- a(X_i,\alpha_i,\mu_0))^2\bar{\bar{u}}_i^2)<\infty$, justifying the first part of condition (ii). The second part is similar to condition (ii) in assumption 1 of battaglia2024. A sufficient condition is $\frac{\partial m(W_i,\alpha_i,\mu_0) }{\partial \alpha_i}- a(X_i,\alpha_i,\mu_0)$ being mean independent of the time varying idiosyncratic shocks $u_{it}$, which guarantees that $\mathbb{E}((\frac{\partial m(W_i,\alpha_i,\mu_0) }{\partial \alpha_i}- a(X_i,\alpha_i,\mu_0))\bar{\bar{u}}_i)=0$. This is the case when ((ref)) is a dynamic panel data model ($X_{it}$ includes $Y_{i,t-1}$): for any $t<T$,

where the second last line follows because $X_i$ includes $Y_{it}$ and $X_{it}$, so $u_{it}$ for $t<T$ is a deterministic function of $X_i$ and $\alpha_i$. If one constructs the estimator $\tilde{\alpha}_{il}$ using only data up to period $T-1$, then the averages will be orthogonal to the deviation $\frac{\partial m(W_i,\alpha_i,\mu_0) }{\partial \alpha_i}- a(X_i,\alpha_i,\mu_0)$. When ((ref)) does not include a dynamic specification but the researcher is willing to assume that $u_{it}$'s are mean independent across $t$, one can simply add $(Y_{i1},\cdots,Y_{it-1})$ into $X_i$. We do this in the empirical application in Section (ref). A similar law of iterated expectation argument can show that condition (ii) still holds.

The last two conditions are only needed when one uses data in fold $l$ to perform shrinkage of $\tilde{\alpha}_{il}$. Condition (iii) requires that there exists non-trivial individual heterogeneity in $\alpha_i$'s, the assumption is needed to ensure that in constructing the EB or SURE type estimators the shrinkage factor $\frac{u_{i}^2}{u_{i}^2+ \sigma_{\alpha}^2}$ will not involve division by a number close to zero. Combined with condition (iii), condition (iv) is a mild moment existence condition that guarantees that $\sup_{i\in I_l}\frac{\hat{u}_{i}^2}{\hat{u}_{i}^2+ \hat{\sigma}_{\alpha}^2} =o_p(N^{-1/2})$. For instance, under the assumption that $u_{it}$ is iid across both time and individuals, an estimator for $u_{i}^2=\mbox{Var}(\frac{1}{T-1} \sum_{t=1}^{T-1} u_{it})$ is $\frac{1}{(T-1)^2} \sum_{t=1}^{T-1} (\hat{u}_{it})^2$ for $\hat{u}_{it}=Y_{it}-X_{it}'\tilde{\beta}_l-\tilde{\alpha}_{il}$. Condition (iv) then requires that \[\mathbb{E}\{[\frac{1}{T-1}\sum_{t=1}^{T-1} (\hat{u}_{it})^2]^{2+\delta}\}<\infty ,\] which holds under general moment existence conditions on $(X_{it},u_{it},\alpha_i)$. A similar argument can be constructed for the case where $u_{it}$ might be serially correlated across time. To control $\sup_{i\in I_l}\hat{u}_{i}^2$, there exists a trade-off between the existence of higher moments and the growth rate of the time dimension $T$ with $N$, as manifested in the condition $g>\frac{1}{2}+\frac{1}{2+\delta}$. In general, a smaller growth rate implies that we are taking supremum over a larger group of individuals for the same level of randomness due to averaging across time. Hence higher-order moments are needed to ensure that the quantity is $o_p(N^{-1/2})$.

In addition to the rate conditions imposed above, we also collect the needed conditions on existence of moments:

Condition (i) imposes a bound on the second-order derivative of $m(W_i,\alpha_i,\mu_0)$ with respect to $\alpha_i$. In Example (ref), the moment function is quadratic in $\alpha_i$, and it will be satisfied automatically. In addition, $\mathbb{E}(u_{iT}^2|X_i,\bar{\bar{Y}}_i,\alpha_i)$ is a restriction on the conditional variance of $u_{iT}$. Conditions (ii)-(iv) guarantee the existence of moments so that suitable Law of Large Numbers can be applied.

Assumption (ref) lists additional moment conditions that are needed when $\tilde{\alpha}_{il}^*$ uses the EB or SURE variant.

Proposition (ref) is the key result that establishes that under the orthogonal moment construction and cross-fitting, the asymptotic distribution is the same as if the true $\alpha_i$ and the true function $a(X_i,\alpha_i,\mu_0)$ are known at $\mu=\mu_0$. The fact that both $\alpha_i$ and the conditional mean function $a(X_i,\alpha_i,\mu_0)$ are estimated will not affect the asymptotic distribution of the estimator $\hat{\mu}$ other than accounting for the addition of the orthogonal correction term $a(X_i,\alpha_i,\mu_0)(Y_{iT}-X_{iT}'\beta_0-\alpha_i)$. For $m^*(W_i,\alpha_i,Z_i,\beta_0,\mu_0,a)= m(W_i,\alpha_i,\mu_0)+ a(X_i,\alpha_i,\mu_0)(Y_{iT}-X_{iT}'\beta_0-\alpha_i)$, the second line of ((ref)) is a normalized sum of iid mean zero random variables, under the condition that $\Omega\equiv \mathbb{E}(m^*(W_i,\alpha_i,Z_i,\beta_0,\mu_0,a) m^*(W_i,\alpha_i,Z_i,\beta_0,\mu_0,a)')$ (with $\mathbb{E}(\|m^*(W_i,\alpha_i,Z_i,\beta_0,\mu_0,a)\|^2)<\infty$), standard Central Limit Theorem will imply \[\frac{1}{\sqrt{N}} \sum_{i=1}^N m^*(W_i,\alpha_i,Z_i,\beta_0,\mu_0,a) \xrightarrow{d} \mathcal{N}(0,\Omega)\]

The following condition ensures that $\Omega$ can be consistently estimated by $\hat{\Upsilon}^{-1}$ in ((ref)).

Condition (i) is not restrictive given $\frac{1}{N} \sum_{i\in I_l} (\alpha_i-\tilde{\alpha}_{il})^4=o_p(1)$ and condition (ii) in Assumption (ref). Conditions (ii) and (v) are suitable regularity conditions in order to apply Law of Large Numbers. Condition (iii) assumes consistency of the preliminary estimator of $\tilde{\mu}_l$, which is established in Appendix (ref). Condition (iv) assumes that the function $m(W_i,\alpha_i,\mu)$ is sufficiently smooth in both $\alpha_i$ and $\mu$, which is satisfied in Example (ref) with $m(W_i,\alpha_i,\mu)$ being quadratic in $\alpha_i$ and linear in $\mu$.

To establish the validity of ((ref)) for estimating the asymptotic variance of $\hat{\mu}$, we still need to show that $\hat{G}$ consistently estimates $G=\mathbb{E}(\frac{\partial m(W_i,\alpha_i,\mu_0)}{\partial \mu})$. The following conditions ensure the convergence of the estimator for the Jacobian matrix $G$, and is similar to assumption 5 of cher2022.

Assumption (ref) is easy to verify in Example (ref). With $m(W_i,\alpha_i,\mu)$ being quadratic in $\alpha_i$ and linear in $\mu$, we can take $F(W_i)=0$ and also $F_{\alpha}(W_i,\alpha_i)= 0$ with $C=1$.

The next theorem combines the above propositions to establish the validity of inference for $\hat{\mu}$:

Theorem (ref) requires consistency of $\hat{\mu}$ for $\mu_0$, which follows under standard assumptions. The conditions required are specified in Appendix (ref).

Monte-Carlo Simulation

We will compare the estimator and inference procedure in (ref) with those in xie2025 and CGK2025 through the lens of example (ref). More specifically, we will consider a panel auto-regression model as the first stage:

The fixed effects $\alpha_i$'s are independent $\mathcal{N}(0,\frac{1}{2})$ draws, and $u_{it}$'s (unobserved) are also $\mathcal{N}(0,\frac{1}{2})$ constructed from the sum of two independent random variables $u_{1it}+u_{2it}$ each distributed as $\mathcal{N}(0,\frac{1}{4})$. Both $u_{1it}$'s and $u_{2it}$'s are also independent across time and individuals. The initial condition $Y_{i0}$ is generated from the stationary distribution as $Y_{i0}\sim \mathcal{N}(\frac{\alpha_i}{1-\beta_0},\frac{1}{2(1-\beta_0^2)})$.

The second stage is specified as in ((ref)) with a slight modification:

where $\lfloor x \rfloor$ denotes the largest integer that is smaller than or equal to $x$. In the notation of example (ref), $\mu_{01}=0$ and $\mu_{02}=1$, and the parameter of interest is $\mu_{02}$. We draw $v_i$ from $\mathcal{N}(0,1)$ independent from all other variables and $c$ is a scaling factor that gauges the level of correlation between $\bar{\bar{u}}_i$ and $v_i-c\frac{1}{ \lfloor \frac{T}{2} \rfloor -1 } \sum_{t=1}^{ \lfloor \frac{T}{2} \rfloor -1 } u_{1it}$. For $c\neq 0$, then $\bar{\bar{u}}_i$ and $v_i-c \frac{1}{ \lfloor \frac{T}{2} \rfloor -1 } \sum_{t=1}^{ \lfloor \frac{T}{2} \rfloor -1 } u_{1it}$ are correlated, which violates the conditions for inference in both xie2025 and CGK2025. In contrast, as explained in Section (ref), we expect the orthogonal moment condition to be insensitive to the potential correlation between them and have null rejection probability close to the nominal size.

As xie2025 and CGK2025 abstracts away from the problem of estimating $\beta_0$, we employ estimators from BB by constructing the following moment conditions: \[\mathbb{E}

=0 \quad for t=4,\cdots,T. \] $\beta_0$ is then estimated via two-step GMM.\footnote{I only included moments $t>3$ for the simulation results with $T=12$ and $t>5$ for simulation results with $T=22$ in order to avoid multicollinearity in constructing the optimal weighting matrix. The cutoffs are chosen based on finite sample performance in the simulations.}

To implement procedures in xie2025, I consider both $\hat{\alpha}_i=\frac{1}{T} \sum_{t=1}^{T} Y_{it}- \frac{1}{T} \sum_{t=0}^{T-1} Y_{it} \tilde{\beta}$ and also its EB ($\alpha_i^{EB}$) and SURE ($\alpha_i^{SURE}$) variants. xie2025 estimates the standard error of $\mu_{02}$ with the usual heteroskedastic robust standard error formula by plugging the EB versions of $\alpha_i^{EB}$.

In addition, I also consider the estimator recommended in CGK2025 which builds on the observation that \[\mu_{02} = \frac{Cov(\alpha_i,W_i)}{\mbox{Var}(\alpha_i)}=\frac{Cov(\hat{\alpha}_i,W_i)}{ \mbox{Var}(\hat{\alpha}_i)} \frac{\mbox{Var}(\hat{\alpha}_i) }{\mbox{Var}(\alpha_i)},\] exploiting independence between $\bar{u}_i$ and $W_i$ in their assumptions (which does not hold when $c\neq 0$). They suggest replacing each population quantity by its sample analogue. In particular, $\mbox{Var}(\alpha_i)$ can be estimated by $\hat{\mbox{Var}}(\hat{\alpha}_i)-\hat{\mbox{Var}}(\bar{u}_i)$, where $\hat{\mbox{Var}}(\bar{u}_i)$ is a consistent estimator of the variance of $\bar{u}_i$ and can be constructed as $\frac{1}{N} \sum_{i=1}^N \frac{1}{T^2} \sum_{t=1}^T \hat{u}_{it}^2 $. Inference can be conducted using percentile bootstrap.

For the orthogonal moment constructions, I followed the procedures in Section (ref) and uses $\tilde{\alpha}_{il}=\frac{1}{T-1} \sum_{t=1}^{T-1} Y_{it}-\frac{1}{T-1} \sum_{t=0}^{T-2} Y_{it}$ so that condition (ii) of assumption (ref) is guaranteed to hold, and also for its EB and SURE variants $\tilde{\alpha}_{il}^{EB}$ and $\tilde{\alpha}_{il}^{SURE}$. A preliminary estimator of $\tilde{\mu}_l=(\tilde{\mu}_{1l},\tilde{\mu}_{2l})$ is then obtained via linear regression. The estimator of conditional mean function $a_1(X_i,\alpha_i,\mu_0)=-\mu_{02}$ for $i\in I_l$ is constructed by $-\tilde{\mu}_{2l}$ and the conditional mean function $a_2(X_i,\alpha_i,\mu_0) = \mathbb{E}(v_i - \mu_{02} \alpha_i|X_i,\alpha_i)$ is constructed via an adaptive elastic net regression of $\hat{v}_j-\tilde{\mu}_{2l} \tilde{\alpha}_{jl} $ on $\tilde{\alpha}_{jl} $ and $(Y_{j0},Y_{j1},\cdots, Y_{jT-1})$ for $j\notin I_l$. Here $\hat{v}_j=W_j-\tilde{\mu}_{1l}-\tilde{\mu}_{2l} \tilde{\alpha}_{jl}$. The penalty factors for the L1 regularization and L2 regularization terms in the Elastic Net regression are chosen using cross-validation. To alleviate potential noise introduced in the sample splitting process, we computed $\hat{\mu}$ and $\hat{V}$ using orthogonal moments for 20 different realizations of sample splits (randomness coming from reshuffling of the sample along the cross-sectional dimension) and then averaged $\hat{\mu}$ and $\hat{V}$ across them.

Table (ref) reports the simulation results for $N=100$, $T=12$ and $\beta_0=0$ and three different choices of $c$ across 1,000 simulations. We also compare performance of the inference procedure with different number of folds $L$.\footnote{As explained in BB, larger $\beta_0$ will make the estimates of $\beta$ very imprecise, resembling a weak instrument problem. Our asymptotic theory essentially ignores any noise in estimating $\beta_0$, and poor estimates of $\beta_0$ will propagate and affect the inference results for $\beta_0$ that is too large.}

The first column reports the results for estimators suggested by CGK2025 and the next three columns report results where heteroskedastic standard error for linear regression is used to compute asymptotic variance with different estimators of $\alpha_i$. Finally, the last six columns employs the orthogonal moment construction with different estimators of $\alpha_i$ and with different number of folds. Comparing the bias of the estimators across different panels, the results show that when the endogeneity increases ($c$ is large), there is significant bias in the estimates using CGK2025 and xie2025. Although these estimators typically has smaller standard error, it does not compensate for the bias. Therefore, estimators using orthogonal moments enjoy smaller Root Mean Squared Error (RMSE) for larger values of $c$.

We also consider testing the null hypothesis of $\mu_{02}=1$ using t-test based on different estimators and their estimated standard errors. The null rejection probability is reported in the last row of each panel with target size 0.05. We see that the endogneity between the measurement error and the residuals in ((ref)) made the estimators of CGK2025 and xie2025 suffer from severe size distortion. In contrast, estimators that uses orthogonal moments have size control within acceptable levels. In addition, using the EB or SURE correction delivers better size control than the original $\tilde{\alpha}_{il}$. Nevertheless, the shrinkage correction inflates the RMSE compared with the case where no correction is applied. The simulation results suggest that in general, using orthogonal moment conditions is more important than using EB or SURE corrections to mitigate the concerns of measurement error especially when endogeneity is a concern. In addition, comparing results from column 6-8 with the results from column 9-11, we see that the number of folds $L$ does not have large impacts on either the estimates or the null rejection probabilities. One can choose a smaller fold $L=5$ for computational reasons in the current DGP.

Figure (ref) plots the power function of the CGK estimator, the estimator of xie2025 using SURE and estimator using orthogonal moment conditions and SURE corrections for the DGP in panel A where all estimators are expected to have null rejection probability close to 0.05. We see that a finite-sample distortion is present for the estimator proposed in xie2025, while the estimator from CGK2025 and estimator based on orthogonal moments have smaller null rejection probabilities. In addition, the estimator proposed in CGK2025 has higher power in rejecting the alternatives compared with the estimator using orthogonal moments. The phenomenon reflects the fact that estimators from CGK2025 achieves the semi-parametric efficiency bound when endogeneity concern is not present, and in this case their estimator should be preferred for power purposes.

Table (ref) reports the results with $N=500$, $T=22$ and $\beta_0=0.5$ with different levels of endogeneity represented by $c$. As adaptive elastic net conducts variable selection, the values of $c$ is set to be 5 and 6 respectively to avoid potential local-to-zero coefficients in $a(X_i,\alpha_i,\mu_0)$. As results in Table (ref) shows that the number of folds has little impact on the estimator performance and size control, we report results only with 5-fold cross-fitting in Table (ref). Table (ref) shows patterns consistent with Table (ref). As $c$ increases, the bias of estimators from both CGK2025 and xie2025 increases and the tests over-reject the null. In contrast, estimators based on orthogonal moment conditions have better size control and smaller RMSE. Comparing column 6 with column 7 and column 8, we see that using EB or SURE shrinkage still delivers null rejection probabilities closer to the nominal size, albeit with a larger RMSE.

Empirical Application

We will consider an empirical application that studies the site selection of a policy experiment on agriculture catastrophe insurance in China. policy built a comprehensive database for past policy experiments conducted by Chinese central government, and they documented consistent patterns of positive selection of experiment sites. Specifically, the $t$-statistics that compare the average fiscal revenue between experimental sites and non-experimental ones have a significant proportion that exceeds the usual 95% quantile of the $t$-distribution. The result is robust when they compare agriculture outputs between sites for agricultural related experiments. policy also listed potential reasons for positive selection of the experimental sites, but did not commit to a particular explanation of the phenomenon. However, understanding the way that experimental sites are positively selected can be important in interpreting the experimental results and determining whether the new policy should be rolled out to the entire country.

Policy experiments are typically implemented when there is substantial uncertainty over the effects of the proposed policy, and as a means to encourage local government to find optimal policy instruments for promoting the central government's policy objective heilmann1.\footnote{For background and case studies on China's policy experimentation before 2008, see heilmann1,heilmann2.} The central government usually sets the policy agenda, but leaves the particular implementation details to the experimental sites. The successful experimental sites will be designated as “model sites”, and the less successful ones are encouraged to learn from the experience of the “model sites”. In such context, positively selecting experimental sites that are better at adjusting policy details and organizing policy reforms can potentially yield valuable experience for non-experimental sites. Therefore, positive selection under this form can be desirable from the view of the central government.

As an illustration, we focus on a particular policy experiment -- the promotion of agriculture catastrophe insurance (ACI) -- that has taken place in 2017. Prior to the experiment, the insured liability of the agriculture insurance typically only covered direct input costs like fertilizer and seeds, which made the farmers vulnerable to natural disasters. In response to the deficiency of the insurance products, the central government would like to develop insurance that has higher coverage to help farmers guard against risks and encourage uptakes of more advanced agricultural production techniques. One major difference under the new ACI is that it covers land costs in addition to direct input costs. Under the new ACI, with the same premium, the insured liability increases from 450 CNY per acre to 900 CNY per acre in Shandong Province for instance insurance. To facilitate adoption of ACI, the premiums that farmers pay were highly subsidized by both the central and local government. The focus of the ACI experiment in 2017 was on moderately scaled farmers who were more susceptible to climate risks and also on strategic agriculture products that consist of rice, corn and wheat. The experiment was conducted at county-level, with in total 200 counties participating in the experiment. However, instead of directly choosing counties, the central government hand-picked 13 major grain production provinces and delegated the choice of specific counties to each provincial government with quotas. The guidelines issued by the central government required the experimental counties to have well-established insurance foundations, a sufficient number of moderately scaled farmers and also to be representative in the province.\footnote{See \url{https://jcs.moa.gov.cn/trzgl/201705/t20170527_5626557.htm} } In addition, policy details like the criteria of moderately scaled farmers, the formula for evaluating land costs and means of introducing the insurance to farmers were delegated to local governments. This suggests that selecting a (random) representative sample was not the sole objective of experimental site selection. Furthermore, the central government would also want to learn the best practice of implementation by enrolling counties that were more specialized in agriculture insurance and agriculture production. We would like to test whether the provincial governments were indeed choosing counties that were more specialized in agriculture production as experimental sites.

Data and measurements

The data comes from the replication package of policy. In addition, the selected counties that participated in the ACI experiment are hand-collected from the provincial governments' websites with the exception of the Inner Mongolia province whose information is not directly available. We restricted the sample to counties that are in the remaining 12 major grain producing provinces and dropped counties whose administrative divisions have changed. The total sample consists of 961 counties, among which 168 of them were selected as experiment sites. For each county, we observe its GDP, population, fiscal income, fiscal expenditure, total grain production, agricultural mechanical power, total rural employment and the total number of policy experiment that the county has from 1997-2016. The variable total rural employment is only observed up until 2012. As a preliminary comparison, Table (ref) reports the means of the above measures in 2016 for both experimental and non-experimental counties. The value of rural employment is computed using 2012 data. In addition, we also compared the GDP and fiscal revenue growth rates over the period 2013-2016.

Compared with the non-experimental counties, the experimental ones tend to be larger in terms of population, GDP and fiscal expense. In addition, they also produced more grains and used more agricultural machinery. This serves as preliminary evidence that when choosing counties, the provincial governments not only aimed at obtaining a representative sample, but also wanted to enroll counties that were more developed and counties that had a larger share of agricultural industry. These counties can potentially lay out a better implementation plan of ACI, and thus establish “model practices” for other non-experimental ones. We now turn to more formal econometric analysis of the problem.

Econometric Specification

As the outcome is binary, we will specify a simple logit model:

where GDP per Capita and number of policy experiments are measured in 2016, and GDP_growth and FR_growth are GDP and fiscal revenue growth rate between 2013 and 2016 for county $i$. $\delta_p$ is province fixed effect that captures heterogeneity across provinces, and $v_i$'s are iid logit error. We also consider a version that omits the province fixed effects. in which case it will be replaced with an intercept. The key variable is $\alpha_i$ which measures the degree of comparative advantage of county $i$ in grain production. If the provincial governments were selecting experimental counties that were more specialized in grain production, then we should expect the coefficient $\mu_{01}$ to be positive. However, the degree of specialization is not directly observable to the researcher, and we propose a panel data model that can help us recover $\alpha_i$'s.

where we regress grain output per rural employment for each county $i$ and year $t$ on that county-year pair's GDP per capita, mechanical power per rural employment and fiscal expenditure per capita, along with county fixed effects $\alpha_i$ and province by time fixed effects $\delta_{pt}$. A larger value of $\alpha_i$ indicates that county $i$ is more efficient in grain production compared with its fellow counties that have similar fiscal support, income level and amount of agricultural machinery. Therefore, $\alpha_i$ captures the level of specialization of grain production for different counties.

Although the presence of $\delta_{pt}$ may cause some difficulties for directly implementing the procedure in Section (ref). In practice, we first transform ((ref)) by demeaning the province by time average of each variable to eliminate the province by time fixed effects.

where $\dot{A_{it}}=A_{it}-\bar{A}_{p(i)t}$ is the demeaned version of $A_{it}$ with $\bar{A}_{p(i)t}$ being the mean of $A_{it}$ across province $p$ to which $i$ belongs and time $t$. The $X_{it}$ in ((ref)) would correspond to $(\dot{\text{GDP\_per\_cap}}_{it},\dot{\text{Mech\_per\_emp}}_{it},\dot{\text{FE\_per\_cap}}_{it})$. Applying procedure in Section (ref) then identifies $\dot{\alpha}_i$. As $\alpha_i$ is time invariant, $\dot{\alpha}_i=\alpha_i-\bar{\alpha}_{p(i)}$. Note that $\dot{\alpha}_i$ and $\alpha_i$ differ by a constant that only varies by province. Therefore, rewriting the logit model in ((ref)), we have

Therefore, $\mu_0$ can still be identified from the demeaned version of $\alpha_i$ as long as the province fixed effects are included in the cross-sectional logit model.

To fit into the framework of ((ref)), we will work with the moment conditions that are associated with the score of the likelihood function define by ((ref)). With some abuse of notation define $W_i =(\text{GDP\_per\_cap}_i, \text{GDP\_growth}_i,\text{FR\_growth}_i,\text{N\_policy}_i,\delta_p)'$ where $\delta_p$'s are understood as 12 dummy variables. Let $\mu_{0,-}=(\mu_{02},\mu_{03},\mu_{04},\mu_{05},\mu_{0,p})'$ where $\mu_{0,p}$ are fixed effects associated with the 12 provinces. Then we are solving the estimating equations

where $\Lambda(x)=\frac{\exp(x)}{1+\exp(x)}$ is the logistic function. With the moment condition \\ $\mathbb{E}[m(W_i,\text{Experiment}_i,\alpha_i,\mu_0)]=0$ defined as in ((ref)), the corresponding $a(X_i,\alpha_i,\mu_0)$ for $X_i=\{X_{it}\}_{t=1997}^{2012} \cup \{\text{grain\_per\_emp}_{it}\}_{t=1997}^{2011} $ is

where $\Lambda_i= \Lambda(W_i'\mu_{0,-}+\alpha_i \mu_{0,1})$. In practice, we estimate $a(X_i,\alpha_i,\mu_0)$ by performing adaptive elastic net regression on $X_i$ and $\alpha_i$.

Estimation of the Jacobian matrix of $\mu_0$ is similar to the classical maximum likelihood theory:

As in the simulation of Section (ref), the estimates of $\mu_0$ and the asymptotic variance ((ref)) are averaged across 20 different sample splits with $L=5$.

Results

Before presenting estimates of $\mu_0$, it is useful to show the results of the panel regression ((ref)) first.

We see that the coefficients on both GDP per capita and fiscal expenditure per capita change significantly when county fixed effects are included, while inclusion of province-by-time fixed effects have more limited impacts on the magnitude and sign of all coefficients. The mechanical power per rural employment turns out to be a salient predictor of grain production per rural employment, while the coefficients on GDP per capita and fiscal expenditure per capita are not significantly different from zero. Now we turn to results on $\mu_0$.

Panel A of Table (ref) report the results of estimating ((ref)) without province fixed effects, while panel B report the ones with province fixed effects. Comparing results across panels, the coefficients on GDP growth rate and fiscal revenue growth rate switch signs. Even though the coefficients associated with the two variables are not statistically significant, it demonstrates the importance of including province fixed effects in ((ref)) when ((ref)) is estimated with province-by-time fixed effects. Now we turn to the comparison between the plug-in estimates and the estimates based on orthogonal moments. The coefficients are similar in sign and magnitude when we compare the plug-in estimates with estimates based on orthogonal moments except the coefficient $\mu_{01}$ on $\alpha_i$. In both panels, $\mu_{01}$ is estimated to be larger (around 0.6) when we use the simple plug-in methods, while the estimates based on orthogonal moments give estimates that range from 0.36 to 0.47. Therefore, correcting errors in estimating $\alpha_i$ could be important in certain applications. Nevertheless, in our setting, the use of orthogonal moments still gives consistent evidence that there is indeed positive selection of experimental counties in terms of their comparative advantage in grain production. The provincial governments were leaning towards more specialized counties to learn from their potentially better implementation practice and make the ongoing experiment more successful.

Conclusion

This paper employs debiased machine learning techniques to study inference of parameters characterized by cross-sectional moments involving latent variables that are fixed effects in an auxiliary panel data regression. Under the asymptotic regime that both the cross-sectional dimension and time dimension go to infinity, the asymptotic variance of the estimator based on orthogonal moments is the same as if the true latent variables are known to the researcher. In addition, shrinkage methods based on Empirical Bayes or Stein's Unbiased Risk Estimation can be readily incorporated for estimating the latent variables. The approach enables researchers to deal with non-linearity in the cross-sectional moments and relaxes independence assumption between the panel data residuals and the cross-sectional moment functions. Simulation results show that using orthogonal moments have better null rejection probabilities compared with existing methods when the panel data residuals are correlated with the cross-sectional moment functions. In an empirical application, we examine the experimental site selection of the agricultural catastrophe insurance experiment in China. Interpreting the score function of a logit model as the targeted moments, we show that specialization in grain production is an important determinant on whether a county is selected into the policy experiment where the level of specialization is captured as county fixed effects in an auxiliary linear panel data model. The estimates using orthogonal moments are smaller in magnitude compared with the simple plug-in approach, while SURE and EB corrections deliver similar results. In general, using orthogonal moments can be important for valid inference especially when researchers are in doubt that the residuals in the panel data model are uncorrelated with the cross-sectional moment functions.