Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
100,442 characters · 23 sections · 62 citation commands
Efficient Estimation of Average Derivatives in NPIV Models: Simulation Comparisons of Neural Network Estimators
\baselineskip=18pt
Deep layer Artificial Neural Networks (ANNs) are increasingly popular in machine learning (ML), statistics, business, finance, and other fields. The universal approximation property of a variety of ANN architectures has been established by hornik1989multilayer and many others. Early on, computational difficulties have hindered the wide applicability of ANNs. Recently, improvements in computing have led to successful applications of deep layer ANNs in computer vision, natural language processing and other areas, with complex nonlinear relations among many covariates and large data sets of high quality.\footnote{By high quality we mean data sets with very high signal-to-noise ratios. Unfortunately, many economic and social science data sets have low signal-to-noise ratios.} Many problems where deep layer ANNs are extremely effective involve prediction problems (i.e. estimating conditional means or densities)---or problems in which nuisance parameters are themselves predictions. Recently, farrell2018deep and athey2019using, among others, have applied multi-layer ReLU ANNs to estimate average treatment effects under unconfoundedness and demonstrated their good performance in estimating unknown conditional means and densities of multivariate covariates. It remains to be seen whether ANNs are similarly effective for structural estimation problems with nonparametric endogeneity.
To that end, we consider semiparametric efficient estimation and inference for a weighted average (partial) derivative (WAD) of a nonparametric instrumental variables regression (NPIV) via ANN sieves. Specifically, we assume an unknown structure function $h$ satisfies the NPIV model: $\E[Y_{1} - h(Y_2) \mid X] = 0$, where $Y_2$ is a continuous random vector of moderately high dimension (including endogenous regressors that are excluded from $X$), and $X$ is a vector of moderately high dimensional conditioning variables. We are interested in efficient estimation and inference for a WAD parameter of the smooth NPIV function $h (Y_2)$, without sparsity assumptions on $h(\cdot)$.\footnote{Of course, in lieu of sparsity, we do require smoothness assumptions.} WADs of structural relationships are linked to elasticities of endogenous demand systems in economics. It is essentially a treatment effect parameter under confounding and endogenous continuous treatment. Although there is a large literature on efficient estimation of the average treatment effect and other causal parameters under unconfoundedness, there are far fewer results on efficient estimation and inference on the average treatment effect in nonparametric models with endogenous continuous treatment.
This paper makes three contributions. First, we present two classes of efficient estimators for WADs of NPIV models where unknown $h_0(Y_2)$ is approximated by ANN sieves: the optimally weighted sieve minimum distance estimators and the efficient score-based estimators. Under some regularity conditions both types of estimators are root-$n$ asymptotically normal, semiparametrically efficient, and hence are first-order equivalent. Second, we detail a {\it practitioner's recipe} that include a step by step guide for implementing these two classes of estimators. Third, and perhaps most importantly, we present a large set of Monte Carlo results on finite-sample performances of various ANN estimators. These are implemented using increasingly complex designs, such as NPIV function containing up to 13 continuous covariates (including endogenous regressors), various nonlinearities and correlations among the covariates.
We now briefly introduce the two classes of efficient estimation procedures that we consider. Both procedures are inspired by the semiparametric efficiency bound characterization in ai2012semiparametric (henceforth AC12) for the WAD of the unknown $h(Y_2)$ in a NPIV model $\E[Y_{1}-h(Y_{2}) \mid X]=0$. The first procedure is based on minimizing an optimal criterion, the optimally-weighted orthogonalized sieve minimum distance (SMD) criterion. This procedure is numerically equivalent to a semiparametric two-step procedure, where the unknown NPIV function $h (\cdot)$ is estimated via an optimally weighted SMD in the first step, and the WAD of $h(\cdot)$ is estimated using a sample analogue of an orthogonalized unconditional moment chamberlain1992jbes in the second step, with the unknown $h$ substituted by the optimally weighted SMD estimator from the first step. This will be denoted as OP-OSMD in our paper. AC12 already introduced this procedure and presented a small Monte Carlo study demonstrating its finite-sample performance using a spline SMD in the first step when the unknown $h(\cdot)$ is a function of a scalar endogenous variable $Y_2$. It is unclear how this procedure will perform when $Y_2$ could be a continuous random vector of higher dimension and when $h(\cdot)$ is approximated via a neural network.
The second procedure is based on the efficient score (equivalently, efficient influence function).\footnote{The efficient score/influence function approach to efficient estimation has a long history in semiparametrics. See, e.g., bickel1993efficient, pfanzagl1982lecture, Section 25.8 of van2000asymptotic and references therein, for an introduction.} AC12 derived a characterization of the efficient influence (or equivalently, efficient score) for the WAD of a NPIV model. It is also the asymptotic influence function of the OP-OSMD estimator.\footnote{This is not surprising since the efficient influence function is unique.} The efficient influence is the sum of the orthogonalized unconditional moment (the one used for the OP-OSMD estimator) and an adjustment term accounting for plugging-in estimated $h(Y_2)$, often referred to as the Riesz representer term. Compared to simpler settings, e.g. estimating average treatment effect under unconfoundedness, the Riesz representer term here has no closed-form expression, but is characterized as one solution to an optimization problem over an infinite-dimensional Hilbert space induced by a norm connected to the optimally weighted minimum distance objective. The components of the efficient influence function can nonetheless be consistently estimated via sieve approximations. The procedure using the sample estimated efficient influence (i.e., efficient score) will be denoted as ES in our paper. To the best of our knowledge, there is no published work on theory or simulation on the finite-sample performance of any ES estimator for the WAD in a NPIV model yet.
In this paper we investigate the finite-sample performance of both efficient procedures when the unknown function $h(Y_2)$ is estimated via various ANN SMDs and when $h (Y_2)$ depends on moderately high dimensional continuous regressors $Y_2$ (some of which are endogenous). We describe some stylized findings from our simulations. Our simulations reveal that the ANN OP-OSMD is more stable and easier to implement than ANN ES for estimation of the average (partial) derivative in a NPIV model with unknown conditional variance $\Sigma(X)\equiv \var(Y_1-h(Y_2) \mid X )$.
In practice, it could be appealing to report simpler inefficient estimators that are still consistent and $\sqrt{n}$-asymptotically normal. It is also possible that computationally simpler inefficient estimators may perform better than the efficient estimators in finite samples, as the efficient estimators often require estimating additional nuisance parameters. For the sake of comparison, we include two first-order asymptotically equivalent inefficient estimators of the WAD of a NPIV function, denoted by P-ISMD and IS. The P-ISMD is a simple plug-in identity-weighted SMD estimator that was proposed in ai2007estimation (henceforth AC07). The IS is what we call “inefficient score” estimator that is based on sample analog of the asymptotic influence function of the P-ISMD estimator (derived in AC07).\footnote{Different inefficient estimators of the WAD can have different asymptotic influence functions and hence different asymptotic variances. That is why we define the IS estimator based on the asymptotic influence function of the \textsf{P-ISMD } estimator of AC07, so that they will have the same asymptotic variance.} We note that both \textsf{P-ISMD } and \textsf{IS } are asymptotically efficient for a WAD of a nonparametric regression $\E[Y_1 \mid Y_2]$, in the absence of endogeneity. However, they are no longer efficient for the WAD of a NPIV function $h (Y_2)$ identified by the conditional moment restriction $\E[Y_{1}-h(Y_{2}) \mid X]= 0$ (for $Y_2\neq X$).
We compare the finite sample performance of these efficient (OP-OSMD, ES) and inefficient (P-ISMD, IS) estimation procedures in four Monte Carlo designs with moderate sample sizes ($n=1000$ to $n=10000$).\footnote{Since both score-based estimators ES and IS are based on orthogonal moments, we also provide comparison with their cross-fitted versions. The cross-fitting orthogonal moments estimators have become very popular following \citep*{chernozhukov2018double,chernozhukov2021locally} and others, although no published work has applied cross-fit to efficient estimation of WAD in NPIV yet.} In (ref), we estimate a simple nonparametric regression and two-stage least squares data-generating process, as a useful baseline. In (ref), we estimate the average partial derivative of a NPIV function $h(Y_2)$ with respect to an endogenous variables using various ANN sieves and spline sieves. In (ref), we calibrate a data-generating process to the gasoline empirical application blundell2012measuring, and repeat the exercises for (ref).
Our Monte Carlo experiments allow for comparisons along several dimensions:
The main, stylized takeaways from our Monte Carlo experiments are as follows:
Lastly, as applications to real data, we apply ANN sieve NPIV to estimate average price elasticity of a gasoline demand using the data set of blundell2012measuring, and to estimate average derivatives of a price-quantity relation in differentiated product markets using the data set of compiani2019market. Both applications involve nonparametric structure functions of multi-dimensional covariates (including endogenous price), and our ANN applications do not impose any semiparametric shape restrictions.
\paragraph{Related literature on ANN NPIVs.} We view various ANNs as examples of nonlinear sieves, which, compared to linear sieves, can have faster approximation error rates for large classes of nonlinear functions of high dimensional regressors. Once after the approximation error rate of a specific ANN sieve is established for a class of unknown functions,\footnote{Different ANN sieves have different approximation error rates for different function classes. See, for example, barron1993universal and chen1999improved for approximation errors rates for single hidden layer ANNs for Barron class; yarotsky2017error, shen2021optimal, shen2021neural for approximation error rates of multi layer ReLU ANNs for typical smooth function class; schmidt2019deep for approximation error rates of deep layer ReLU ANNs for composition function classes.} the asymptotic properties of estimation and inference based on the ANN sieve could be established by applying the general theory of sieve-based methods. The nonparametric convergence rates in ai2003efficient,ai2007estimation (henceforth, AC03, AC07) explicitly allow for nonlinear sieves such as ANNs to approximate and estimate the unknown structure functions of endogenous variables. They establish the root-$n$ asymptotic normality of regular functionals of nonparametric conditional moment restrictions with smooth residual functions. Due to the small sample size and computational limitation, earlier applications in econometrics have focused on single-hidden layer ANNs. For instance, chen2009land applied single hidden layer sigmoid ANN SMD to estimate the unknown habit function in a semi-nonparametric asset pricing conditional moment model with a time series sample size of about 200 quarterly observations. To the best of our knowledge, hartford2017deep is the first paper to apply multi-layer (2 hidden layer) ANNs to estimate an NPIV structural function. Since then, numerous studies have followed up, see dikkala2020minimax and the references there in.\footnote{However, as documented in our simulation results, the WAD parameter estimated via plugging in the estimated $h (\cdot)$ via adversarial GMM dikkala2020minimax can be biased. This is not surprising since the tuning parameter choice for nonparametric estimation of $h(\cdot)$ may be different from that for the efficient estimation of the WAD.} In a project that started after our first draft, chen2020nn2 established rate of convergence for multi-layer ANN optimally weighted SMD estimation of general nonparametric conditional moment restrictions for time series data, and proposed ANN sieve quasi-likelihood ratio inference for possibly slower-than-root-$n$ estimable linear functionals. However, they do not consider efficient estimation for root-$n$ estimable linear functionals of NPIV such as the WAD parameter, which is our parameter of interest.
Our simulation studies and empirical applications indicate that, although multi-layer ANNs can perform well after careful choice of tuning parameters, they have no clear advantage over single hidden layer ANNs or spline sieves for efficient estimation of WAD in a NPIV model when the unknown structure function $h(Y_2)$ is a relatively smooth function of multi-dimensional $Y_2$, which is likely the case in economic endogenous demand estimation. Just like the simulation paper by lee1993testing about the performance of single-hidden layer ANNs on testing nonlinear regression models, our paper documents that ANNs can also be one promising tool in efficient estimation and inference for causal
The rest of the paper is organized as follows. (ref) introduces the model, and the two classes of efficient estimation procedures. (ref) provides implementation details for all the estimators considered in the Monte Carlo studies. (ref) contains three simulation studies and detailed Monte Carlo comparisons of various ANN and spline based estimators. (ref) presents two empirical illustrations and (ref) concludes.
We first present the model and recall the semiparametric efficiency bound characterization. We then present two classes of efficient estimation procedures.
We are interested in semiparametrically efficient estimation of the average partial derivative:
where $a(\cdot)$ is a known positive weight function, $\nabla_1$ is the partial derivative w.r.t. the first argument and the unknown real-valued function $h_0 \in \mathcal{H}$ is identified via a conditional moment restriction\footnote{See, e.g., NP2003instrumental, BCK2007, andrews2017examples for identification of a NPIV model.}
Previously, ai2007estimation (AC07) presented a root-$n$ consistent asymptotically normally distributed identity-weighted SMD estimator of $\theta_0$, nonlinear sieves such as single hidden layer ANN sieve is allowed for in their sufficient conditions. AC12 presented the semiparametric efficiency bound of $\theta_0$ and an efficient estimator based on orthogonalized optimally weighted SMD (see their section 4.2).\footnote{ai2012semiparametric derived the efficiency bound via the \textquotedblleft orthogonalized residual\textquotedblright\ approach, which extends the earlier work of chamberlain1992jbes to allow for unknown functions entering a system of sequential moment restrictions.} ECO-019 presented efficiency bound calculation for average weighted derivatives of a NPIV model without assuming point identification of the NPIV function, but pointed out that the $\sqrt{n}$-asymptotically normal estimator of linear functionals of NPIV in santos fails to achieve the efficiency bound. chenpouzopowell proposed efficient estimation of weighted average derivatives of nonparametric quantile IV regression via penalized linear sieve GEL procedure, without providing any simulation results on how their procedure performs in finite samples.
Since weighted average treatment effects under confounding and endogenous continuous treatments can be regarded as an example of the WAD in a NPIV model, it is important to conduct some detailed Monte Carlo studies to compare finite-sample performance of various efficient estimators of $\theta_0$ when $h_0(Y_2)$ depends on multi-dimensional covariates $Y_2$. In this paper we present large scale simulation studies focusing on the performance of several estimators of $\theta_0$ when $h_0(Y_2)$ is approximated via various ANN sieves and $Y_2$ is up to $13$-dimensional vector of continuous covariates.
In this section, we specialize the general efficiency bound result of AC12 to our setting. We rewrite our model using their notation. Denote the full parameter vector as $\alpha_0 \equiv (\theta_0 ,h_0)\in \Theta \times \mathcal{H}\equiv \mathcal{A}$. The model can be written as the following sequential moment restriction
We define the orthogonalized residual as \[\varepsilon_1(Z, \alpha) \equiv \rho_1 (Z, \alpha) -\Gamma(X)\rho_2(Z,h) = a(Y_2) \nabla_1 h(Y_2) - \theta - \Gamma(X) \cdot (Y_1 - h(Y_2)),\] which is the residual from a projection of $\rho_1$ on $\rho_2$ conditional on $X$, where $ \Gamma(X)$ is the orthogonal projection coefficient: \[ \Gamma(X)\equiv \frac{\cov(\rho_1(Z,\alpha_0) \rho_2(Z,h_0) \mid X)}{\var(\rho_2 (Z, h_0) \mid X)}. \] Orthogonalizing the two moment conditions makes an efficiency analysis tractable---the same technique is used in, e.g., chamberlain1992jbes.
We now specialize the results in AC12 to the plug-in model:
where $\theta $ is a scalar and $h$ is a real-valued function of $Y_2$, and $\alpha = (\theta,h)$. Define the following variances:
We recall the efficiency bound characterization for WAD of a NPIV model from AC12 (see their Example 3.3) for the sake of easy reference, and compute
where $\overline{\mathcal{W}}=\{r:\E[\Sigma(X)^{-1}(\E\{r(Y_{2})|X\})^{2}]+\left( E\{a(Y_{2})\nabla_1 r(Y_{2})+ \Gamma(X) r(Y_2)\}\right) ^{2}<\infty \}$. Let $r_{0}\in \overline{\mathcal{W}}$ be one solution (not necessarily unique) to the optimization problem ((ref)). We note that such a solution always exists since the problem is convex, and we have:
In the rest of the paper we shall assume that $J_{0}>0$ and hence $\theta_0$ is a $\sqrt{n}$-estimable regular parameter. We note that by definition, the efficient score (indeed any moment condition proportion to an influence function) automatically satisfies the orthogonal moment condition.
From (ref), the semiparametric efficient influence function for $\theta_0$ takes the form
Denote \[ \alpha_{e} (X)\equiv (J_{0})^{-1}\frac{\E[r_{0}(Y_2)|X]}{\Sigma (X)}. \] It is clear that $\theta_0$ is the unique solution to the efficient IF equation $\E[ \psi^*(Z,\theta_0)]=0$, that is \[ \E\left[a(Y_2) \nabla_1 h_0(Y_2) - \theta - [\Gamma(X) -\alpha_{e} (X) ](Y_1 - h_0(Y_2))\right]=0 \iff \theta=\theta_0. \] One efficient estimator, $\hat{\theta}_{ES}$, for $\theta_0$ is simply based on the sample version of the efficient IF equation with plug-in consistent estimates of all the nuisance functions: \[ \hat{\theta}_{ES}=n^{-1}\sum_{i=1}^n\left( a(Y_{2i}) \nabla_1 \hat{h}(Y_{2i}) -[ \hat{\Gamma}(X_i) -\hat{\alpha}_{e}(X_i) ] (Y_{1i} - \hat{h}(Y_{2i}))\right). \] In this paper $\hat{h}(Y_2)$ can be various ANN sieve minimum distance estimators (see below), but, for simplicity, the nuisance functions $\hat{\Gamma}(X)$ and $\hat{\alpha}_{e}(X) $ are estimated by plug-in linear sieves estimators.
Another efficient estimator for $\theta_0$ can be found by optimally-weighted sieve minimum distance, where the population criterion is (see AC12):
The discrepancy measure is the optimally weighted quadratic distance of the expectation of the two moment conditions \[ m(X, \alpha) = \colvecb{2}{\E[\varepsilon_1(Z, \alpha)]}{\E[Y_1 - h(Y_2) \mid X]} = \colvecb{2}{\E[a(Y_2) \nabla_1 h (Y_2)-\theta - \Gamma(X) (Y_1 - h(Y_2))]}{\E[Y_1 - h(Y_2) \mid X]} \] from zero, where the optimal weight matrix $W_0(.)$ is diagonal and proportional to the inverse variance of each moment condition: \[ W_0(X) =
\] Two remarks are in order. First, note that the optimal weight matrix $W_0(X)$ is diagonal because $\varepsilon_1$ and $\rho_2$ are uncorrelated by design. Second, since the optimal weight matrix is diagonal and $\theta$ is a free parameter, we can view the minimization as sequential: \[ h_0 = \argmin_{h\in\mathcal{H}} \E\bk{\frac{1}{\Sigma(X)} (\E[Y_1 - h(Y_2) \mid X])^2}, \quad \theta_0 = \E[a(Y_2) \nabla_1 h_0 (Y_2) - \Gamma(X) (Y_1 - h_0(Y_2))]. \] This is important because solving the model sequentially while maintaining efficiency suggests a simple way to compute the estimators.
A sieve minimum distance estimator for $\alpha_0=(h_0, \theta_0)$ may be constructed by (i) replacing expectations with sample means, (ii) replacing conditional expectations with projection onto linear sieve bases, (iii) replacing the optimal weight matrix with a consistent estimator, and (iv) replacing the infinite dimensional optimization with finite dimensional optimization over a sieve space for $h$. This paper focuses on approximating $h$ by ANN sieves. In particular, a sample analogue of the above objective function is
where $\widehat m(.;.)$ and $\widehat W_0(.)$ are estimators of $m(\cdot, \cdot)$ and $W_0(\cdot)$ respectively; see (ref) below for examples of different estimators. Let $\mathcal{H}_n$ be a sieve parameter space for $h$ (e.g., in this paper we focus on various ANN sieves). We define the optimally weighted SMD estimator $\hat{\alpha}= (\hat{\theta},\hat{h})$ as an approximate solution to \[ \min_{ h\in\mathcal{H}_n,\theta \in \Theta} \widehat{Q}^0_{n}(h,\theta). \] This is an estimator proposed in AC12.
We may analyze the asymptotic properties of this estimator. Since we may view the optimally weighted SMD problem as either a minimum distance program or a sequential GMM estimator, we may carry out two separate analyses of the asymptotic properties. The analysis of the estimator as a minimum distance problem is a specialization of ai2007estimation,ai2012semiparametric,ai2003efficient,CP2015sieve, while (ref) presents a heuristic review of the analysis as a sequential moment restriction, which specializes CL2015gmm. Either approach will lead to the following asymptotic efficient influence function expansion:
\ \
{\bf Riesz Representer.} Lastly, we need to characterize the Riesz representer $v^\star$. The argument in AC03 parametrizes $v^\star = v_\theta^\star (1, -w^\star)$ as a “scale times direction” coordinate. For a fixed scale $v_\theta^\star$, the minimum norm property of Riesz representers implicitly defines the optimal direction $w^\star$ as the following:
Solving the condition \[ \frac{1}{\sigma_0^2}\E[-v_\theta^\star +a(Y_2) \nabla_1 v_h^\star + \Gamma(X) v_h^\star] = -1 \] by plugging in $v_h = -w^\star v_\theta^\star$ then yields \[ v_\theta^\star = \frac{\sigma_0^2}{\E[1 +a(Y_2) \nabla_1 w^\star + \Gamma(X) w^\star]} \quad v_h^\star = \frac{-w^\star \sigma_0^2}{\E[1 +a(Y_2) \nabla_1 w^\star + \Gamma(X) w^\star]} \] as the solutions for the representers where $w^\star$ is defined in (ref) above. If we assume completeness condition then $w^\star = r_0$ as the unique solution to (ref) or (ref) and $v_\theta^\star =(J_0)^ {-1}$.
The consistency, root-$n$ asymptotic normality, consistent variance estimation can all be obtained by directly applying AC03, AC07 for single hidden layer ANN sieves. chen2020nn2 results can be applied for multi-layer ANN sieves.
In this section, we describe in broad strokes the implementation of the eventual estimators for the average derivative of a NPIV, which often involves estimation of nuisance parameters and functions. These nuisance parameters---which often take the form of known transformations of conditional means and variances---require further choice of estimation routines and tuning parameters, details of which are relegated to (ref).
\paragraph{A note on notation.} Recall that we use $Y_1$ to denote the outcome, $Y_2$ to denote variables (endogenous or exogenous) that are included in the structural function, and $X$ to denote exogenous variables that are excluded from the structural function. Certain entries of $X$ and $Y_2$ may be shared. Again, the NPIV model is:
Let $Z = [Y_1, Y_2, X]$ collect the observable random variables (in the population). The parameter of interest is $\theta_0 = \E\bk{\nabla_1 h_0(Y_2)}$, where $\nabla_1 h_0(Y_2)$ is the partial derivative of $h_0$ with respect to its first argument, evaluated at $Y_2$.
We also set up notation for objects related to the sample. Let there be a random sample of $n$ observations. We denote $y_1 \in \R^n, y_2 \in \R^{n \times p}, x \in \R^ {n\times q}$ as vectors and matrices respectively of realized values of the random vector $(Y_1, Y_2, X)$. We will slightly abuse notation and write $f(y_2)$, for a function $f : \R^p \to \R^d$, to be the $(n \times d)$-matrix of outputs obtained by applying $f$ row-wise, and similarly for expressions of the type $f(x)$.\footnote{This notation conforms with how vector operations are broadcast in popular numerical software packages, such as Matlab and the Python scientific computing ecosystem (NumPy, SciPy, PyTorch, etc.).} For a vector valued function $f$, we let $P_f = f(x)(f(x)'f(x))^+ f(x)'$ be the projection matrix onto the column space of $f(x)$.
\paragraph{Quick map of estimation procedures.} We provide a simple map that connects the above model and estimation approaches to the estimators we implement below.
Additionally, we describe the estimator when the analyst is willing to assume more semiparametric structure (e.g. partial linearity) on the structural function $h_0(\cdot)$. We also conclude the section with a brief discussion of software implementation issues.
Consider a linear sieve basis $\phi(\cdot)$ for $X$, where $\phi(X) \in \R^k$. In sample, let $P_\phi$ be the projection matrix projecting onto the column space of $\phi(x)$. For a sample of realizations $v \in \R^n$ of $V$, let $P_\phi v$ be the sample best mean square linear predictor (that approximates the conditional mean) of $v$, since it returns the fitted values of a regression of $v$ on flexible functions of $x$: \[P_\phi v \approx [\E[V_1 \mid X_1],\ldots, \E[V_n \mid X_n]]'.\] Under the NPIV restriction (ref), taking $V = Y_1 - h_0(Y_2)$ and $v = y_1 - h_0(y_2)$, we should expect \[ P_\phi (y_1 - h_0(y_2)) \approx 0. \] This motivates the analogue of the SMD criterion (ref) in the sample, where we choose $h$ so as to minimize the size of the projected residual $P_\phi(y_1 - h (y_2))$:
When the norm chosen is the usual Euclidean norm $\norm{\cdot } = \norm{\cdot}_2$, we obtain the identity-weighted SMD estimator for $h_0$, $\hat h_{\mathrm{ISMD}}$.
Given a preliminary estimator $\tilde h$ for $h_0$, we may form an estimator of the residual conditional variance $\Sigma(X) \equiv \E[ (Y_1-h_0(Y_2))^2 \mid X]$ by forming the estimated residuals $y_{1} - \tilde h (y_2)$ and then projecting $(y_{1} - \tilde h(y_2))^2$ onto $x$, e.g. via the linear sieve basis $\phi (x)$ or via other nonparametric regression techniques such as nearest neighbors. With such an estimator of the heteroskedasticity, we can form a weight matrix $\hat W = \diag(\hat \Sigma(x))^{-1}$. Using the norm $\norm{z}_W^2 \equiv z'Wz$ in (ref) yields the optimally-weighted SMD estimator for $h_0$, $\hat h_{\mathrm{OSMD}}$.
With an estimated $\hat h$ of the structural function $h_0$, we can form two plug-in estimators of $\theta$. The first is the simple plug-in estimator: \[ \hat \theta_{\mathrm{SP}}(\hat h) = \frac{1}{n} \sum_{i=1}^n \nabla_1 \hat h(y_ {2i}). \] See AC (2007) for the root-$n$ asymptotic normality of this estimator, and its asymptotic linear expansion is of the form: \[ \sqrt{n} (\hat \theta_{\mathrm{SP}}(\hat h) - \theta) = \frac{1}{\sqrt{n}} \sum_ {i=1}^n \bk{\nabla_1 h(Y_{2i}) - \theta + \E[v_{h, \text{id}}^\star \mid X] (Y_{1i} - h(Y_{2i})) } + o_p(1). \] where
The simple plug-in estimator does not take into account the covariance between the two moment conditions, $Y_1 - h(Y_2)$ and $\nabla_1(Y_2) - \theta$. The second estimator, the orthogonalized plug-in estimator, orthogonalizes the second moment against the first: \[ \hat \theta_{\mathrm{OP}}(\hat h, \hat \Gamma) = \frac{1}{n} \sum_{i=1}^n [\nabla_1 \hat h(y_ {2i}) - \hat \Gamma(x_i) (y_{1i} - \hat h(y_{2i}))], \] where $\hat \Gamma$ is an estimator of the population projection coefficient of the second moment $\nabla_1 h_0(Y_2) - \theta_0$ onto the first moment condition $Y_1 - h_0(Y_2)$:
One choice of $\hat \Gamma$ is to plug in sample counterparts---plugging in $\hat h$ for $h_0$, plugging in a preliminary $\hat \theta$ (which could be the $\hat\theta_{\mathrm{SP}}(\hat h)$) for $\theta_0$, and plugging in an estimator $\hat \Sigma$ for $\Sigma$---and finally approximate $\E[\cdot \mid X]$ via a linear sieve regression, say with the basis $\phi(\cdot)$.
\ \ \
To summarize, the SMD estimator can be implemented as follows.
\ \ \
{\bf Identity Weighted SMD Estimator of $h(.)$}
{\bf Optimal SMD Estimator of $h(.)$}
{\bf Estimators for $\theta_0$}
Combining simple plug-in with identity-weighted SMD yields the estimation procedure that we term P-ISMD, and combining orthogonal plug-in with optimally weighted SMD yields the estimation procedure that we call OP-OSMD.
We also implement influence function based estimators. As we highlighted in the previous section, one influence function estimator for $\theta_0$ takes the following form
with $\kappa(\cdot)$ defined below. Moreover, given an estimator $\hat h$ for $h$ and $\hat \kappa$ for $\kappa$, we can form the influence function estimator: \[ \hat\theta(\hat h, \hat \kappa) = \frac{1}{n}\sum_{i=1}^n \bk{\nabla_1 \hat h(y_{2i}) - \hat\kappa(x_i)\pr{y_{1i} - \hat h(y_{2i})}}. \]
\paragraph{Identity score estimator (IS)} One influence function, which corresponds to the influence function of the P-ISMD estimator has $\kappa$ taking the following form. We refer to the resulting influence function estimator as IS, for identity score.
\paragraph{Efficient score estimator (ES)} On the other hand, the efficient influence function (ES) uses a different $\kappa(\cdot)$:
where $\Gamma(\cdot)$ is as in (ref), and
are the same as (ref), which are also weighted analogues of the identity weighted $w^\star$, (ref).
The above formulation writes $v^\star$ as a function of $w^\star$; alternatively, we may follow the strategy in (ref) and estimate $v^\star$ directly. One way to estimate the above representer and hence get a feasible score is as follows. Recall that the definition of $v^\star$ is \[ \E[\E[v^\star \mid X] \Sigma(X)^{-1} \E[v \mid X]] = \E[\nabla_1 v + \Gamma(X) v] \quad \norm{v^\star}_{\rho_2}^2 = \sup_v \frac{(\E [\nabla_1 v + \Gamma(X) v])^2}{\E[ \Sigma(X)^{-1} \E[v \mid X]^2 ]}. \] Let $\nu(Y_2)$ be the basis approximating $Y_2$. Suppose we view that $v^\star$ is well approximated by $\nu(Y_2)'\beta$, and that $\E[\cdot \mid X]$ is well approximated by projection onto a basis $\lambda(x)$, then the above definition of $v^\star$ yields a finite-dimensional problem that we may solve in closed form to obtain the following. Consider the following quantities \[ F = \E[\nabla_1\nu(Y_2) + \Gamma(X) \nu(Y_2)] \text{ and } R = \E\bk{\Sigma(X)^{-1} \E[\nu \mid X]\E[\nu \mid X]'}. \] This then implies that $v^\star = \nu'R^{-} F$. In sample, this amounts to
These can then be used to obtain $\hat v^\star$ and the influence function correction term \[ \kappa_{\text{EIF}}(X) = \Gamma(X) - \E[v^\star(Y_2) \mid X] \Sigma(X)^ {-1}. \]
\ \ \
We now discuss how to compute standard errors and confidence intervals---again in broad strokes---for the estimating algorithms P-ISMD, OP-OSMD, IS, and ES. In a nutshell, for the score estimators IS and ES, the estimator $\hat\theta$ is a sample mean of \emph{estimated} influence functions, and its sample variance is directly the properly normalized variance of the influence functions. As a result, under appropriate conditions, a sample variance of the estimated influence functions is consistent for the variance of the influence functions, leading to consistent estimation of standard errors. For the estimators \textsf{IS } and \textsf{ES}, practitioners can therefore compute the standard errors without adjusting for the estimation of the nuisance parameters.
Similarly, estimating the standard errors for the P-ISMD and OP-OSMD estimators amounts to estimating the variance of the influence function. One approach is to simply use the influence function estimates from IS and ES, and leverage the fact that (P-ISMD, IS) and (\textsf{OP-OSMD}, \textsf{ES}) are respectively asymptotically equivalent.
Another approach is to estimate the variance of the influence functions directly, without necessarily estimating the influence functions themselves. The details are stated in (ref), and we may turn the theory into estimators by “putting hats on parameters”: replacing unknown functions with their finite-dimensional sieve approximations, conditional expectation with sieve projections, and expectations and variance with their sample counterparts. For convenience, we reproduce the calculation here:
A third approach, which in our experience seems more accurate than analytic standard errors, is a multiplier bootstrap for the SMD estimators. The bootstrap simply replaces the residual $y_1 - h(y_2)$ in (ref) with the weighted residuals $\omega(y_1 - h(y_2))$ where $\omega = \diag(\omega_1,\ldots, \omega_n)$ are such that $\omega_i \iid F_\omega$, independently of data, for some positively supported distribution $F_\omega$ with unit mean and variance (e.g. the standard Exponential distribution). Given a realization of the bootstrap weights $\omega$, the estimation routines P-ISMD and OP-OSMD would yield an estimate for $\theta$. Repeating this procedure a large number of times would generate a large number of bootstrapped estimates, whose percentiles form confidence interval boundaries.
Assume $h_0$ is partially linear in its first argument, or, additionally, partially additive in subsets of its arguments. Since $h_0$ is linear in its first argument, the slope on that argument is the average derivative $\theta_0$. Therefore, under such a restriction, $h_0$ can be identified with the pair $(\theta_0, \vartheta_0)$ where $\vartheta_0$ is some nuisance parameter governing the rest of the function.
As in the case with SMD estimators in the nonparametric case, we solve the SMD problem (ref), while constraining $\mathcal H$ to conform to the functional form assumptions made. The parameter $\theta_0$ is estimated via direct plug-in, since a solution $\hat h = (\hat\theta, \hat \vartheta)$ for (ref) naturally produces an estimator $\hat\theta$ for $\theta_0$ ai2003efficient.
We now provide a brief recipe on working with neural networks. A feedforward neural network is a composition of layers of the form\footnote{For instance, a ReLU layer is a function of the form \[ x \mapsto \max(0, Wx + b) \] for $W$ a conformable matrix and $b$ a conformable vector.} \[ f_{\sigma, W, b} : \R^m \to \R^n \quad x \mapsto \sigma(W x + b) \quad \text{$\sigma : \R \to \R$ is applied entry-wise.} \] for some conformable matrix $W$, vector $b$, and nonlinear activation function $\sigma$:, i.e. a $k$-hidden-layer neural network has the representation \[ h_\eta : \R^m \to \R^n \quad h = b_{k+1} + W_{k+1} \cdot (f_{\sigma_k, W_k, b_k} \circ \cdots \circ f_ {\sigma_1, W_1, b_1}) \] where we collect the {learnable parameters} $\{W_j, b_j : j =1,\ldots, k+1\}$ as $\eta$. The gradient $\nabla_\eta h_\eta(y_2)$ can be computed efficiently using the celebrated backpropagation algorithm, and, as a result, in practice, neural networks are often optimized via first-order methods such as (stochastic) gradient descent or its variants, such as the popular Adam algorithm kingma2014adam in the machine learning community. Optimization with neural networks is easiest with an unconstrained, differentiable objective, for which numerous computational frameworks exist. We use PyTorch paszke2017automatic in this paper.\footnote{See \url{https://pytorch.org/}} In particular, (ref) is an unconstrained, differentiable objective function, and we may optimize over $\eta$ since the overall gradient may be decomposed into components that are efficiently computed: By the chain rule,\[ \nabla_\eta L(h, y_1, y_2, x) = \nabla_h L \cdot \nabla_\eta h, \] where $L(\cdot,\cdot,\cdot,\cdot)$ denote the objective function (ref).
Compared to conventional numerical linear algebra packages such as NumPy or MATLAB, PyTorch offers two computational advantages particularly suited for deep learning: automatic differentiation and GPU integration. PyTorch tracks the history of computation steps taken to produce a certain output, and automatically computes analytic gradients of the output with respect to its inputs (See (ref) for an example). Autodifferentiation allows gradient descent methods to be carried out conveniently, without the user supplying analytical or numerical gradient calculations manually.
PyTorch also allows arithmetic operations to be computed on GPUs, which have computing architecture that allows for large-scale parallelization of simple operations. For instance, multiplying two $k\times k$ matrices is of order $O(k^3)$ with a naive algorithm, which can be viewed as $k^2$ dot products of size $k$; GPUs allow for parallelized computing of the $k^2$ dot product operations, in contrast to CPUs, where the level of parallelism is determined by the number of CPU cores. For optimization, we use the Adam algorithm kingma2014adam, which is an enhancement of basic gradient descent by estimating higher order gradients.
\Copy{linearsieves}{ We note that even in our ANN implementation, a few nuisance functions are estimated with linear sieves. For instance, the instrument projection, the conditional covariance function $\hat \Gamma(X)$, and $w^\star$ in (ref) are all approximated by linear sieves. It should be in principle possible to use nonlinear sieves, including neural networks, for all of them, and we consider that to be important open work. Here, we detail computational and conceptual difficulties that we have encountered.
First, many nuisance functions take the form of a conditional expectation $\E[r(Y_2) \mid X]$, for some known or unknown function $r(\cdot)$. This is the case, for instance, with $\hat\Gamma$ in (ref), as well as the conditional variance $\hat \Sigma$. In such cases, we can in principle use neural networks to minimize the empirical squared error loss\[ \min_{h_\eta \text{ is an ANN}} \frac{1}{n} \sum_{i=1}^n (\hat r(Y_{2i}) - h_\eta(X))^2 \] and use $\hat h_\eta$ as an estimate. At least computationally, such a procedure makes sense, though its theoretical properties may be delicate. In this work, we avoided using neural networks for $\hat \Gamma$ and $\hat \Sigma$ for computational convenience.
Replacing the instrument projection with neural networks is considerably more challenging for sieve minimum distance. In this case, we are not approximating a function, so much as approximating an operator that projects onto $L^2(X)$. For a given estimate $\hat h$ with estimated structural residuals $Y_{1} - \hat h(Y_2)$, it is not difficult to project it onto $X$ and obtain the estimated projected residuals $\hat r (X; \hat h)$ (by minimizing squared error empirical risk), as well as its squared sample mean $\frac{1}{n} \sum_i \hat r^2(X_i, \hat h)$. However, it becomes challenging to update the neural network weights on $\hat h$. Since the neural network $\hat r$ is trained based on $Y_1 - \hat h (Y_2)$, its weights depend on the weights of $\hat h$ in a complex fashion, and the gradient of $\frac{1}{n} \sum_i \hat r^2(X_i, \hat h)$ with respect to $h$'s weights become computationally intractable. As a result, we could not easily devise a scheme that replaces the instrument projection with nonlinear sieves. Recently, dikkala2020minimax do propose a different method that allows for using nonlinear sieves for the instruments. The performance of this method is compared later in (ref).
Lastly, there are nuisance parameters which are defined through an optimization problem that includes its gradient with respect to its input. The nuisance parameters in the Riesz representer display this property, for instancec $w^\star$ in (ref). To approximate such a parameter, e.g. $w^\star$, with neural networks, we would minimize some criterion function that includes both $w$ and $\nabla_1 w$. To use current off-the-shelf gradient-based training procedures, we would then require the gradient of $\nabla_1 w(\cdot)$ with respect to the $w$'s neural network weights, as well as that of $w(\cdot)$. The former is a niche use-case in deep learning, and so is not well-supported by the autodifferentiation methods in PyTorch. }
\ \ \
We present four Monte Carlo designs in the first subsection. We then describe exactly how we estimated the various components that are needed for the estimators in the next subsection. The last subsection discusses some Monte Carlo results.
We consider a set of Monte Carlo experiments that combine simple but relevant designs that include high dimensional regressors. These designs are also relevant to the kinds of empirical models that are of interests to economists. We describe the four Monte Carlo designs below. A preponderance of our empirical results are based on (ref).
\ \ \
\ \ \
\ \ \
Additionally, we provide a Monte Carlo calibrated to the empirical application in (ref).
Next, we provide a step by step guidance on how to implement the estimators.
\ \ \
We explain here the exact choices of estimators that we used for these Monte Carlo designs. A detailed overview is presented in (ref). Various ANN SMD estimators for $h$ have additional tuning parameters regarding nonlinear optimization, which are described in (ref).
Below, we describe the procedures underlying (ref), which are representative of the procedures in (ref). We also describe the procedures for (ref), which are in turn representative of the procedures in (ref).
Due to the length of the paper, we report representative simulation results in a sequence of figures and tables below.\footnote{\Copy{time}{As a note on computational difficulty, for a single run in estimating OP-OSMD on a sample size of 5000 on a Mac Mini (2020, Apple M1, 8GB RAM), the neural network procedures takes about 33 seconds for (ref). Spline estimation usually takes about a second, as it is equivalent to solving linear IV problems in closed form.} }
\Copy{condexpresults}{For (ref)(a), we compare the performance of a ReLU network (one hidden layer, 40 neurons) with the performance of a spline basis in (ref).\footnote{Like the spline basis we use in other designs, it is a two-knot cubic spline for each variable with pairwise variable interactions of the form $X_i X_j$.} The performance metric we choose is the scaled integrated MSE: \[R^2 = 1- \frac{\sum (\hat f(X_i) - f(X_i))^2}{\min_c \sum (f(X_i) - c)^2}, \] on 1000 out-of-sample data points. $R^2 = 1$ indicates perfect estimation of $f$, and $R^2 = 0$ indicates estimation quality on par with using a constant prediction. We find that for low and moderate noise, neural networks perform better than splines, presumably since it captures more complex interaction patterns in $f$. In the high noise regime, neural network underperforms, as neural networks overfit. In the very-high-noise regime, both estimators overfit and are in fact worse than simply using a constant. We caution that these performance of neural networks results from very minimal tuning, in particular, without validation samples.} \Copy{mcsimpledisc}{Next, for (ref)(b), we show the results in (ref) for P-ISMD estimators, varying over the dimension $p$. We see that, despite the DGP involving a neural network, using neural network estimators only attains a modest improvement over spline estimators, in terms of slightly lower bias, consistent with our findings in the main text.}
The rest of the figures correspond to more difficult (ref) where the first element of $Y_2$ is endogenous ($R_1$). (ref) reports the performance of various ANN SMD estimators for $\theta$ in (ref). The top display plots the results for $n=1000$ and the bottom for $n=5000$. Note here that the columns correspond to various assumptions we maintain on what the econometrician knows about the true structure of $h_0(.)$ in the model $\E[Y_1 - h_0(Y_2)|X]=0. $ The true design is partially additive, and the first column, NP, assumes that the econometrician has no knowledge of the true structure. As we can see, across all implementations (the rows), most of the ANN SMD estimators perform well, which indicates that ANNs seem able to adapt to the unknown structure of $h_0$. The second column labeled PL (for partially linear) assumes that $h_0 (Y_2)$ is partially linear (i.e., $h_0(Y_2) = \theta R_1 + h(R_2, X_2, \tilde X)$) while the third column labeled PA assumes the correct additive structure (i.e., $h_0(Y_2) = \theta R_1 + h_1(R_2)+h_2(X_2)+h_3(\tilde X)$) in the Monte Carlo design is known to the econometrician (but the functions $h_1, h_2, h_3$ within it are of course not known). PA column corresponds to the case where we use ANN sieves to learn all the unknown functions $h_1, h_2, h_3$ although $h_1,h_2$ are functions of scalar random variable. Its performance slightly deteriorates as compared to the NP and PL columns. Notice here that for comparison, the last column for the \textsf{PA} case uses splines to approximate the two scalar valued unknown functions $h_1$ and $h_2$ while $h_3$ is always estimated via ANN (since it is of higher dimensions (at least when $\dim(\tilde X) >0$). We see that the spline results are in line with the \textsf{PL} and \textsf{NP} results, and are adequate here.
In (ref), we report results for the various estimators for (ref), where the unknown function $h_0$ is now nonlinear in the endogenous $R_1$ (the first element of $Y_2$). In the top panel (a), we report results for the case with $R^2/2$ and panel (b) reports results for the case where the unknown function is $ R_1^2/2 + R_1f (X_2)$, where now the derivative depends on the regressor $X_2$ nonlinearly (as the function $f$ is highly nonlinear). Both results are for $n=1000, 5000.$ For panel (a) we see that the spline estimator remain well behaved across all designs (across rows), the single-hidden layer (1L) sigmoid ANN estimators remain adequate while both versions of the AGMM estimators dikkala2020minimax exhibit some bias. In panel (b), spline remains well behaved and so are the ANN estimators. In (ref) we show estimates of the partial derivative evaluated at various fixed values for some regressors. Though the estimators do not track the function well, especially in the tails in the bottom display, the average derivative is estimated well. Interestingly, 3L relu ANN\footnote{This fact seems to be robust to architectural choices.} seems to estimate the derivative function marginally better than splines, perhaps since ANNs are able to automatically generate rich interaction behavior, whereas specifying tensor products for spline sieves is somewhat onerous.
We now examine various implementation choices for (ref). In (ref), we compare various implementations of ANN estimators and spline estimators in (ref). In (ref), we compare identity-weighted estimators (IS, P-ISMD, AGMM, IS-X). Note that P-ISMD and OP-OSMD are the plug in and optimal plug in SMD estimators. In (ref), we compare optimally-weighted estimators that are semiparametrically efficient under suitable regularity conditions (\textsf{ES}, \textsf{OP-OSMD}, \textsf{ES}-X).\footnote{As a reminder, we consider the following estimators: \textsf{IS } or identity weighted score estimator, \textsf{ES } or the efficient score estimator, while \textsf{IS}-X and \textsf{ES}-X are score estimators with two-fold cross fitting. } It is important at the outset to keep in mind that all ANN implementations require some non-negligible tuning as the optimization problem is non-convex and the problem itself with endogeneity, correlation among the regressors, and high dimensions is not easy to tune. Also, currently and for NPIV models, there is no theory for data driven approaches to picking width, depth, or activation functions and finite sample behavior in our design varied \citep*[For linear splines, there are data-driven choice of sieve terms, see][]{chen2021adaptive}.\footnote{although we have not implemented any data-driven choice of spline sieve terms in our paper.} The results across various combinations of $\dim(\tilde X)$ and correlations for $n=1000,5000$ indicate first that ANN OP-OSMD and especially spline estimators seem to behave best. In particular, spline estimators require little tuning and are more stable than all ANN based estimators we use. The SMD ANN estimators are adequate with slight bias for the single-layer, varying-width case. \textsf{IS } and \textsf{ES } ANN estimators are generally less biased and slightly higher variance than \textsf{P-ISMD } and \textsf{OP-OSMD } ANN estimators, but we note that good performance of \textsf{ES } (in the ANN case) is very sensitive to the choice of $\hat\Sigma(X)^{-1}$ in the “optimally-weighted” Riesz representer estimation. (ref) compare the performances of \textsf{ES } in a variety of choices for $\hat\Sigma(X)$. It is interesting that the poor choice of $\hat\Sigma(X)^{-1}$ leads to biased estimation of \textsf{ES } and its cross-fitted versions.
\Copy{calibrateddisc}{Lastly, (ref) examines various estimators on the empirical calibration (ref). Consistent with the findings in other Monte Carlo settings, we generally find that all estimators perform adequately, with similar performance across a variety of neural architectures. We also continue to find that SMD estimators P-ISMD , OP-OSMD have slightly better mean-squared error performance than the score-based estimators in exchange for slightly higher bias. The one exception is ES with variance estimation with only five nearest neighbors, which performs best across the specifications. We conjecture that this is due to how we constructed the residuals in (ref). In particular, we multiply standard Gaussians with estimated residuals $R_i^P, R_i^Y$, which may result in a conditional variance function that is highly non-smooth, and, as a result, a low number of nearest neighbor estimates performs well. In any case, the performance of ES and ES-X continues to show that it is sensitive to $\hat \Sigma(X)$, as is shown in (ref). }
See (ref) for additional Monte Carlo results.
(ref) provides various inference statistics for the ANN SMD estimators P-ISMD and OP-OSMD for (ref), without assuming any semiparametric structure on $h(\cdot)$ beyond smoothness. In particular, we report bootstrapped confidence intervals for ReLU and sigmoid and for depths 1 and 3 when the dimension of the nuisance variables $\tilde{X}$ ranges from 0 to 10. The results are also given for sample sizes $n=1000$ and $n=5000.$ Across all specifications, the two ANN estimators perform adequately.
In (ref), we examine various standard error approaches for a set of estimators in (ref). For each of these, we compute the MC standard deviation, a feasible estimator based on the estimator variance derived from theory, and a bootstrapped standard error. Overall, the theory and bootstrapped standard errors are adequate. In unreported results, criterion (SMD) based bootstrap confidence intervals showed reasonable coverage performance.
Overall, it seems that ANN methods are useful in approximating potentially high dimensional functions in NPIV models. Also, in the class of models we investigated, choices of layers, widths or activation functions are not very consequential in terms of finite sample performance. On the other hand, ANN based estimators in these non-standard NPIV models are hard to tune, and a researcher needs to choose many smoothing parameters. These ANN estimators are also unstable in some runs as they are based on highly complex (and non-convex) optimization programs. In addition, ANNs are not as effective in estimating univariate functions. Finally, to our surprise, We find that various plug-in spline SMD estimators appear stable, less biased generally and can outperform ANNs for NPIV models even in high dimensional cases with 13 continuous regressors.
We present two empirical applications of estimating average derivatives with respect to endogenous price of a nonparametric demand $h_0(Y_2)$ for some non-durable goods. We apply ANN sieves to approximate $h_0(\cdot)$ nonparametrically when its argument $Y_2$ consists of 7 covariates (for gasoline demand) and 6 covariates (for strawberry demand). In the existing literature researchers have used both data sets to estimate unknown $ h_0(\cdot)$ in the model $\E[Y_1 - h_0(Y_2) \mid X] = 0$ by assuming $h$ takes some parametric or semiparametric (such as partially linear) form to avoid the “curse of dimensionality” of $Y_2$. Although served as illustrations, our applications below are the first to estimate the endogenous demand function $ h_0(\cdot)$ fully nonparametrically when $\dim(Y_2)>5$.
We use data on gasoline demand from the 2001 National Household Travel Survey \citep*{blundell2012measuring}. The sample we use include 4,812 observations in the full sample as in chen2018optimal. We estimate an NPIV analogue of the model (11) in blundell2012measuring, $\E[Y_1 - h_0(Y_2) \mid X] = 0$ where $Y_1$ is the log gasoline demand, and $Y_2$ is a vector of 7 random variables consisting of the log gasoline price (possibly endogenous) and the other included covariates following Column (3) in Table 2 of blundell2012measuring. The instrument is the distance from the Gulf coast. We define the estimand as the average price derivative of the unknown structural function $h_0(\cdot)$, which has an average elasticity interpretation. {blundell2012measuring} via OLS, and
(ref) shows our estimates for the average price elasticity (and bootstrapped $95\%$ confidence intervals). Broadly speaking, these estimates point to a similar range of values and are similar to a parametric two-stage least-squares specification. Across estimator classes, the ANN SMD estimates are slightly larger in magnitude than the spline SMD estimates and the ANN IS estimates. Within the ANN SMD estimator class, architecture choices of the networks do not appear to matter much for the result.
We also consider a setting where consumers choose two substitutable goods. We use the Nielsen dataset from compiani2019market,\footnote{Our results do not necessarily represent the views of the Nielsen Company.} where consumers in California choose from strawberries, organic strawberries, and an outside option.\footnote{For a detailed description of the data, see Appendix G of \url{https://www.tse-fr.eu/sites/default/files/TSE/documents/sem2019/eee/compiani.pdf}.} We observe the market share of each type of product, their prices, and a variety of covariates at the market (store-week) level. In the analysis, we consider NPIV model $\E[Y_1 - h_0(Y_2) \mid X] = 0$ where $Y_1$ is the log market share of a type of good (non-organic or organic strawberries) and $Y_2$ is a vector of 6 random variables, including endogenous prices for both types of strawberries and the outside good, and other market-level covariates. The instruments $X$ include Hausman instruments as well as cost shifters such as measurements of consumer taste and income at the market level. We focus on the target parameter $\theta_0=\E[\nabla_1 h_0]$, which is the average derivative of $h$ with respect to the own-price in logs, which we interpret as a version of price elasticity.\footnote{Under a model of the demand where the NPIV condition $\E[Y_1 - h_0(Y_2) \mid X] = 0$ defines the demand function $h_0$, we can understand $\theta_0$ as a price elasticity. However, this model---which implicitly assumes that endogeneity is additive---may not be consistent with microfoundations of consumer behavior berryhaile2016, and so care should be taken in interpreting $\theta_0$ as an elasticity. Nevertheless, for purposes of our illustration here, we may continue to view $\theta_0$ as some well-defined function of the distribution of the data. For a more detailed implementation of demand in this setting, see compiani2019market where in principle one can also use the neural networks based implementation in this paper in a natural way.}
We present the results in (ref). As is perhaps expected from a casual intuition, estimates of $\theta_0$ are negative across both products, and more negative for the more price-sensitive product (organic strawberry). Moreover, results are broadly similar across estimation methods (SMD vs. score) and sieve choices (spline vs. neural net), with perhaps more variability for neural networks in organic strawberries. The estimates for non-organic strawberries hover around $-1.5$, and are reasonably stable across choices of tuning parameters and estimators (IS vs. SMD estimators). The estimates for organic strawberries are more variable across specification of nuisance parameters and neural architectures, but seem to be around $-2$ and $-3$, and larger in magnitude than the own-price elasticity estimate for non-organic strawberries.
These estimates are qualitatively similar to compiani2019market's estimates, which reports median own-price elasticities of $-1.4$ (0.03) for non-organic strawberries and $-5.5$ (0.7) for organic strawberries.\footnote{Interestingly, our estimates are closer to estimates from BLP that compiani2019market reports in Figure 4, which are also around -2 to -3.} Our estimates are more dissimilar for organic strawberries, for which we offer a few conjectures. First, compiani2019market reports estimates following berryhaile2016's approach to demand estimation, that accounts for price endogeneity differently. Under his assumptions, it is possible that our estimator is consistent for a different parameter than his. Second, organic strawberry market shares are very small, and hence fluctuates more on a log scale, thereby resulting in worse estimation precision.
\Copy{conc}{ In this paper, we present two classes of semiparametric efficient estimators for weighted average derivatives (WADs) of nonparametric instrumental variables regressions (NPIV) of moderate and high dimensional endogenous and exogenous regressors. We have conducted detailed Monte Carlo comparisons of finite sample performance of various inefficient and efficient estimators of the WADs using various ANN sieves. The simulation studies and empirical applications confirm the theoretical advantage of ANN approximation of unknown continuous functions of moderately high-dimensional variables, after some tuning of hyper-parameters. Perhaps the most practical findings from our large amount of reported and unreported simulation studies using moderate sample sizes are as follows: the ANN efficient SMD estimators have smaller biases than those of the ANN inefficient SMD estimators, and are less sensitive to the tuning parameters than those of the ANN efficient score estimators. In addition, simple spline based estimators of WADs of NPIVs perform very well in terms of finite sample biases and variances. More research is needed to close the gap between approximation theory and finite sample computational performance in applying flexible ANNs to nonparametric models with endogeneity. }