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.
55,340 characters · 12 sections · 97 citation commands
One-step smoothing splines instrumental regression
\setcounter{equation}{0} We consider the prototypical model
where $Y \in {\mathbb{R}}$ is the dependent variable, $Z \in {\mathbb{R}}$ is the endogenous continuous explanatory variable, and $W\in {\mathbb{R}}^{p}$ are the instrumental variables (IVs). The goal is to estimate nonparametrically $g_0$, the causal effect of the variable $Z$ on $Y$, using $W$ to account for endogeneity. If we assumed linear relationships, we could use the two-stage least squares estimator: in a first stage, one obtains the linear projection of $Z$ given $W$, then in a second stage one linearly regresses $Y$ onto the previously estimated linear projection.
Considering a nonparametric function $g_0$ allows estimating the causal relationship of $Y$ and $Z$ in a more flexible manner. Existing nonparametric estimators of $g_0$ typically rely on two steps. newey2003instrumental develop a nonparametric equivalent to the two-stage least squares estimator: they use linear-in-parameter series expansions of $E [Y | W ]$ and $E[ g(Z) |W ]$ in a generalized method of moments framework, see also ai2003efficient, hall2005nonparametric, blundell2007semi, Johannes2011, horowitz2014adaptive for other series-based methods. Alternatively, Florens2003, hall2005nonparametric, darolles2011nonparametric, and gagliardini2012tikhonov rely on kernel methods to estimate the unknown conditional expectations. As is well-known, backing up a nonparametric estimate of $g_0$ is an ill-posed inverse problem. Hence, one needs some kind of regularization, such as hard thresholding, see horowitz2011applied, chen2012estimation, Tikhonov or ridge-type regularization, see newey2003instrumental, darolles2011nonparametric, florens2011identification, gagliardini2012tikhonov, singh2019kernel, or a Landweber-type iterative method, see Johannes2013, Dunker2014. A general exposition of some of these methods is given by carrasco2007linear. A recent machine learning literature considers solving a saddle point problem that is dual to a generalized method of moments criterion. Here one first maximizes an objective function with respect to a function of the instruments $W$, then one minimizes with respect to a function of $Z$ to obtain $g_0$, see, e.g., bennett2019deep. If the set of functions upon which one optimizes is large, then one has in addition to introduce some penalization in the optimization problem, see dikkala2020minimax, Liao2020. muandet2020dual consider a related but different saddle point problem.
We here develop a smoothing splines instrumental regression estimator for $g_0$ that fully avoids nonparametric first stage estimation. {\color{black} For a random sample $ \left(Y_i,Z_i,W_i\right), i=1, \ldots n$, it is obtained as
where the coefficients $\widehat a_0$, $\widehat a_1$, and $\widehat \delta_i$ $(\text{with }i=1,\ldots,n)$ have a closed-form expression, see Section (ref) for details. Hence, due to its spline nature, our estimator is computationally simple and should be easily appropriated by practitioners. Besides its simplicity,} our estimator has several characteristics that should make it appealing for empirical work. First, our approach is particularly attractive because it is one-step. {\color{black} A key benefit of the one-step nature of our estimator is that it depends upon one regularization parameter only. In existing two-step methods, each stage relies on a particular choice of a smoothing or regularization parameter, whose fine-tuning may be difficult in practice, while affecting the final results. Two-step procedures typically lead to further difficulties, besides choosing one tuning parameter at each step: one needs to estimate in a first step an object that may be more complex than the final object of interest, while first-stage estimation typically affects the second-stage small sample and asymptotic properties.} In some methods, a third parameter is introduced to deal with the ill-posed nature of the inverse problem. To choose our unique regularization parameter, we devise a practical cross-validation method that yields good performances in simulations. {\color{black} {\color{black} Second,} as we show,} our estimator is a natural generalization of the popular smoothing splines estimator, see Wahba1990, green1993nonparametric. The appeal of splines lies in their simplicity together with their excellent approximations of smooth functions, see schumaker_2007. Splines-based methods have been extensively studied, see, e.g., Hall2005, Li2008, Claeskens2009, Schwarz2016, and have been found to have excellent performances in practice. {\color{black} Third, as an additional advantage, one obtains straightforward estimators of derivatives, that are typically of practical interest.}
We also propose some extensions to our method. First, we show how to impose monotonicity constraints by relying on a method proposed by Hall2001. The constrained estimator is simple to implement in practice. Second, to illustrate the versatility of our method, we extend our results to a partly linear model. When used for estimating Engel curves, our smoothing splines estimator and its monotone constrained version yield comparable results, that are reasonable from an economic viewpoint.
The paper is organized as follows. In Section (ref), we detail the {\color{black} construction} of our estimator. We exhibit a global quantity that accounts for all the information contained in Model ((ref)), and that is minimized by the true function $g_0$. We then consider an empirical equivalent, and we set up a minimization problem penalized by a roughness measure of the function to regularize the solution. We show that our estimator extends smoothing splines to the instrumental variables context, and we give a closed form formula for its computation. The asymptotic properties of our estimator are analyzed in Section (ref), where uniform rates of convergences are derived for the function itself and its derivative. {\color{black} We do not however address the issue of rate optimality, see chen2011rate.} {\color{black} Section (ref) focuses on the choice of the regularization via cross-validation.} Section (ref) deals with estimation under monotonicity constraints. In Section (ref), we report selected simulation results, where our estimator exhibits excellent finite sample performance compared to some existing two-step methods, and we illustrate our method for Engel curves estimation. Finally, we extend our estimator to the partly linear model in Section (ref). Concluding remarks are given in Section (ref). Details and supplementary results of simulations, as well as the proof of subsidiary results are included in the online appendix.
\setcounter{theorem}{0} \setcounter{equation}{0}
We assume that $g_{0}$ belong to some space of functions $\mathcal{G}$ on which identification holds, that is,
{\color{black} For a discussion of this condition called {\em completeness}, see newey2003instrumental, DHaultfoeuille2011, and Freyberger2017.} {When $Z$ is continuous, as we assume here, $W$ should typically have at least one continuous component for completeness to hold. Some of the instruments, however, could be discrete, and this will not affect further our exposition and reasoning.}
Instead of dealing directly with ((ref)), as done by most previous work, we consider an equivalent formulation that does not require estimating a conditional expectation given the instruments $W$. {By the results of Bierens1982,
Consider now
where $\mu$ is a symmetric probability measure {\color{black} with support ${\mathbb{R}}^p$}. Then it is straightforward to see that $M(g)\geq 0$ for all $g \in \mathcal{ G}$, and that under ((ref)) \[ M(g)= 0 \Leftrightarrow g = g_{0} \quad \mbox{ a.s.} \] With a random sample $\left\{ \left(Y_i,Z_i,W_i\right), i=1, \ldots n \right\}$ at hand, a natural estimator of $M(g)$ is
{where} \[ \omega(z) = \int_{{\mathbb{R}}^{p}} \exp(i t'z) \, d\mu(t) = \int_{{\mathbb{R}}^{p}} \cos(t'z) \, d\mu(t) \, , \] due to the symmetry of $\mu$: $\omega$ is (up to a constant) the Fourier transform of the density of $\mu$. The above formulation as a V-statistic will be used in practice for computational purposes.
The condition for $\mu$ to have support ${\mathbb{R}}^{p}$ translates into the restriction that $\omega$ should have a strictly positive Fourier transform almost everywhere. Examples include products of triangular, normal, logistic, see JKB95, Student, including Cauchy, see DK2002, or Laplace densities. To achieve scale invariance, we recommend, as in Bierens1982, to scale the exogenous instruments by a measure of dispersion, such as their empirical standard deviation. Note that the function $\omega$ is not similar to a typical kernel used in nonparametric estimation, as there is no smoothing parameter entering $\omega$, which is thus a fixed function that does not vary with the sample size. Hence, our estimation procedure introduces no smoothing on the instruments.}
If $W$ has bounded support, results from Bierens1982 yield that the equivalence ((ref)) holds when $t$ is restricted to lie in a (arbitrary) neighborhood of $0$ in ${\mathbb{R}}^p$. Hence, $\mu$ can be taken as any symmetric probability measure that contains $0$ in the interior of its support. As noted by Bierens1982, there is no loss of generality assuming a bounded support, as his equivalence result equally applies to a one-to-one transformation of $W$, which can be chosen with bounded image.
{\color{black} Our statistic accounts for an infinity of moment conditions, as stated in ((ref)). It is different from a generalized method of moments (GMM) criterion that accounts for an increasing but finite number of moment conditions, as used by e.g. ai2003efficient and chen2012estimation. One can nonetheless relate our estimator to GMM, in a way similar to Carrasco2000. As can be checked, under suitable conditions, the function $\omega$ is a positive definite kernel, that is for any sequence $z_1, \ldots z_n$, and all real numbers $a_1, \ldots a_n$, \[ \sum_{1\leq i, j \leq n }{ \omega(z_i-z_j) a_i a_j} = \int_{}^{}{\left| \sum_{i=1}^{n}{a_i \exp(i z_i^\top t )}\right|^2 \, d\mu(t) }\geq 0, \quad \, . \] Using Mercer's theorem, one can thus write $\omega(t_i - t_j) = \sum_{k = 1}^{\infty}{\lambda_k e_k(t_i) e_k(t_j)} $ for some infinite sequence of nonnegative eigenvalues $\lambda_k$, and an orthonormal basis of functions $\left\{ e_k, k = 1, \ldots \right\}$ in the space of square integrable functions. Hence, \[ M_n(g) = \sum_{k = 1}^{\infty} \lambda_k \left[ \frac{1}{n} \sum_{i=1}^{n} \left( Y_i - g (Z_i)\right) e_k(Z_i) \right]^2 \, , \] which is a GMM criterion that accounts for the countable sequence of moment conditions $E \left[\left(Y - g(Z) \right) e_k(W) \right] = 0$ for which $\lambda_k \neq 0$. As shown in the Supplementary Appendix, if all $W_i$'s are different, then all eigenvalues $\lambda_k$ are strictly positive, and our criterion accounts for an infinitely countable number of moment conditions.}
Minimizing $M_n(g)$ would lead to interpolation. We regularize the problem by assuming some smoothness for the function $g$. We assume that $Z$ has compact support, say $[0,1]$ without loss of generality, and that $\mathcal{ G}$ is the space of differentiable functions on $[0,1]$ with absolutely continuous first derivative.\footnote{{\color{black}For the derivation of our asymptotic results, when $Z$ has a compact support, say $[a,b]$, we can always normalize it on $[0,1]$ by taking the affine transformation $z\mapsto (z-a)/(b-a)$ and obtain the new variable $\widetilde{Z}:=(Z-a)/(b-a)$ supported on $[0,1]$. Then, we can consider the regression function $z\mapsto g_0((b-a)z+ a)$ of $\widetilde Z$ defined on $[0,1]$, and our theoretical analysis remains the same.}} That is, if $g \in \mathcal{G}$, there is an integrable function $g''$ such that $\int_0^z g''(t) \, dt = g'(z) - g'(0)$. We then estimate $g_{0}$ as a minimizer of a penalized version of $M_n(g)$ on $\mathcal{G}$. {\color{black} Specifically,
where $\lambda>0$ is a regularization parameter.}
A recent approach we became aware of when preparing this paper is the “kernel maximum moment loss" approach proposed by zhang2023instrumental. While it does not smooth on the instruments, it assumes that the regression of interest belongs to a Reproducing Kernel Hilbert Space (RKHS), and solves a minimization problem by penalizing by the norm on such a space. The estimator thus depends on the chosen RKHS. Differently, we assume that the regression of interest belongs to a space of smooth functions, and we penalize by the integral of the squared second derivative of the regression, which is a very intuitive measure of roughness, but not a RKHS norm.
We here show that ((ref)) has a unique solution, a natural cubic spline, that we characterize in Proposition (ref) below. We begin with some definitions. {\color{black} For $a < t_1 < \ldots < t_n < b$, a function $g$ on $[a,b]$ is a cubic spline if two conditions are satisfied: on each of the intervals $(a, t_1), \ldots (t_n, b)$, $g$ is a cubic polynomial; the polynomial pieces fit together at each $t_i$ in such a way that $g$ and its first and second derivatives are continuous. The points $t_i$ are called knots.} A cubic spline on $[a,b]$ is said to be a {\em natural cubic spline} if its second and third derivatives are zero at $a$ and $b$. Without loss of generality, we consider hereafter that $[a,b] = [0,1]$.
Given any values $(g_i, Z_i), \ i=1, \ldots, n$, $n \geq 2$, there is a unique interpolating natural cubic spline, that is, a natural cubic spline $g$ with knots {\color{black} at observations} $Z_i$ such that $g(Z_i) = g_i, i = 1,\ldots n$. For details, see e.g. green1993nonparametric. A key result for our analysis is the following.
This result allows us to restrict our attention to natural cubic splines, when studying the potential minimizers of {\color{black} $S_n(g)$}. Indeed, suppose $\widetilde{g}$ is any function in ${\cal G}$ that is not a natural cubic spline with knots at $Z_i$. Let $g$ be the natural cubic spline interpolant to the values $\widetilde{g}(Z_i)$. Then $M_n\left({g}\right) = M_n\left(\widetilde{g}\right)$. Because of the above optimality property of the natural cubic spline interpolant, ((ref)) holds with strict inequality, and thus $ S_n(\widetilde{g}) > S_n(g)$. This means that, unless $\widetilde{g}$ itself is a natural cubic spline with knots at $Z_i$, we can find a natural cubic spline with knots at $Z_i$ that attains a smaller value of $S_n(g)$. It follows at once that a minimizer $g$ of $S_n(g)$, if it exists, must be a natural cubic spline. It is key to notice that we have not forced $g$ to be a natural cubic spline. This arises as a mathematical consequence of the choice of the roughness penalty. Now, as detailed below, we only need to minimize $S_n(g)$ over a finite-dimensional class of functions.
{Assuming the $Z_i$'s are all different, which happens with probability one for a continuous $Z$, a natural cubic spline with knots at $Z_i$ can be written as}
{see green1993nonparametric}. The function $g$ is uniquely defined by the coefficients $a_0$, $a_1$, and $\delta_i, i = 1, \ldots n$, or equivalently by its value at the knots, see Proposition (ref)'s proof for details.
It will be useful for what follows to use matrix notations. Let \[ \bm{Z}= \left(
\right) \, , \] $\bm{E}(Z_1, \ldots Z_n) = \bm{E} = \left[ \frac{1}{12} |Z_{i}-Z_{j}|^{3} , i,j=1,\ldots n\right]$, and $\bm{g} = \left(g(Z_1), \ldots g(Z_n)\right)^T$. Then $\bm{g} = \bm{Z a} + \bm{E \delta}$ with constraints $\bm{Z^T \delta} = 0$. Also, one can check that \[ \int g'' (z)^{2} \, dz = \bm{\delta^T {E} \delta} \, , \] see green1993nonparametric. Let $\bm{Y}$ be the vector $\left(Y_1, \ldots Y_n\right)^T$, then
where $\bm{\Omega}$ is the matrix with generic element $n^{-2} \omega(W_i-W_j)$. Hence, we want to minimize a quadratic function in parameters under the constraints $\bm{Z^T \delta} = 0$. This yields a unique solution under the usual requirements. The following proposition gives a precise characterization.
Our estimator is obtained directly by solving the linear system of equations ((ref)). It does not necessitate estimation of other nonparametric quantities, and relies on only one regularization parameter $\lambda$. It also directly provides an estimator of the first derivative of $g$ as
There are alternative ways to ((ref)) for expressing a natural cubic spline. We focus on this formulation as it does not rely on a particular support of $Z$, nor on the fact that the knots $Z_i$ are arranged in increasing order. In particular, the closed-form expression in Proposition (ref) is valid regardless of the support of $Z$ and therefore it can be used without first transforming $Z$ into $(0,1)$. We also found this formulation to be convenient for practical implementation. For large samples, where the above formula may not be computationally efficient, one can adapt to our context the Reinsch algorithm, see green1993nonparametric.
{\color{black} We note that it would be possible to generalize our estimator ((ref)), say by penalizing by the integral of the square of the third derivative, so as to obtain a natural quartic spline, see Wahba1990. Such a generalization would entail more technicalities without changing the wide picture. In practice, computational methods for storage of and computation based on natural cubic splines are well developed, see e.g. green1993nonparametric and references therein. Therefore, in our work, we focus on this case.}
\setcounter{theorem}{0} \setcounter{equation}{0}
The formal study of our estimator is based on a reformulation of $M(g)$ in ((ref)). Let $\mathcal{D}^2$ be the set of twice weakly differentiable functions. Consider \[ \mathcal{G}=\left\{ g:[0,1] \to {\mathbb{R}}, g \in \mathcal{D}^2 : \int_{0}^1 |g''(t)|^2 \, dt < \infty \right\}, \, \qquad \mathcal{H}=\left\{ h \in \mathcal{G}: h(0) = h'(0) = 0 \right\}, \, \] and the inner product $\left<h_1,h_2\right>_{\mathcal{H}}=\int_{0}^1 h_1''(z)\,h_2''(z)\, dz$ on $\mathcal{H}$.\footnote{\color{black} The considered inner product on $\mathcal{H}$ is standard in the spline literature, see Wahba1990. Other work where penalization on derivatives is used also considers inner products implying derivatives, see e.g. florens2011identification or florens2018nonparametric. Note there is a one-to-one correspondence between a function $h \in \mathcal{H}$ and its second derivative.} Each $g\in \mathcal{G}$ can be uniquely written as $g(z)=(1,z)\beta+h(z)$, where $\beta=(g(0), g'(0))\in \mathbb{R}^2$, $h(z)=g(z)-g(0)-g'(0)z$, $h\in\mathcal{H}$. Denote by $L^2_\mu $ the space of complex functions $l$ from $\mathbb{R}^q$ onto $\mathbb{C}$ such that \[ \|l\|^{2}_\mu = \int_{}^{}{|l(t)|^2 \, d\mu(t)} < \infty \, . \] Consider the operators $\operatorname{A}:\mathcal{H}\mapsto L^2_\mu $ and $\operatorname{B}:{\mathbb{R}}^2\mapsto L^2_\mu $ such that
The minimization problem ((ref)) identifying $g_0$ can be expressed as
for $(\beta,h)\in\mathbb{R}^2\times\mathcal{H}$. The above quantity reaches its minimum zero at $(\beta_0,h_0)$, with $g_0(z)=(1,z)\beta_0+h_0(z)$. A key advantage of this formulation for theoretical analysis is that using orthogonal projections, we can profile ((ref)) to first determine ${h}_0$, then ${\beta}_0$ as a function of ${h}_0$. In our proofs, we will also consider the penalized empirical counterpart of ((ref)) and use a similar profiling method to obtain $\widehat g(z)$. The following assumption ensures that $E[Y\exp(\mathbf{i}W^T\cdot)]\in L^2_{\mu}$, and that $\operatorname{A}$ and $\operatorname{B}$ are valued in $L^2_{\mu}$.
Our assumptions on the support of $Z$ is without much loss of generality, since we can always use a one-to-one transformation that maps $Z$ into $[0,1]$. We then formalize the completeness assumption, under which the problem ((ref)) admits a unique solution $(\beta_0,h_0)$.
We now introduce a {\em source condition}, which is common in the literature on inverse problems. While it is not needed to establish the consistency of $\widehat g$ and its first derivative, it is necessary to obtain convergence rates.
\textcolor{black}{ As noted in darolles2011nonparametric, the $\gamma$ in our source condition is related to (i) the smoothness of the function $g_0$, as measured by the rate of decay of the Fourier coefficients $(\left<g_0,\varphi_j\right>_{\mathcal{H}})_j$; and (ii) the degree of ill-posedness of the inverse problem, as measured by the rate of decay of the singular values $(\sigma_j)_j$. Again following darolles2011nonparametric, the singular values $(\sigma_j)_j$ are related to the link between the endogenous regressor and the instruments. For example, in the extreme case where the endogenous regressor and the instruments are independent, there is at least one null singular value. }
We obtain consistency of our estimator and its derivative under mild assumptions, that only involve a standard condition on the regularization parameter $\lambda$. \color{black}The role of the condition $n\lambda\to \infty$ is to ensure that the variance term vanishes, while the restriction $\lambda\to 0$ guarantees that the bias goes to 0.\color{black} \ {\color{black}By contrast, in} two-step estimation methods that smooth over the instruments, one has to ensure that first-step estimation is consistent, and one typically needs conditions that relate the different smoothing parameters, see e.g. ai2003efficient, chen2012estimation. In some instances, consistency may further necessitate regularization parameters, see chen2012estimation. A general discussion can be found in carrasco2007linear.
Turning now to our rates of convergence, we do not claim that these are sharp. {\color{black} As noted by a referee, our rates cannot be minimax, in the sense of chen2011rate, since we are not bounding away from infinity the norm of functions in ${\cal G}$. hall2005nonparametric and chen2018optimal also obtain estimators with minimax rates. It is, however, unclear how to compare our results to these optimal rates, which are in $L^2$ norm and hold under a different set of assumptions. In the Supplementary Appendix, we study the $L^2$ norm of our estimator and its derivatives under a set of assumptions that resemble the ones of chen2011rate, namely, we rely on Hilbert scales and a link condition.
By contrast to previous results in this literature, our rates depend upon only one smoothing parameter. Our assumption, though, assumes that the problem is mildly ill-posed, while some previous work also considers the case of severely ill-posed inverse problems. As apparent from our proofs,} the rate $1/\sqrt{n\lambda}$ corresponds to a standard deviation term, while the second rate $\lambda^{\frac{\gamma \wedge 2}{2}}$ corresponds to a bias term. If $\lambda$ is chosen to balance these two rates, we obtain the convergence rate $n^{- \frac{\gamma \wedge 2}{2\left(1+ \gamma \wedge 2\right)}}$. For $\gamma = 2$ or 1, this respectively yields $n^{-1/3}$ and $n^{-1/4}$.
{\color{black}
We propose to choose our single penalty parameter $\lambda$ via cross-validation, where the criterion to be minimized is akin to $M_n(g)$. The simplest procedure is based on random sample splitting, where an estimator $\widehat{g}_\lambda$ is obtained from the first fold of size $n_1 = \lfloor n/2\rfloor$ (the largest integer smaller than $n/2$), and the second fold of size $n_2= n - n_1$ is used for validation. The criterion
is thus computed and minimized over a fixed set $\Lambda=[\underline{\lambda}_n,\overline{\lambda}_n]\subset \mathbb{R}_+$ of penalty parameters, where $(\underline{\lambda}_n)_n$ and $(\overline{\lambda}_n)_n$ are two sequences such that $\overline{\lambda}_n\rightarrow 0$ and $n \underline{\lambda}_n\rightarrow \infty$. A more elaborate procedure consists in letting each fold play the role of the training fold in turn. Specifically, we compute $\widehat{g}_{k, \lambda}$ for each fold $k = 1,2$, and we create the cross-validated criterion \[ \int \left| \frac{1}{n} \left( \sum_{i=1}^{n_1}{ \left( Y_i - \widehat{g}_{2,\lambda} (Z_i) \right) \exp(\mathbf{i}W_i^\top t) } + \sum_{i=n_1 + 1}^{n}{ \left( Y_i - \widehat{g}_{1,\lambda} (Z_i) \right) \exp(\mathbf{i}W_i^\top t) } \right) \right|^2 \, d\mu(t) \, . \] This two-fold cross-validation method can easily be extended to k-fold cross-validation, or repeated k-fold cross-validation. In empirical implementations, we used the above two-fold cross-validation method with good results.
For such a data-driven procedure, we here provide some theoretical rationale, similar to the ones put forward by, e.g., darolles2011nonparametric. Technical details can be found in the supplementary material. To simplify exposition, let us consider sample splitting (similar arguments could be applied to k-fold cross-validation schemes). Then
Under suitable assumptions, using simple arguments, one can show that the first term is $O_p (n^{-1/2})$ in $\|\cdot\|_\mu$ norm. For the second term, one can show that \[ \left| \frac{1}{n_2} \sum_{i=n_1+1}^{n}{ \left( {g}_0 (Z_i) - \widehat{g}_\lambda (Z_i) \right) \exp(\mathbf{i}W_i^\top t)} \right| \leq \sup_z |g_0(z) - \widehat{g}_\lambda(z)| = O_P \left( \frac{1}{\sqrt{n\lambda}} + \lambda^{\frac{\gamma \wedge 2}{2}}\right) \] uniformly over $\Lambda$. Hence, $M_{n_2}(\widehat{g}_\lambda) = O_P \left(\frac{1}{{n\lambda}} + \lambda^{\gamma \wedge 2}\right)$ uniformly over $\Lambda$, which is the same order as $\sup_z |g_0(z) - \widehat{g}_\lambda(z)|^2$.
Minimizing our criterion with respect to $\lambda$ should thus give a penalty parameter that yields the appropriate balance between bias and variance. This argument is similar to the one put forward in previous works on data-driven choices of penalty parameters in Tikhonov nonparametric instrumental regressions, see feve2010practice,feve2014non,Carrasco2014, and also engl1996regularization for a Tikhonov estimator with a known operator. These works similarly show that their criterion used to select the penalty parameter is an $O_P$ of the same order as the rate of convergence of their estimator, and argue that optimizing their criterion with respect to the regularization parameter should be roughly equivalent to optimizing the convergence rate of their estimator. A more detailed analysis would focus on the asymptotic properties of the estimator that uses the above data-driven regularization parameter, but this is outside the scope of the current paper.}
{
\setcounter{theorem}{0} \setcounter{equation}{0}
In some instances, we may expect the function of interest $g_0$ to be monotonic. If $g_0$ is the Engel curve that relates the proportion of consumer expenditure on a good as a function of total expenditure, we typically expect this function to be increasing for a “normal” good and decreasing for an “inferior” good. Accounting for monotonicity in estimation is expected to improve accuracy in small and moderate samples, {\color{black}see mammen1999smoothing} and chetverikov2017nonparametric.
To implement such a monotonicity constraint into estimation, we note that since our smoothing splines estimator is linear, the derivative estimator ((ref)) is linear as well. Let us express it in matrix form. We can write $ \bm{g'} = \bm{O a} + \bm{D \delta} $, where $\bm{g'} = \left(g'(Z_1), \ldots g'(Z_n)\right)^T$, $\bm{D} = \left[ \frac{1}{4} \operatorname{sign}(Z_i - Z_j) |Z_{i}-Z_{j}|^{2} , i,j=1,\ldots n\right]$, and \[ \bm{O} = \left(
\right) \, . \] From Proposition (ref),
We rely on a method proposed by Hall2001, that is based on the same linear estimator but reweights the observations $Y_i$ to impose monotonicity at observations points. It adjusts the unconstrained estimator by tilting the empirical distribution to make the least possible change, in the sense of a distance measure, subject to imposing the constraint of monotonicity at observation points. Specifically, if $g_0$ is assumed to be monotonically increasing, we consider the constrained optimization program
where $\bm{p} \circ \bm{Y} = (p_1 Y_1,\ldots, p_n Y_n)^T$ is the Hadamard product between vectors. If $g_0$ was assumed to be monotonically decreasing, we would modify the last inequalities. Hall2001 considered more general optimization problems based on a family of Cressie-Read divergences, but we focus on the above program for convenience. It is strictly convex, so it admits a unique solution $\bm{p^*}$, and it is computationally fast to solve. The final estimator ${\widehat{g}^*}$ is the natural cubic spline with coefficients $\bm{a}^*$ and $\bm{\delta}^*$ defined as in ((ref)), with $\bm{p^*} \circ \bm{Y}$ in place of $\bm{Y}$.
{\color{black} To obtain the asymptotic properties of our constrained smoothing splines estimator, we state the following supplementary assumption.}
{\color{black} The results from Theorem (ref) then directly apply to our constrained estimator.} Indeed, as $\widehat{g}'$ is uniformly consistent, the constraint in the optimization problem ((ref)) becomes asymptotically irrelevant from Assumption (ref). Accordingly, $\widehat g^*=\widehat g$ with probability approaching one, and our results readily follow. While the monotonicity constraints become asymptotically irrelevant, they can matter in finite samples, as shown by chetverikov2017nonparametric and illustrated by our empirical results.}
\setcounter{equation}{0} \setcounter{theorem}{0}
We used a DGP in line with Equation ((ref)), where
and $(W,V,\eta)$ are independent standard Gaussians. This yields standard Gaussian marginal distributions for $\varepsilon$ and $Z$ whatever the values of the parameters. The correlation $\rho_{\varepsilon V}$ measures the level of endogeneity of $Z$. The correlation $\rho_{WZ}$ measures instead the strength of the instrument $W$.
We implemented our smoothing splines estimator with $\omega$ equal to the density of a Laplace distribution with mean zero and variance 1. The choice of the penalty parameter $\lambda$ was based on two-fold cross-validation, as detailed {\color{black}in Section (ref)}. We considered $\lambda$ within the grid $\{p/(1-p),p=10^{-5}+k*(0.7-10^{-5})/399,k=0,\dots,399\}$.
\color{black}
We compared our estimator to two existing methods for which a data-driven procedure has been proposed for the choice of smoothing {or} regularization parameters. We considered first \textcolor{black}{the Tikhonov estimator penalized by the Sobolev norm of degree $2$ of gagliardini2012tikhonov}, hereafter referred as Tikhonov. We also considered the series estimator of horowitz2014adaptive based on a basis of Legendre polynomials. The implementation details of both methods are given in the online supplement, together with supplementary results.
We first considered two functional forms for $g_0$, each normalized to have unit variance: a quadratic function $g_{0,1}(z) = {z^2}/{\sqrt{2}}$, and a non-polynomial function $g_{0,2}(z) = \sqrt{3\sqrt{3}}\,z\,\exp(-z^2/2)$. We ran 2000 Monte Carlo simulations with sample sizes $n=200$ and $400$. We consider three couples of values for $(\rho_{\varepsilon V},\rho_{WZ})$: (a) $(0.5,0.9)$, a setting with low endogeneity and a strong instrument, (b) $(0.8,0.9)$, corresponding to high endogeneity and a strong instrument, (c) $(0.8,0.7)$, a more complex setting with high endogeneity level but a weaker instrument. To evaluate the gains of imposing monotonicity, we then considered a third function $g_{0,3}(z) = (\sqrt(10/3) \log(|z-1|+1) \operatorname{sign}(z-1) - 0.6 z+ 2 z^3)/8$. The regularization parameter $\lambda$ was chosen before the monotonizing step, and we used the {\tt R} package CVXR to solve ((ref)), see Fu.
Table (ref) reports our results. {\color{black} The series estimator is severely biased in all cases, but for the non-polynomial $g_{0,2}$ and strong instruments, while our estimator is always almost unbiased. The Tikhonov estimator mostly lies in between, but with large differences depending on the functional form of $g_0$. In terms of variance, the series estimator does slightly better than smoothing splines in most cases, that itself does better than Tikhonov, but for the non-polynomial $g_{0,2}$. Smoothing splines performs best in terms of MSE in almost all cases. Exceptions are cases corresponding to the non-polynomial $g_{0,2}$ with $n=400$ and strong instruments, where the series estimator is close to unbiased. Overall, the severity of endogeneity does not affect much the estimators' performances. When the instrument's strength decreases, our smoothing splines estimator appears to perform best. Finally, imposing monotonicity to our smoothing spline estimator does not affect much bias, but yields a substantial decrease in variance, as expected. Depending on the particular setup, it can be more than halved.}
We applied the smoothing spline estimator to the estimation of Engel curves, which relate the proportion of spending on a given good as a function of total expenditures. We used the “Engel95” dataset engel95, from the {\tt R} package {\tt np}, see hayfield2008. This dataset is a random sample from the 1995 British Family Expenditure Survey and contains data for 1655 households of married couples for which the head-of-household is employed and between 25 and 55 years old. We focused on the subsample of 628 households with no kids. We report results for two Engel curves, pertaining to the expenditure shares on {leisure and fuel}. Economic theory suggests that the Engel curve for leisure is increasing and the one for fuel is decreasing. Following blundell2007semi, we instrumented the logarithm of total household's expenditure, which is likely endogenous, by the logarithm of total earnings before tax. We consider the four estimators used in our simulations, and implementation details remain the same.
The estimated nonparametric functions are reported in Figure (ref). The Tikhonov estimate exhibits a non-monotonic and quite irregular behavior, while the series estimate is mainly monotonic and very regular. Since our smoothing splines estimates are monotonic, but at the boundaries of the data, our constrained and unconstrained estimates are very close. Both are in line with the findings of blundell2007semi.
\setcounter{equation}{0} \setcounter{theorem}{0}
We have proposed a generalization of regression smoothing splines to the context where there is endogeneity and instrumental variables are available. While we detail our estimator and its properties in the simple univariate context, a multivariate extension could be considered. However, including more covariates in a fully nonparametric way would submit us to the curse of dimensionality typical of functional estimation. Hence, we focus here on a partly linear model, as considered by e.g. heckman1986Spline, Robinson1988, blundell2007semi, Chen2009, and florens2012instrumental. This model provides a simple and economical way to include additional controls. We thus consider
where ${X} \in\mathbb{R}^q$ is a vector of exogenous covariates, whose components are thus included in $W$, while, as earlier, $Z\in\mathbb{R}$ is the endogenous variable. The following condition ensures identification of $({\gamma}_0,g_0)$ in ((ref)).
Our identification assumption is similar to ones imposed in other work, see, e.g., chen2012estimation or florens2012instrumental. First, it excludes collinearity between the components of $X$. Second, it requires that the distribution of $Z$ given $W$ must be complete so that the mapping $g\mapsto E[g(Z)|W]$ is injective. Third, it rules out the presence of an intercept in $X$, since an intercept can always be absorbed by the nonparametric function $g_0$. Fourth, it requires that no function $g$ is such that $E[g(Z)|W]$ is a linear function of the variables in $X$. Under Assumption (ref),
We proceed as in the benchmark model and estimate $(\gamma_0,g_0)$ by minimizing the empirical counterpart of the above criterion penalized by a roughness measure of the nonparametric component, that is
The estimators $(\widehat \gamma,\widehat g)$ can be computed in the same way as in the benchmark model. To show this, let $\bm{L}$ be the $n\times (q+2)$ matrix whose row $i$ is $(1,Z_i,X_i^T)$. The following is a direct extension of Proposition (ref) to the partly linear model.
We now focus on the consistency and the convergence rates of the regression function. Let us define the operator $\operatorname{D}:\mathbb{R}^{q+2}\mapsto L^2_\mu$ such that
The following is a direct extension of Theorem (ref).
\color{black}
Further extensions to our method could be considered, such as a nonparametrically additive model, see Linton1995, \color{black} a more detailed theory for the choice of the regularization parameter, and inference results on the nonparametric function of interest or some functional of it.
Concerning the data-driven selection of $\lambda$, breunig2016adaptive,breunig2020adaptive, and chen2024adaptive have derived theoretical results for a Lepski-type selection procedure for series nonparametric instrumental variables estimators without penalization. An interesting avenue for research is to extend such results to Tikhonov estimators such as ours. jansson2020towards provides high-level conditions for the Lepski's method for regularized estimators. One could try to verify these conditions for our approach.
Regarding inference results, most of the literature has focused on non-penalized estimators for nonparametric instrumental variable regressions. Two exceptions are chen2015sieve, which contains inference results on functionals of a penalized series estimator, and babii2020honest, which develops honest and uniform confidence bands for a Tikhonov estimator. The results in chen2015sieve need to be significantly adapted. Indeed, the latter paper considers a finite but growing number of moments while we consider infinitely many moments. Moreover, the penalized sieve minimum distance estimator studied chen2015sieve combines two types of regularizations: slowly growing finite-dimensional sieves and penalization, while we only rely on the latter. Concerning the approach of babii2020honest, we consider a different Tikhonov estimator with other operators, so the proofs must also be adapted. \color{black}