Identification, Estimation and Inference Based on Structural Error Projection

Introduction

Identifying structural parameters when regressors are correlated with unobserved disturbances is a central problem in econometrics. Instrumental variable (IV) methods provide a widely used solution, but they rely on the availability of valid instruments that are correlated with the endogenous regressors but uncorrelated with the structural errors. In many empirical settings such instruments are difficult to obtain or justify. This paper proposes a semiparametric instrumental variable approach that constructs instruments directly from the observed regressors. The key observation is that endogeneity arises from the component of the structural error that is correlated with the regressors. By projecting the structural equation onto the sigma--field generated by the regressors, we decompose the structural error into a conditional expectation component and an orthogonal residual. This decomposition yields an exogenous regression representation that allows us to construct valid instruments as functions of the observed covariates.

Suppose that we have a linear structural model

where $\beta _{0}$ are the parameters of interest and $\mathbf{x}$ is the observed $d$-dimensional vector of variables related to the outcome $y$. The main issue we face regarding identification and estimation is that ${\mathbb{ E}}[\mathbf{x}\varepsilon ]\neq 0.$ We write $\varepsilon ={\mathbb{E}} [\varepsilon |\mathbf{x}]+e,$ where $m(\mathbf{x)}={\mathbb{E}}[\varepsilon | \mathbf{x}]$ is of unknown functional form and ${\mathbb{E}}[e|\mathbf{x}]=0$ . We suppose that $m\in S,$ a subspace of the Hilbert space of square integrable functions of $\mathbf{x}$. Letting $\mathbf{P}_{\mathcal{S}}$ denote the projection operator on to the subspace $\mathcal{S}$, we define instruments $\mathbf{z}=\mathbf{x}- \mathbf{P}_{\mathcal{S}} \, \mathbf{x}$, and then apply the IV method to estimate the structural parameters. Because $ \mathbf{z}$ is by construction orthogonal to the space containing $m,$ it follows that ${\mathbb{E}}[\mathbf{z}\varepsilon ]=0,$ so that $\mathbf{z}$ serves as a valid instrument for $\mathbf{x}.$ Identification then follows from the condition that

under which the structural parameter satisfies $\bm{\beta} _{0}=\bm{\Sigma} _{{x}}^{-1} \bm{\Sigma} _{xy},$ where $\bm{\Sigma} _{xy}={\mathbb{E}}[\left( \mathbf{x}- \mathbf{P}_{\mathcal{S}}\mathbf{x}\right) y].$

The first condition $m\in {\mathcal{S}}$ specifies the admissible sources of confounding. The second condition $\bm{\Sigma} _{x}>\mathbf{0}$ requires sufficient variation in the structural regressors after those sources have been removed. We develop estimation technology and inference methods for the structural parameters. A key consideration is the choice of ${\mathcal{S}}$, which could be justified by economic, institutional, or design information. Plausible examples include: (a) a known collection of fixed effects or seasonal components; (b) a symmetry or invariance restriction; (c) specified interactions that exclude the structural linear index; (d) a finite-dimensional nuisance space motivated by the empirical design. We explore some examples below. In the absence of such concrete choices, we propose to use the LASSO method applied to a dictionary of basis functions to parsimoniously represent this space.

Our approach is related to several strands of the literature that develop identification strategies derived from instruments internally generated from features of the data generating process. RR2003 shows that heteroskedasticity across regimes can identify structural parameters in simultaneous equations models even in the absence of conventional instruments, see also sp2001. In that framework, changes in the variance of structural shocks generate additional moment conditions that identify the structural coefficients. AL2012 proposes another approach for simultaneous equations based on heteroskedasticity, showing that covariance restrictions between regressors and structural disturbances can generate valid instruments even when external instruments are unavailable.

These methods exploit variation in second moments of the data generating process to produce internal instruments. lewis2025 points to a related literature in macroeconometrics that exploits non-Gaussianity in the structural shocks to obtain identification, such as the independent component analysis (ICA). The literature on endogeneity in general is of course very extensive. There are several well known survey articles by swy2002, anderson2005, sy2005, Imbens2014, ass2019, cr2020, Angrist2023, and lewis2025, and many results have been written into textbooks, such as jmw2001, wooldridge2016, greene2018, and jsmw2018.

The approach developed in this paper follows a similar philosophy but relies on a different source of identification. Instead of exploiting heteroskedasticity or regime variation, we use a projection-based decomposition of the conditional mean of the structural error. This decomposition allows us to construct instrumental variables directly as functions of the regressors. Identification therefore arises from the structure of the conditional expectation rather than from variance shifts or heteroskedasticity--based moment conditions.\footnote{ However, under joint Gaussianity, the endogeneity will be linear and our identification strategy fails as well.}

Our approach is also related to the control function literature discussed, for example, in npv1999, np2003, blundell2004, bck2007, IN2009, newey2013, and JMW2015. These existing methods address endogeneity in semiparametric models by introducing auxiliary instruments or structural restrictions. By contrast, the proposed SP method constructs instruments internally through an orthogonal decomposition of the structural error. This yields a closed--form instrumental variable in each case, which is termed as internally “Constructed Instrumental Variable" (CIV), that depends only on observed regressors and allows identification in a number of classes of linear and nonlinear models as discussed in Appendices (ref) and A.3 below.

In summary, the main distinctions and contributions of this paper are given follows.

\noindenta: \, We propose using a semiparametric projection (SP) method to derive the so--called “CIV" to correctly identify and consistently estimate the unknown parameters and functions involved in a class of correctly identified models under study;

\noindentb: \, We then show that the SP method also offers an alternative way to those methods available from the existing literature about how to deal with endogeneity issues involved in certain classes of weakly identified models;

\noindentc: \, The proposed SP method is invariant to the degree of endogeneity, including a wide range of weak endogeneity, involved in many classes of nonlinear models;

\noindentd: \, We establish the corresponding estimation theory to show that the proposed SP method produces consistent and unbiased estimates as well as asymptotically normal distributions;

\noindente: \, We develop a Hausman--type of nonparametric statistic for testing a full level of strong exogeneity versus a wide range of weak endogeneity with a near--optimal rate;

\noindentf: \, The proposed SP method offers a closed--form expression for an IV function of the data for the correctly identified model under study, and it is easily implementable; and

\noindentg: \, We propose a simple LASSO selection method to evaluate and implement the estimation method by simulated and real data examples.

The organization of this paper is as follows. Sections (ref) and (ref) respectively propose and discuss how to identify and estimate the parameters of interest for a class of correctly identified linear models. Section (ref) then shows that the proposed SP method can also be applied to a class of weakly identified linear models, including such linear models associated with weak endogeneity. Section (ref).1 proposes a LASSO selection method for an optimal set of orthonormal series before it is demonstrated in an analysis of the return to schooling case study in Section (ref).2. Section (ref) concludes and discusses the applicability of the SP method to several classes of extended models.

Appendix (ref) provides some heuristics about the identification and estimation proposed in Sections (ref) and (ref). Appendix (ref) shows that the proposed SP method can be generalized to deal with a wide class of nonlinear models. Appendix A.3 outlines how to extend the proposed SP method to two classes of non-- and semi--parametric regressions, and one class of nonlinear and non--separable binary models. Appendix A.4 briefly discusses other estimation issues including a type of GMM estimation.

Appendix B.1 discusses some implementation issues. Appendix B.2 evaluates the finite--sample properties of the SP method using a number of simulated datasets. Appendix B.3 further investigates the empirical example considered in Section (ref).2 and then adds an extra empirical dataset. Appendix B.4 proposes a jackknife bias--correction method. Appendix C collects the proofs of the main results. Mathematical proofs for Appendix (ref), and additional simulations and extensions, are given in Appendices D-F of an online supplementary document.

For notational consistency, we use $\mathbf{x}$ to stand for a vector regardless of whether it is stochastic or deterministic, and $x$ to stand for either a univariate variable or a given point. When there are no confusions, $x$ can also be used as a generic variable for integration. The usage of several other notation and symbols is standard, such as using $ \mathbf{A}$ as either a matrix or a vector, $\|\mathbf{A}\|^2$ as the standard Euclidean norm, $a$ as a real number, and $[a]\leq a$ as the largest integer part of $a$. There are two modes of convergence involved in the rest of this paper, with “$\rightarrow_{\mathcal{D}}$ and $\rightarrow_P $ representing convergence in distribution and convergence in probability, respectively.

Linear Model Identification

Economic Interpretation and Ex Ante Specification of the Nuisance Space

The identifying content of the model comes from the specification of the nuisance space $\mathcal{S}$. The condition $m(\cdot )\in \mathcal{S}$ should be interpreted as a restriction on the mechanisms through which the endogenous component of the disturbance may depend on the observed state. An ex ante specification of $\mathcal{S}$ should satisfy three principles. First, its defining features should be motivated by economic, institutional, or design information available independently of the outcomes. Second, the restriction should have substantive content: it must exclude at least some functions of the observed state that could otherwise be confused with the structural index $\mathbf{x}^{\top }\bm{\beta }_{0}$. Third, the restriction must leave enough residual variation in $\mathbf{x}$ for ${\mathbf{Q}}_{ \mathcal{S}}=\mathbf{x}-\mathbf{P}_{\mathcal{S}}\mathbf{x}$ to be nonsingular. Ex ante specification does not require every element of the projection operator to be known before seeing the sample. The abstract space $\mathcal{S}$ may be fixed while its inner products are estimated from the observed regressors. For example, the researcher may specify in advance that $\mathcal{S}$ is generated by a set of interactions, calendar functions, or invariant functions, and then use the empirical distribution to orthonormalize those functions or estimate $\mathbf{P}_{\mathcal{S}}\, \mathbf{x}$. Because these operations use only the regressors, they implement the maintained restriction rather than select it according to its implications for the outcome. We next give several suggestive examples.

Example 1. Group and Time Effects. Suppose that observations are indexed by unit $i$ and period $t$, and consider the structural panel model $ y_{it}=\mathbf{x}_{it}^{\top }\bm{\beta}_{0}+\varepsilon _{it}.$ Institutional knowledge may suggest that the conditional mean of the disturbance consists of a unit-specific component and a common time component ${\mathbb{E}}[\varepsilon _{it}\mid x_{it}]=\alpha _{i}+\lambda _{t}.$ The nuisance space is then

This space is specified by the panel structure rather than selected from the outcome data. The residualized regressor is the familiar within-transformed variable, $\mathbf{z}_{it}=\mathbf{x}_{it}-\mathbf{P}_{\mathcal{S}}\,\mathbf{ x}_{it},$ and identification requires nondegenerate variation in $\mathbf{x} _{it}$ after removing the unit and time components: ${\mathbb{E}}[\mathbf{z} _{it}\,\mathbf{z}_{it}^{\top }]>\mathbf{0}$. The economic restriction is not simply that unit and time effects should be included because they improve fit. It is that all systematic dependence between the disturbance and the observed state operates through permanent unit heterogeneity and aggregate time shocks. If, for example, the conditional mean disturbance also contains unit-specific trends that are correlated with $\mathbf{x}_{it}$, the proposed space is misspecified. Such departures can be accommodated through the approximate restriction developed in the paper. A similar type of setting arises in high-frequency demand, production, energy, or financial data, where the institutional environment may imply that predictable omitted shocks follow a known calendar structure described by dummy variables or trigonometric functions, which can similarly be handled in this setting.

Example 2. Symmetry and Invariance Restrictions. In some cases, economic reasoning implies that the confounding component is invariant under a known transformation, whereas the structural component is not. Suppose first that $x$ is scalar and that the selection mechanism depends on the magnitude of $x$, but not on its sign: $m(x)=m(-x).$ Then one may take

The structural function $x\,\beta _{0}$ changes sign under $x\mapsto -x$, whereas every element of $\mathcal{S}_{\mathrm{even}}$ is invariant. Hence, $ \mathrm{span}\{x\}\cap \mathcal{S}_{\mathrm{even}}=\{0\},$ provided the support contains sufficiently rich positive and negative values. The substantive interpretation could be that participation, selection, or measurement error depends on the absolute magnitude of an exposure, while the structural effect depends on its direction. For example, monitoring intensity might respond to the absolute size of a position, whereas the structural return depends on whether the position is long or short. Note that the distribution of $x$ need not itself be symmetric. If it is asymmetric, $P_{\mathcal{S}_{\mathrm{even}}}\,x$ need not be zero, so the residualization is nontrivial and $E[x\,m(x)]$ may be nonzero. Identification nevertheless follows from the separation between the even nuisance functions and the odd structural index.

A multivariate version could be based on permutation invariance. Let $ x=x_{1}-x_{2},$ and suppose that $m(x_{1},x_{2})=m(x_{2},x_{1}).$ The nuisance depends only on the unordered pair of characteristics, while the structural regressor is directional and changes sign when the two components are exchanged. The relevant space is $\mathcal{S}_{\mathrm{sym}}=\left\{ s:s(x_{1},x_{2})=s(x_{2},x_{1})\right\} .$ This structure may arise in matched-pair, bilateral, ranking, or relative-performance settings in which common selection depends symmetrically on the characteristics of the pair, but the parameter of interest multiplies a directional difference.

More generally, if a group of transformations $\mathcal{G}$ acts on the observed state, one may let $\mathcal{S}$ consist of functions invariant under that group: $\mathcal{S}_{\mathcal{G}}=\left\{ s:s(gx)=s(x)\text{ for every} \ g\in \mathcal{G}\right\} .$ Identification is possible when the structural regressors contain components that transform differently from the invariant nuisance functions.

Example 3. Specified Nonlinearities and Interactions. Economic theory may imply that the omitted component operates through known nonlinear features of the observed state but does not contain an unrestricted additive linear component. For example, with two centred regressors, suppose that $m_{0}(x_{1},x_{2})=\gamma _{1}q_{1}(x)+\gamma _{2}q_{2}(x)+\gamma _{3}q_{3}(x),$ where $q_{1}(x)=x_{1}x_{2},$ $ q_{2}(x)=x_{1}^{2}-E[x_{1}^{2}],$ and $q_{3}(x)=x_{2}^{2}-E[x_{2}^{2}].$ The nuisance space is $\mathcal{S}=\mathrm{span}\{q_{1},q_{2},q_{3}\},$ whereas the structural index is $x_{1}\beta _{01}+x_{2}\beta _{02}.$ A possible economic interpretation is that the omitted state reflects complementarities, congestion, dispersion, or adjustment costs that depend on products and squared deviations, while the coefficients of interest measure additive marginal effects. In a production setting, for example, an omitted utilization component might be known to respond to the interaction of capital and labour or to deviations from normal operating levels, rather than to unrestricted linear combinations of the inputs. The basis functions may be orthogonalized with respect to the distribution of $x$. Such orthogonalization can be estimated using the regressors alone and does not change the substantive restriction.

Example 4. Design-Based Strata and Known Assignment Rules. Suppose treatment or exposure varies within strata determined by a known assignment rule. Let $A$ be a predetermined score or collection of baseline characteristics, and suppose the conditional mean disturbance is constant, or follows a specified low-dimensional form, within predefined strata: $m(A)=\sum_{j=1}^{J}\alpha _{j}\,1\{A\in \mathcal{A}_{j}\}.$ The sets $\mathcal{A}_{1},\ldots ,\mathcal{A}_{J}$ might be administrative eligibility groups, sampling strata, geographic zones, cohorts, or categories fixed by the treatment-assignment process. Then $\mathcal{S}= \mathrm{span}\left\{ 1\{A\in \mathcal{A}_{j}\}:j=1,\ldots ,J\right\} .$ The coefficient is identified by variation in $x$ within these predetermined strata: $z=x-{\mathbb{E}}[x\mid \mathcal{A}_{j}]$ for $A\in \mathcal{A}_{j}.$ The identifying assumption is that all systematic selection on unobservables is captured by the known strata.

A related case arises when an administrative assignment algorithm is known to depend on a finite set of basis functions $b_{1}(A),\ldots ,b_{K}(A)$. If economic reasoning implies that the conditional mean of the disturbance depends on the assignment variables through the same features, one may specify $\mathcal{S}=\mathrm{span}\{b_{1}(A),\ldots ,b_{K}(A)\}.$

Example 5. Low-Rank Institutional, Network, or Spatial Components. In some applications, the institutional structure suggests that confounding is generated by a small number of common components. Suppose, for example, that observations belong to known markets or network communities and that $m(x_{i})=\sum_{k=1}^{K}\gamma _{k}f_{k}(x_{i}),$ where the functions $f_{k}$ are fixed before the outcome analysis. They might consist of market indicators, predetermined network eigenvectors, geographic basis functions, or exposure to known aggregate shocks. Then $\mathcal{S}= \mathrm{span}\{f_{1},\ldots ,f_{K}\}.$ Identification comes from the component of the structural regressor that is not explained by these common institutional factors.

Motivations and heuristics

We return to the general structural linear model

Appendix (ref) shows that the proposed SP method is applicable to several classes of nonlinear models, including the general nonparametric regression: $y = g(\mathbf{x}) + \varepsilon$.

Noting that ${\mathbb{E}}[\mathbf{x}\varepsilon]\ne 0$ implies ${\mathbb{E}} [\varepsilon|\mathbf{x}]\ne 0$, we define

where $\bm{\beta}_0$ can be identifiable and consistently estimable as discussed in the rest of this section.

For the nonlinear model ((ref)), we consider the projection process of ${ \varepsilon} = {\mathbb{E}}[{\varepsilon}|\mathbf{x}] + {e}$ in ((ref)) as the proposed the SP method, which has a similar spirit to the so--called “Control Function Approach" discussed in newey1990, npv1999, np2003, JMW2015 and others, but the proposed SP method addresses endogeneity through an internally constructed IV system for such classes of correctly identified models, without assuming the existence and validity of an external IV structure.

To explain the main ideas and for notational simplicity, the main sections of this paper then focus on a semiparametric linear regression of the form:

A diagram is given in Appendix (ref).1 to illustrate model ((ref)) from a geometric point of view. Before we discuss how to identify $\bm{\beta} _0$, we make the following comments.

To address the endogeneity issue, a simple and commonly used IV method considers combining model ((ref)) with a linear decomposition of $\mathbf{ x}$ of the form:

where $\mathbf{z}_{}$ is assumed to be a valid IV. Under the invertibility of ${\mathbb{E}}\left[\mathbf{z} \, \mathbf{x}^{\top}\right]$, $\bm{\beta}_0$ can be identifiable by $\bm{\beta}_0 = \left({\mathbb{E}}\left[\mathbf{z} \, \mathbf{x}^{\top}\right]\right)^{-1} {\mathbb{E}}[\mathbf{z}_{} \, y]$.

As shown in equations ((ref))--((ref)) below, the proposed SP method consequently enables us to construct an explicit form of $\mathbf{z}$ for it to be the so--called “CIV" as a valid IV.

We now outline the main heuristics about how to construct instrumental variables by the so--called “SP method” to enable that $(\bm{\beta}_0, m(\cdot))$ is identifiable in the rest of this section.

Assume that $m(\mathbf{x})$ can be written as $m(\mathbf{x}) = \mathbf{v}( \mathbf{x})^{\top} \bm{\gamma}$, where $\bm{\gamma}=\left(\gamma_1, \cdots, \gamma_k\right)^{\top}$ is a vector of unknown coefficients, and

is a vector of series functions, where $k$ is a fixed and finite truncation parameter at this stage. The case of varying $k$ will be discussed from Section (ref).2. There is a substantial literature about series estimation, such as gallant1981, andrews1991, newey1997 , ac2003, chen2007, and BCCK2015. A recent book by dg2025 provides an update about nonparametric series estimation.

Letting $\mathbf{v}=\mathbf{v}(\mathbf{x})$ satisfy ${\mathbb{E}}[\mathbf{v} \mathbf{v}^{\top}]>0$, we rewrite model ((ref)) as

where ${\mathbb{E}}[e|\mathbf{x}]=0$. Let $w = y - \mathbf{v}^{\top} \left({ \mathbb{E}}[\mathbf{v} \mathbf{v}^{\top}]\right)^{-1} {\mathbb{E}}[\mathbf{v} y]$. Model ((ref)) can then be written as

where $\mathbf{z} = \mathbf{x} - {\mathbb{E}}[\mathbf{x} \mathbf{v}^{\top}] \left({\mathbb{E}}[\mathbf{v} \mathbf{v}^{\top}]\right)^{-1} \mathbf{v}$ satisfies ${\mathbb{E}}[e|\mathbf{z}] = {\mathbb{E}}[{\mathbb{E}}[e|\mathbf{x }]|\mathbf{z}]=0$, which, along with ${\mathbb{E}}[\mathbf{z} \mathbf{z} ^{\top}]>\mathbf{0}$, ensures that $\mathbf{z}$ can be chosen as an IV. It will be shown in Section (ref).2 that

is identifiable. We therefore have constructed a version of equation ((ref)) of the form:

which implies that $\mathbf{z}$ is the so--called CIV as a valid IV for the linear model: $y = \mathbf{x}^{\top} \bm{\beta}_0 + \varepsilon$.

Assume that $\{(\mathbf{x}_i, y_i)\}$ is a sequence of observed random variables. We introduce $\mathbf{v}_i = \mathbf{v}(\mathbf{x}_i)$, $ w_i = y_i - \mathbf{v}_i^{\top} \left(\sum_{l=1}^n \mathbf{v}_l \, \mathbf{v} _l^{\top}\right)^{-1} \, \sum_{j=1}^n y_j \,\mathbf{v}_j^{\top}$ and $ \mathbf{z}_i = \mathbf{x}_i - \sum_{j=1}^n \mathbf{x}_j \, \mathbf{v} _j^{\top}\, \left(\sum_{l=1}^n \mathbf{v}_l \, \mathbf{v}_l^{\top} \right)^{-1} \mathbf{v}_i$. A sampling version of model ((ref)) is given by:

where $e_i$ is the error term satisfying ${\mathbb{E}}[e_i|\mathbf{x}_i]=0$ and $0<{\mathbb{E}}[e_i^2|\mathbf{x}_i] = \sigma_i^2<\infty$ for $1\leq i\leq n$. Equation ((ref)) implies that $\bm{\beta}_0$ of ((ref)) can be estimated by

where $\mathbf{Z} = \left(\mathbf{z}_1, \cdots, \mathbf{z}_n\right)^{\top}$ and $\mathbf{W} = \left(w_1, \cdots, w_n\right)^{\top}$.

It is pointed out that $\widehat{\bm{\beta}}_{\mathrm{SP}}$ is the same as the following 2sLS estimator:

defined by the conventional two--stage LS method when choosing $\mathbf{Z}$ as the IV function, where $\widetilde{\mathbf{Y}} = \left(\mathbf{I}_n - \mathbf{P}_v\right) \, \mathbf{Y}$ with $\mathbf{Y} = \left(y_1, \cdots, y_n\right)^{\top}$, $\widetilde{\mathbf{X}} = \left(\mathbf{I}_n - \mathbf{P} _v\right) \, \mathbf{X}$ with $\mathbf{X} = \left(\mathbf{x}_1, \cdots, \mathbf{x}_n\right)^{\top}$, $\mathbf{P}_v = \mathbf{V} \, \left(\mathbf{V} ^{\top} \mathbf{V}\right)^{-1} \mathbf{V}^{\top}$ and $\mathbf{V} = \left( \mathbf{v}(\mathbf{x}_1), \cdots, \mathbf{v}(\mathbf{x}_n)\right)^{\top}$.

As noted in ((ref)) above, it is reasonable to require that

to ensure that $\bm{\beta}_0$ is correctly identifiable, as discussed rigorously in Section (ref).3 below.

Model identification

To present the main ideas about model identification, we introduce a Hilbert space. Let ${\mathcal{L}}^2=\{g(\mathbf{x}): \ {\mathbb{E}}[g^2( \mathbf{x})]<\infty\}$, where the inner product is given by $\langle g_1( \mathbf{x}), g_2(\mathbf{x})\rangle={\mathbb{E}}[g_1(\mathbf{x})g_2(\mathbf{x })]$ from which the induced norm is $\|g(\mathbf{x})\|=\sqrt{{\mathbb{E}} [g^2(\mathbf{x})]}$ for any given $g_1(\mathbf{x}), g_2(\mathbf{x}), g( \mathbf{x})\in {\mathcal{L}}^2$. Equipped with this $\|\cdot\|$, ${\mathcal{L }}^2$ becomes a Hilbert space.

Appendix (ref).2 provides a geometric illustration of Assumption (ref).

Given the closed subspace ${\mathcal{S}}$, we have an orthogonal decomposition for ${\mathcal{L}}^2$, that is, ${\mathcal{L}}^2={\mathcal{S}} \oplus{\mathcal{S}}^\bot $ where ${\mathcal{S}}^\bot$ is the orthogonal complement subspace of ${\mathcal{S}}$. This decomposition means that for every element $\xi \in {\mathcal{L}}^2$, it can be uniquely decomposed as $ \xi=\xi_1+\xi_2$ where $\xi_1\in{\mathcal{S}}$, $\xi_2\in {\mathcal{S}}^\bot$ , and they are orthogonal, ${\mathbb{E}}(\xi_1\xi_2)=0$. Accordingly, we have two projection mappings ${\mathcal{P}}_{\mathcal{S}}$ and ${\mathcal{M}} _{\mathcal{S}}$ that map every element in ${\mathcal{L}}^2$ into ${\mathcal{S }}$ and ${\mathcal{S}}^\bot$, respectively, defined by the decomposition, that is, ${\mathcal{P}}_{\mathcal{S}}(\xi)=\xi_1$ and ${\mathcal{M}}_{ \mathcal{S}}(\xi)=\xi_2$.

We establish the following lemma; its proof is given in Appendix C.

Lemma (ref) involves the orthonormality under a probability space. Appendix (ref).3 discusses that all mathematical moments in population can be replaced by their sampling versions, such as replacing ${\mathbb{E}} \left[\psi_j(\mathbf{x}) \, \psi_l(\mathbf{x})\right]$ by $\frac{1}{n} \sum_{i=1}^n \psi_j(\mathbf{x}_i) \, \psi_l(\mathbf{x}_i)$ for all $(j,l)$, in practice, when there is no prior knowledge about the distributional structure of the data under study.

For notational simplicity of our derivations in the rest of this section and the derivations in Appendices A, C and D, we assume that $\{\psi_j(\cdot): j\geq 0\}$ is a sequence of orthonormal functions by definition: ${\mathbb{E} }[\psi_j^2(\mathbf{x}_1)]=1$ and ${\mathbb{E}}[\psi_j(\mathbf{x}_1) \, \psi_l(\mathbf{x}_1)]=0$ for $j\neq l$. Appendices B and E indicate that, moreover, the proposed SP method works well numerically in finite--sample simulations without necessarily requiring orthonormality on $ \{\psi_j(\cdot): j\geq 1\}$.

As the closed subspace of ${\mathcal{L}}^2$, ${\mathcal{S}}$ itself, along with the norm in ${\mathcal{L}}^2$, is a Hilbert space too. After specifying the orthonormal basis, we have the expression of the projection mapping ${ \mathcal{P}}_{\mathcal{S}}$; and therefore the expression of ${\mathcal{M}}_{ \mathcal{S}}=\mathcal{I}-{\mathcal{P}}_{\mathcal{S}}$ where $\mathcal{I}$ is the identity operator. More importantly, the positive definiteness of ${\bm \Sigma}_x$ is the corner stone for the identifiability and estimation of $\bm \beta_0$.

It is noteworthy that the requirement ${\mathbb{E}}[\psi_j(\mathbf{x})]=0$ is to tailor ${\mathcal{S}}$ to cater for ${\mathbb{E}}[m(\mathbf{x})]={ \mathbb{E}}[\varepsilon]=0$. Recalling that ${\mathcal{M}}_{\mathcal{S}}( \mathbf{x})\in {\mathcal{S}}^\bot$, it is orthogonal with any element in ${ \mathcal{S}}$, in particular $m(\mathbf{x})$. Multiplying ${\mathcal{M}}_{ \mathcal{S}}(\mathbf{x})$ on both sides of equation (ref) and taking expectation imply

which yields

which is the expression of $\bm{\beta}_0$ in population.

Note that ${\mathcal{M}}_{\mathcal{S}}(\mathbf{x})$ plays an essential role in the identification of $\bm{\beta}_0$ and satisfies the necessary conditions in the construction of an IV as follows:

Lemma (ref) below shows that the expression (ref) is invariant to the choice of the basis in ${\mathcal{S}}$; its proof is given in Appendix C.

We have established a rigorous treatment about the proposed SP method for the identification of $\bm{\beta}_0$ before we discuss about how to estimate $\bm{\beta}_0$ in Section (ref) below.

In the rest of this paper, we consider the case where $d$, the dimensionality of $\mathbf{x}$, is finite and fixed. In Sections (ref)-- (ref), we focus on the case where $d$ is small. Appendix B.1 of the supplemental document discusses how to choose the series function involved for the case of $d\geq 2$ can be large but still fixed in practice.

Estimation and Inference

Estimation method and theory

Let us start to consider estimation and inference issues for model ((ref) ). Suppose that $\{\psi_j(\mathbf{x}), j\ge 1\}$ is an orthonormal basis for ${\mathcal{S}}$, where ${\mathbb{E}}[\psi_j(\mathbf{x})]=0$ and ${\mathbb{E}} [\psi_j(\mathbf{x})\psi_{\ell}(\mathbf{x})]=\delta_{j\ell}$ for all $(j,\ell) $. Hence, we have an orthogonal series expansion for $m(\mathbf{x})$ of the form:

Given a truncation parameter $k$, let $\mathbf{V}_k(\mathbf{x})=(\psi_1( \mathbf{x}), \cdots, \psi_k(\mathbf{x}))^\top$ and $\bm{\gamma}=(\gamma_1, \cdots, \gamma_k)^\top$, and then define $\delta_k(\mathbf{x} )=\sum_{j=k+1}^\infty \psi_j(\mathbf{x}) \, \gamma_j$. Model ((ref)) is written as $y=\mathbf{x}^\top {\bm \beta}_0 +\mathbf{V}_k(\mathbf{x})^\top { \bm \gamma}+\delta_k(\mathbf{x})+e$.

Moreover, given a sample $\{(y_i,\mathbf{x}_i), 1\leq i \leq n\}$, we have $ y_i=\mathbf{x}_i^\top{\bm{\beta}_0}+\mathbf{V}_k(\mathbf{x}_i)^\top {\bm \gamma}+\delta_k(\mathbf{x}_i)+e_i$ for $1\leq i\leq n$. Accordingly, we have a matrix model of the form:

where $\mathbf{y}=(y_1, \cdots, y_n)^\top$, $\mathbf{X}=(\mathbf{x}_1, \cdots, \mathbf{x}_n)^\top$, $\mathbf{V}=(\mathbf{V}_k(\mathbf{x}_1), \cdots, \mathbf{V}_k(\mathbf{x}_n))^\top$, $\mathbf{e}=(e_1, \cdots, e_n)^\top$, and $\bm{\delta}=(\delta_k(\mathbf{x}_1), \cdots, \delta_k( \mathbf{x}_n))^\top$.

Let $\mathbf{P}_v=\mathbf{V}\left(\mathbf{V}^{\top} \mathbf{V}\right)^{-1} \mathbf{V}^{\top}$ and $\mathbf{M}_v= \mathbf{I}_n - \mathbf{P}_v$. We then have

Assuming that the matrix involved is invertible, the SP estimator is then given by

Note that if ${\mathcal{S}}$ is of finite dimension, then ${\bm \delta}=0$ and hence $\widehat{\bm{\beta}}$ is unbiased. For the theoretical generalization, we treat the dimension of ${\mathcal{S}}$ as infinity in the rest of this paper.

In order to establish asymptotic consistency and normality results for $ \widehat{\bm{\beta}}$, we need to introduce the following assumption.

The proof of Theorem (ref) is given in Appendix C. As shown in part of the proof, we have $\frac{1}{n}\mathbf{X}^\top \mathbf{M}_v \mathbf{\Omega}\mathbf{M}_v\mathbf{X}\rightarrow_P \bm{\Sigma}_x \, \overline{\sigma}_e^{2}$ and $\frac{1}{n}\mathbf{X}^\top \mathbf{M}_v\mathbf{X}\rightarrow_P \bm{\Sigma}_x$ as $(k, n)\to (\infty, \infty)$. Thus, the asymptotic covariance matrix is $\bm{\Sigma}_x^{-1} \, \overline{\sigma}_e^{2}$.

To start our discussion on the weak endogeneity setting, we propose a simple statistic to test for a full level of exogeneity versus a wide range of weak endogeneity in Section (ref).2.

Testing for weak endogeneity

Consider $y =\mathbf{x}^{\top}\bm{\beta}_0 + m(\mathbf{x}) + e$ with $m( \mathbf{x})={\mathbb{E}}[\varepsilon|\mathbf{x}]$ under the following null hypothesis:

We approximate $m(\mathbf{x})$ by $\sum_{j=1}^k \psi_j(\mathbf{x}) \gamma_{j} \equiv \mathbf{V}_k(\mathbf{x})^{\top} \bm{\gamma}$, where $ \mathbf{V}_k(\cdot) = \left(\psi_1(\cdot), \cdots, \psi_k(\cdot)\right)^{\top}$, and ${\bm\gamma} =\left(\gamma_{1}, \cdots, \gamma_{k}\right)^{\top}$ can be estimated by $\widehat{\bm\gamma}$ in the same way as in Section (ref).1.

In order to test $H_0: \, {\mathbb{P}}(m(\mathbf{x}) =0)=1$, it suffices to test ${\bm\gamma}=0$, or equivalently $(\mathbf{V}^{\top} \mathbf{V}) \, {\bm \gamma}=0$. Since the ${\bm\gamma}$ parameter can be estimated by $\widehat{ \bm\gamma} = (\mathbf{V^{\top} V})^{-1} \mathbf{V}^{\top} (\mathbf{y} - \mathbf{X}\widehat{\bm\beta})$, we introduce $\widetilde{\bm\gamma}_n = ( \mathbf{V^{\top} V}) \, \widehat{\bm\gamma}$ and then define for each given point $x$:

This motives us to propose a simple test statistic (a nonparametric version of the Hausman test proposed in hausman1978) of the form:

by orthonormality: $\int_{\mathcal{X}} \mathbf{V}_k({x}) \mathbf{V}_k^{\top}( {x}) \, f_{\mathbf{x}}(x) d {x}=\mathbf{I}_k$, where $f_{\mathbf{x}}(\cdot)$ denotes the density of $\mathbf{x}$.

In order to avoid dealing with robustness issues regarding the choice of individual $k$ values in practice, we propose using a summarized version of $ L_n(k)$ of the form:

where $k_{\min}$ and $k_{\max}$ are the respective smallest and largest integers satisfying $1\leq k_{\min} < k_{\max}<[\sqrt{n}]$, and $k_{\max} - k_{\min}\rightarrow \infty$ and $\frac{k_{\max}^2}{n}\rightarrow 0$ as $ n\rightarrow\infty$.

Note that $\frac{S_n^2}{\sigma_n^2}\rightarrow_P 1$ with $\sigma_n^2 = \frac{ 2 \,{\sigma}_e^4}{3} \cdot n^2 \, k_{\max}^3\, (1+ o(1))$ as shown in the proof of Theorem (ref) in Appendix C. To show the consistency of our testing statistic under a sequence of local alternatives, we test

where positive sequence $a_n\to 0$ with certain rate while $m(\mathbf{x})\in {\mathcal{L}}^2$ and ${\mathbb{E}}[m^2(\mathbf{x})]>0$.

It is pointed out that such a sequence of local alternatives naturally covers a wide range of weak endogeneity. The following theorem establishes the consistency of the proposed test under $H_1$; its proof is outlined in Appendix C.

Note that Theorem (ref) shows that the proposed test is capable to detect certain types of weak endogeneity of an order of $a_{n}$, as long as $ a_n$ satisfies $a_{n}^2 \, n \, k_{\max}^{-1/2} \rightarrow \infty$.

Note also that the fast possible rate of $a_n = n^{-\frac{1}{2}} \, k_{\max}^{c_1}$ is near--optimal when $k_{\max} \rightarrow \infty$ at the slowest possible rate, such as $k_{\max} = \left[ c_2 \, \log(\log(n))\right] $, for $c_1>\frac{1}{4}$ and $c_2>0$ to be chosen by the user, although the conventional optimal rate of $a_n = n^{-1/2}$ (see, for example, ss1997) chosen in the parametric setting is not achievable under the proposed SP framework.

In Appendix (ref) below, we show how to extend the proposed SP method to a class of nonlinear models. Appendix A.3 of the supplementary document outlines the main ideas about how to extend the SP method to deal with some other classes of nonlinear models, including two classes of non-- and semi--parametric regression models, and one class of binary models associated with nonlinearity and endogeneity.

Linear Models under Weak Endogeneity

The existing literature pays particular attention on the case where $ (\varepsilon, \mathbf{x})$ in ((ref)) follows a joint Gaussian distribution, which implies that $m(\mathbf{x}) = \left(\mathbf{x} - { \mathbb{E}}[\mathbf{x}]\right)^{\top} \bm{\gamma}_0$ and model ((ref)) becomes $y =\mathbf{x}^{\top} \bm{\beta}_0 + \left(\mathbf{x} - {\mathbb{E}}[ \mathbf{x}]\right)^{\top} \bm{\gamma}_0+ e$ with ${\mathbb{E}}[e|\mathbf{x} ]=0$. Appendix A.4.3 of the main supplementary document discusses that the SP method is applicable to consistently estimate $\bm{\beta}_0$ in certain linearity cases. We now focus on several weak endogeneity settings. When $m( \mathbf{x}) = \left(\mathbf{x} - {\mathbb{E}}[\mathbf{x}]\right)^{\top} \bm{\gamma}_0$, we rewrite model ((ref)): $y = \mathbf{x}^{\top} \bm{\beta}_0 + \varepsilon$ as

which implies that $\left(\bm{\beta}_0 + \bm{\gamma}_0\right) = {\mathbb{E}} ^{-1}\left[\left(\mathbf{x} - {\mathbb{E}}[\mathbf{x}]\right) \left(\mathbf{x } - {\mathbb{E}}[\mathbf{x}]\right)^{\top}\right] \, {\mathbb{E}}\left[\left( \mathbf{x} - {\mathbb{E}}[\mathbf{x}]\right)\, (y - {\mathbb{E}}[y])\right]$ can be correctly identified collectively, rather than $\bm{\beta}_0$ individually.

When there are some additional restrictions imposed on $\bm{\gamma}_0$, such as those weak endogeneity scenarios discussed below, however, $\bm{\beta}_0$ may still be asymptotically identifiable. In such cases where $\bm{\gamma}_0 \equiv \bm{\gamma}_{n0} \rightarrow \mathbf{0}$ as $n\rightarrow \infty$, $ \bm{\beta}_0$ can then be estimated consistently.

We start with model ((ref)) associated with a type of weak endogeneity of the form:

where ${\mathbb{E}}[e_i|\mathbf{x}_i]=0$ for $1\leq i\leq n$, $m_n(\mathbf{x} )={\mathbb{E}}[\varepsilon_i|\mathbf{x}_i = \mathbf{x}] = a_n \, m(\mathbf{x} )$ with ${\mathbb{E}}[m(\mathbf{x})]=0$, ${\mathbb{E}}[\mathbf{x} \, m( \mathbf{x})] \neq 0$, ${\mathbb{E}}[m^2(\mathbf{x})]>0$, $m(\mathbf{x}) \in { \mathcal{S}}$, $a_n\rightarrow 0$ and $a_n \, n\, k^{-\frac{2s}{d} }\rightarrow 0$ as $n\rightarrow \infty$, in which $s$ is the smoothness order of $m(\cdot)$ as in Assumption (ref)(iv). We have as $ n\rightarrow \infty$

Probably because of the projection in ((ref)) and then the definition of weak endogeneity in ((ref)), we are able to show that $\bm{\beta}_0$ can still be consistently estimated.

Since we assume that $m(\mathbf{x}) \in {\mathcal{S}}$, we approximate $m( \mathbf{x})$ in the same way as in Section (ref). In view of model ((ref)), the `endogenous component' represented by $\mathbf{v}^{\top} \bm{\gamma}_n$, with $\bm{\gamma}_n= a_n \, \bm{\gamma}$, has been eliminated. Consequently, equations ((ref)) and ((ref)) show that $ \bm{\beta}_0$ of ((ref)) is correctly identified and therefore it can still be estimated by $\widehat{\bm{\beta}} =\widehat{\bm{\beta}}_{ \mathrm{SP}}$ of ((ref)) consistently.

The establishment and the proof of Theorem (ref) imply Theorem (ref) below.

Theorem (ref), along with Theorem (ref), shows that the proposed SP method is therefore invariant to the degree of endogeneity under study.

Let us return to the linear weak endogeneity case where $m_n(\mathbf{x}) = \frac{\bm{\mu}^{\top} \left(\mathbf{x} - {\mathbb{E}}[\mathbf{x}]\right)}{ \sqrt{n}}$ with $\|\bm{\mu}\|^2<\infty$. There is a substantial literature as reviewed in a recent survey by ass2019 about treatments of weak endogeneity involving weak IVs.

Let $\bm{\beta}_{n, \bm{\mu}} \equiv: \bm{\beta}_0 + s_n^{-1} \, \sigma_n \, \frac{\bm{\mu}}{\sqrt{n}}$. It follows that the standard OLS estimator becomes

where $m_n(\mathbf{x}) =\frac{\bm{\mu}^{\top} \left(\mathbf{x} - {\mathbb{E}} [\mathbf{x}]\right)}{\sqrt{n}}$, $s_n = \sum_{i=1}^n \mathbf{x}_i \, \mathbf{ x}_i^{\top}$ and $\sigma_n = \sum_{i=1}^n \mathbf{x}_i \, \left(\mathbf{x}_i - {\mathbb{E}}[\mathbf{x}_i]\right)^{\top}$.

We then establish the following theorem; its proof follows trivially from verifying Lindeberg conditions needed for a central limit theorem for an i.i.d. data setting.

Equation ((ref)) implies that $\widehat{\bm{\beta}}_{\mathrm{LS}} = \bm{\beta}_0 + \widehat{\bm{\Sigma}}_n^{-1} (\mathbf{x})\, {\mathbb{E}}^{-1} [\mathbf{x} \, \mathbf{x}^{\top} ] \, \mathrm{var}[\mathbf{x}] \, \frac{ \bm{\mu}}{\sqrt{n}} + O_P\left(\frac{\widehat{\bm{\Sigma}}_n^{-1}(\mathbf{x}) }{\sqrt{n}}\right) \rightarrow \bm{\beta}_0$, indicating that $\widehat{ \bm{\beta}}_{\mathrm{LS}}$ is a consistent estimator of $\bm{\beta}_0$.

However, $\bm{\mu}$ itself cannot be identified correctly. This becomes clearer in the special case where $\bm{\beta}_0= \mathbf{0}$, $ m^{\ast}(\cdot) =0$ and $y = \left(\frac{\bm{\mu}}{\sqrt{n}}\right)^{\top} \left(\mathbf{x} - {\mathbb{E}}[\mathbf{x}]\right) + e$. In this case, we have as $n\rightarrow \infty$

which shows that $\widehat{\bm{\mu}} = \sqrt{n} \, \left(\sum_{i=1}^n \left( \mathbf{x}_i - \overline{\mathbf{x}}_n\right) \, \left(\mathbf{x}_i - \overline{\mathbf{x}}_n\right)^{\top}\right)^{-1} \, \sum_{i=1}^n \left( \mathbf{x}_i - \overline{\mathbf{x}}_n\right) \, y_i$ is not consistent to $ \bm{\mu}$, although $\frac{\widehat{\bm{\mu}}}{\sqrt{n}}$ is consistent to $ \frac{\bm{\mu}}{\sqrt{n}}$, where $\overline{\mathbf{x}}_n = \frac{1}{n} \sum_{i=1}^n \mathbf{x}_i$.

This is similar in spirit to those discussed in the weak IV literature, see equation (2.5) of ss1997, for example. We demonstrate in Appendix B.4 of the supplementary document that one may eliminate $\bm{\mu}$ by a simple jackknife method.

Furthermore, we show that the SP method offers a simple way to deal with a mixture of strong and weak endogeneity of the form:

which implies

where ${\mathbb{E}}[m^{\ast}(\mathbf{x})]=0$, ${\mathbb{E}}|m^{\ast}(\mathbf{ x})|^2>0$ and $m^{\ast}(\mathbf{x}) \in {\mathcal{S}}$. Our approach allows either ${\mathbb{E}}[\mathbf{x} \, m^{\ast}(\mathbf{x})] = 0$ or ${\mathbb{E} }[\mathbf{x} \, m^{\ast}(\mathbf{x})] \neq 0$.

We use the same $\{\psi_j(\mathbf{x}): j\geq 1\}$ as an orthonormal series to approximate $m^{\ast}(\mathbf{x})$ by $m_{k}^{\ast}(\mathbf{x}) \equiv:\sum_{j=1}^{k} \psi_j(\mathbf{x}) \, \gamma_{j}^{\ast}$, where $ \bm{\gamma}^{\ast}_k = \left(\gamma_{1}^{\ast}, \cdots, \gamma_{k}^{\ast}\right)^{\top}$ is an array of coefficients.

In the same spirit as in Section (ref).1, we estimate $\bm{\beta}_{n, \bm{\mu}}= \bm{\beta}_0 + s_n^{-1} \, \sigma_n \, \frac{\bm{\mu}}{\sqrt{n}}$ by $\widehat{\bm{\beta}}_{\mathrm{SP}}$ of ((ref)). Obviously, $\bm{\beta} _0$ is only asymptotically identifiable. We finally establish the following theorem; its proof follows trivially from that of Theorem (ref).

Equation ((ref)) implies $\widehat{\bm{\beta}}_{\mathrm{SP}} = \bm{\beta}_0 + \widehat{\bm{\Sigma}}_n^{-1}(\mathbf{x}, k) \, {\mathbb{E}} ^{-1} [\mathbf{x} \, \mathbf{x}^{\top} ] \, \mathrm{var}[\mathbf{x}] \, \frac{\bm{\mu}}{\sqrt{n}} + O_P\left(\frac{\widehat{\bm{\Sigma}}_n^{-1}( \mathbf{x}, k)}{\sqrt{n}}\right) \rightarrow \bm{\beta}_0$, which indicates that $\widehat{\bm{\beta}}_{\mathrm{SP}}$ is still a consistent estimator of $\bm{\beta}_0$.

Section (ref) proposes a simple LASSO selection method to realise the proposed SP method in practice before we employ it for a real dataset analysis in Section (ref) and then in Appendix B.2 of the supplementary document for the finite--sample evaluation of simulated datasets and another real dataset.

LASSO Selection and Empirical Analysis

Construction of $\mathcal{S}$

Constructing the subspace ${\mathcal{S}}$ involved in Assumption 2.1 is crucial in our theory that affects both the asymptotic consistency and unbiasedness of the estimator proposed. We focus on the linear model discussed in Sections 2 and 3, and a similar discussion for the nonlinear setting in Appendix A.2 follows accordingly.

As can be seen, since the limit covariance matrix of $\widehat{\bm \beta}_{ \mathrm{SP}}$ of (3.2) is proportional to ${\bm \Sigma}_x^{-1}$ where ${\bm \Sigma}_x={\mathbb{E}}[\mathbf{x}\mathbf{x}^\top]-\sum_{j=1}^\infty {\mathbb{ E}}[\mathbf{x}\psi_j(\mathbf{x})]{\mathbb{E}}[\psi_j(\mathbf{x})\mathbf{x} ^\top]$, in which the smaller the number of basis functions we use, the more efficient $\widehat{\bm \beta}_{\mathrm{SP}}$ can be. Therefore, it is not unreasonable that we choose such $\mathcal{S}$ that forms the smallest subspace containing $m(\mathbf{x})$. We then propose a data driven approach to selecting the basis functions that spans the subspace $\mathcal{S}$.

Note that, as a closed subspace of the Hilbert space $\mathcal{L}^2$, $ \mathcal{S}$ is a Hilbert space too. Since $\mathcal{L}^2$ is separable, it possesses an orthonormal sequence as its basis (see, Page 169 of dudley2003, for example); some subsequence of this basis will be the basis of $\mathcal{S}$ by construction. Precisely, let $\{\phi_j(\mathbf{x}), j\ge 0\}$ be the basis of $\mathcal{L}^2$ and the subsequence $\{\phi_{s_j}(\mathbf{x}), j=1, 2, \cdots\}$ be the basis of $\mathcal{S}$ that we denote as $\{\psi_j(\mathbf{x }), j=1, 2, \cdots\}$ in previous sections, that is, $\psi_j(\mathbf{x} )=\phi_{s_j}(\mathbf{x})$ for $j=1,2 \cdots$. The choice of $\{\phi_j( \mathbf{x}), j\ge 0\}$ should avoid orthogonal polynomial sequence because the endogeneity implies ${\mathbb{E}}[\mathbf{x}m(\mathbf{x})]\ne 0$; otherwise, this condition cannot be fulfilled. Consequently, with the help of Lasso approach and the regression equation, along with Assumption (ref), we are able to select $\{\psi_j(\mathbf{x}), j=1, 2, \cdots\}$ from its parent sequence $\{\phi_j(\mathbf{x}), j\ge 0\}$, and therefore, $ \mathcal{S}=$span$\{\psi_j(\mathbf{x}), j=1, 2, \cdots\}$.

{To elaborate this data driven method,} without loss of generality, we consider the function space ${\mathcal{L}}^2$ where the support of $\mathbf{x }$ is $[0,\pi]$. Accordingly, we introduce an orthonormal basis in ${ \mathcal{L}}^2$ of the form: $\{\phi_0(\mathbf{x}) =1, \phi_j(\mathbf{x})= \sqrt{2/\pi}\cos( j \mathbf{x}) \text{ for } j\ge 1\}$, which has been employed for the numerical exercises conducted in Section (ref) and Appendix B.2. In Appendix E of the online supplementary document, we also consider some other basis functions to check the robustness of the SP method. In the rest of this subsection, our goal becomes to find the basis of the subspace $\mathcal{S}$ from the basis $\{\phi_j(\mathbf{x}), j\ge 0\}$ in ${\mathcal{L}}^2$.

{Since $m(\mathbf{x})\in \mathcal{S}\subsetneqq{{\mathcal{L}}^2}$,} we may expand $m(\mathbf{x})$ in terms of $\{\phi_j(\mathbf{x}), j\ge 0\}$ by $m( \mathbf{x})=\sum_{j=1}^\infty \phi_j(\mathbf{x}) \, \xi_j$, where we use $ \{\xi_j: j\geq 1\}$ as the coefficients {in order to distinguish the $ \{\gamma_j: j\geq 1\}$ in the expansion of $m(\mathbf{x})$ in terms of $ \{\psi_j(\mathbf{x}), j\ge 1\}$. Here, we abandon the constant term $\phi_0( \mathbf{x})$ due to ${\mathbb{E}}[m(\mathbf{x})]=0$. More importantly, many $ \gamma_j$ are either zero or statistically insignificant. This is because $m( \mathbf{x})$ belongs to the subspace $\mathcal{S}$, and by Parseval's identity in Hilbert space, ${\mathbb{E}}[m(\mathbf{x})^2]=\sum_{j=1}^\infty \xi_j^2<\infty$ implying the attenuation of $\xi_j$.}

{Meanwhile, from the sieve estimation point of view, we could not estimate an infinite dimensional function; chen2007 (see, Page 5561 for example). Instead, we use a “sieve”, $\cdots \subset {\mathcal{L}}^2_k\subset {\mathcal{ L}}^2_{k+1}\subset \cdots \subset {\mathcal{L}}^2$ where ${\mathcal{L}}^2_k=$ span$\{\phi_j(\mathbf{x}), j=1, \cdots, k\}$ that approach ${\mathcal{L}}^2$ ; correspondingly, ${\mathcal{S}}_{(k)}=$span $\{\phi_{s_j}(\mathbf{x}), s_j\le k\}$ form a sieve for ${\mathcal{S}}$. Once we identify each element in ${\mathcal{S}}_{(k)}$ with high probability for each $k$, we obtain $ \widehat{\mathcal{S}}_{(k)}$ by the data driven method.}

{To proceed,} given a sample $\{(\mathbf{x}_i, y_i), i=1,\cdots,n\}$, we can rewrite the model as follows:

where $k$ is the truncation parameter, and $\delta_{k}(\mathbf{x} _i)=\sum_{j=k+1}^\infty \phi_j(\mathbf{x}_i) \, \xi_j$. The parameters $ \{\xi_j\}$ need to meet certain sparsity conditions that will be specified soon.

Based on (ref), we propose the following LASSO selection method.

where $\widehat{\bm\xi}^*=(\widehat{\xi}_1^*,\ldots, \widehat{\xi}_k^*)^\top$ , and $k$ and $\lambda$ are the user--chosen truncation and tuning parameters, respectively.

where $\{\zeta_j\}$ is a set of pre-determined weights such as $1/\widehat{ \xi}_j^*$ with $\widehat{\xi}_j^*$ being the $j^{th}$ element of $\widehat{ \bm\xi}^*$. \ The LASSO estimate $\widehat{\bm\xi}^\dag=(\widehat{ \xi} _{1}^\dag,\ldots, \widehat{ \xi}_{k}^\dag)^\top$ enables us to obtain

for which we show that $\Pr(\widehat{\mathcal{S}}_{(k)}=\mathcal{S} _{(k)})\rightarrow 1$ as $(n, k) \rightarrow (\infty, \infty)$ in the following Theorem (ref).

Assumption (ref)(i) is a typical sparse restriction and infers that $\mathscr{C}\cap \bar{\mathscr{C}}=\emptyset$ and $\sharp \bar{ \mathscr{C}}=k-k_0$. Assumption (ref)(ii) regulates the truncation parameter and the sample size further. {In some special cases, e.g. $m(\cdot) $ is even function or odd function, $k_0\le k/2$, this condition is fulfilled.}

We establish the following asymptotic consistency; its proof is given in Appendix C.

The additional conditions in Theorem (ref) can be easily fulfilled. It is worth mentioning that $\Pr(\operatorname*{\normalfont\textrm{sgn}}(\widehat{\pmb{\xi}} _k^{\dag}) =\operatorname*{\normalfont\textrm{sgn}}(\pmb{\xi})) \to 1$ implies that $\Pr(\widehat{\mathcal{S}} _{(k)}=\mathcal{S}_{(k)})\rightarrow 1$ as $(k, n)\to (\infty, \infty)$. Consequently, we have chosen the subset $\mathcal{S}_{(k)}$ asymptotically.

The return to schooling case study

\noindentExample 5.1: In this example, we examine the returns to schooling (see, card2001, for example) using data from the 1979 National Longitudinal Survey of Youth\footnote{ The data is collected from \url{https://davidcard.berkeley.edu/data_sets.html}.}. After removing those individuals having missing values, we have $n=2639$ observations left in the dataset. We specifically consider the following variables:

{\

}

In Table (ref), $y$ is the $\log$(wage) (i.e., the response), and $x$ is the endogenous regressor (i.e., the number of years of child's education received). We follow Chapter 15 of wooldridge2016 as well as agmr2023 and cff2025 to consider the following set of IVs: $z_1$ stands for the number of years of dad's education; $z_2$ equals to 1 if the individual grew up near an accredited four-year college, 0 otherwise; $z_3$ equals to 1 if the individual grew up near an accredited two-year college, 0 otherwise; $z_4 $ equals to 1 if the individual grew up near an accredited four-year public college, 0 otherwise; $z_5$ equals to 1 if the individual grew up near an accredited four-year private college, 0 otherwise.

We consider the following linear model:

and our focus is to infer $\beta_0$. To be consistent with the simulation designs considered in Appendix B, we rescale $x_i$'s as follows:

where $\underline{x} =\min_{i}x_i$ and $\overline{x}=\max_i x_i$. Similarly, we normalize $z_{1i}$'s. This is in the same spirit of the data transformation process as in cff2025. Thus, $\widetilde{x}_i\in [0,\pi]$. After the data transformation (ref), the model becomes

where $\widetilde{\beta}_0 = \frac{1}{\pi}(\overline{x} - \underline{x} )\beta_0$, $\widetilde{\alpha}_0 =\frac{1}{\pi} \underline{x}\beta_0 + \alpha_0$, and $\varepsilon_i =\widetilde{m}(\widetilde{x}_i)+e_i$ with $ \widetilde{m}(\widetilde{x}_i) = m(\widetilde{x}_i \cdot \frac{\overline{x} - \underline{x}}{\pi}+ \underline{x}) =m(x_i)$. By doing so, we maintain the property of $\varepsilon_i$ in order to compare with the 2SLS method widely adopted in the literature. For the purpose of comparison, we consider the following methods:

(i) SP method; \, (ii) LS method: (a) regressing $y$ on $(1, \widetilde{x})$ ; (b) regressing $y$ on $(1, \widetilde{x}, z_{1}, z_{1}^2)$; \, (iii) 2sLS method for the model (ref) under two scenarios: (a) Just identified: using $z_j$ as IV only for $j=1,2,\ldots, 5$ respectively; (b) Over identified: using two IVs $z_{1}$ and $z_{j}$, where $j=2,3,4,5$; (c) Over identified: using $z_{1}$ and $\phi_2(\cdot),\phi_3(\cdot),\phi_8(\cdot)$ together as IVs, where $\phi_2(\cdot)$, $\phi_3(\cdot)$, $\phi_8(\cdot)$ are selected by the SP method; and (d) Over identified: using $\phi_2(\cdot)$, $ \phi_3(\cdot)$, $\phi_8(\cdot)$ as IVs.

In what follows, we consider the case where the level of the dad's education is included in the following model:

where $x$ and $z_1$ denote the number of years of the education level received for the child and dad, respectively. We estimate the model using the SP method as well. We then select $\mathcal{S}$ from $\{\phi_0( \widetilde{x}),\phi_1(\widetilde{x}),\ldots \}\otimes \{\phi_0(z),\phi_1(z),\ldots \}$. Using the proposed LASSO method, we select $\phi_1(\widetilde{x})\phi_6(z)$, $\phi_2(\widetilde{x})$, and $\phi_3( \widetilde{x})\phi_1(z)$. The selection of the basis functions is somewhat consistent with those in the main text, as we also identify $\phi_2( \widetilde{x})$ in the LASSO procedure. The estimated $m(x,z)$ is $\widehat{m }(x,z)= 0.0680\phi_1(\widetilde{x})\phi_6(z)-0.0678 \phi_2(\widetilde{x}) -0.1232 \phi_3(\widetilde{x})\phi_1(z)$. We refer to this model as SP* in what follows, and still focus on reporting $\widetilde{\beta}_0$.

For each method, we report the estimated value of $\widetilde{\beta}_0$ and it 95% confidence interval (referred to CI in Table (ref)). Moreover, we report the in-sample root mean squared errors of each method: $ \text{RMSE} =\sqrt{\frac{1}{n}\sum_{i=1}^n (y_i-\widehat{y}_i)^2}$.

{

}

We summarize the results in Table (ref). All methods confirm a positive relationship between wage and education, and the SP method offers better model fitting in terms of RMSE. The LASSO method selects $ \{\phi_2(x), \phi_3(x), \phi_8(x)\}$. It is noteworthy that Table (ref) provides the estimates of both $\widetilde{\beta}_0$ and $\beta_0$ . After accounting for the relationship $\widetilde{\beta}_0 = \frac{1}{\pi}( \overline{x} - \underline{x})\beta_0$, the estimates of $\beta_0$ are reported in the last column. The result from the SP method shows that an additional year eduction suggests around 7.1% increase in wage on average. The finding from the SP method is very close to 2sLS method when the IVs are selected as $z_1$ only, as $(z_1, z_3)$, or as $(z_1, z_5)$.

Meanwhile, Table (ref) also shows that the proposed SP* model accommodates the case where $\mathbf{x}$ includes externally observed variables, such as the level of the dad's education. Among the SP, SP* and LS(b) estimates, SP* performs better than SP and LS(b), while LS(b) is a commonly used one in such empirical analysis. By using different IVs, obviously, the outcomes of 2sLS vary from case to case, which does create some uncertainty for practical interpretation. By contrast, the SP and SP* estimation offers robust estimates regardless of the choice of the series functions and truncation parameters.

Appendix B.3.2 examines possible source of endogeneity involved in Example 5.1. In Appendix B.3.3, additionally, we discuss possible endogeneity for a stock return dataset and offer a simple and robust way to deal with the case where the structural model error term may be endogenously correlated with the market portfolio variable under study.

Conclusions and Discussion

This paper has proposed the SP method to deal with certain types of endogeneity involved in linear and nonlinear models. The corresponding estimation and testing theory has been established, and evaluated by both simulated and real datasets. Our finite--sample results established in Section (ref).2, Appendix B of the supplementary document and Appendix E of the online appendix have shown that the proposed SP method works well numerically with the support of the sound large--sample theory established in Sections (ref)--(ref).1.

Appendix (ref) below discusses about how to extend the proposed SP method to the identification, estimation and testing for a class of nonlinear regression models. As alluded in the introduction and then outlined in Appendix A.3 of the supplementary document, the proposed SP method can also be extended in a unified way to a number of non-- and semi--parametric regression models associated endogeneity issues.

Since such models have their own identifiability issues, existing developments in the relevant literature, and there are estimation and testing properties to be established for each of them, we wish to write up details for each of them separately.

{

Linear and Nonlinear Models

Linear regression

A diagrammatic illustration about model ((ref))

In model ((ref)): $y={\bf x}^\top{\bm\beta}_0+\varepsilon$, the regressor ${\bf x}$ has an effect on $y$ in two ways. One is a direct effect that adds ${\bf x}^\top {\bm \beta}_0$ to $y$; and the other is an indirect effect that affects $y$ through $\varepsilon$, that is, when ${\bf x}$ changes, the error term $\varepsilon$ changes as well due to correlation with ${\bf x}$, as illustrated below.

More precisely, the projection $m({\bf x})={\mathbb E}(\varepsilon|{\bf x})$ splits the contribution of $\varepsilon$ to $y$ into two parts, $\varepsilon=m({\bf x})+e$, where the first contribution $m({\bf x})$ stems from ${\bf x}$ while the second contribution $e$ is due to all factors orthogonal with ${\bf x}$. The effect of $m({\bf x})$ on $y$ is of all factors implicitly included in $\varepsilon$ that are correlated with ${\bf x}$. When ${\bf x}\sim {\bf x}+\Delta {\bf x}$, we have

For the endogenous model, because of these two--fold effect of ${\bf x}$ on $y$, OLS could not identify the direct effect of ${\bf x}$. Our method extracts the “energy” of ${\bf x}$ in $\varepsilon$, exposing the effect $m({\bf x})$ in the equation, so that under certain conditions we can identify the direct and indirect effects of ${\bf x}$ on $y$.

A geometric illustration about Assumption (ref)

We have defined the Hilbert space ${\mathcal L}^2=\{g({\bf x}):\; {\mathbb E}[g^2({\bf x})]<\infty\}$, in which ${\mathbb E}[\xi_1\xi_2]$ is an inner product for $\xi_1, \xi_2\in {\mathcal L}^2$; in Assumption (ref), we stipulate ${\mathcal S}$ as the closed proper subspace of ${\mathcal L}^2$. Accordingly, we have direct sum decomposition ${\mathcal L}^2={\mathcal S}\oplus {\mathcal S}^\bot$ where ${\mathcal S}^\bot$ is the orthogonal complement subspace of ${\mathcal S}$ that is the set of all elements from ${\mathcal L}^2$ orthogonal with each element in ${\mathcal S}$. Notice that both ${\mathcal S}$ and ${\mathcal S}^\bot$ are Hilbert space too, with the norm defined in ${\mathcal L}^2$.

At element level, for any $\xi\in{\mathcal L}^2$, one can uniquely decompose $\xi=\xi_1+\xi_2$ with $\xi_1\in{\mathcal S}$ and $\xi_2\in{\mathcal S}^\bot$. Accordingly, two projection mappings can be defined as ${\mathcal P_S}(\xi)=\xi_1$ and ${\mathcal M_S}(\xi)=\xi_2$. Certainly, ${\mathcal M_S}\equiv {\mathcal I-P_S}$ where ${\mathcal I}$ is the identity operator.

In Assumption (ref), the condition $\lambda^\top {\bf x}\not \in {\mathcal S}$ for any $\lambda\neq 0$ implies that ${\mathcal M_S}(\lambda^\top{\bf x})$ does not degenerate, as illustrated in Figure (ref).

{ \tikzset{global scale/.style={ scale=#1, every node/.append style={scale=#1} } } {

} }

Consequently, as shown in Lemma 2.1 above, we have the following identifiability condition:

ensures that $\bm{\beta}_0 = {\bm\Sigma}_x^{-1} \, \bm{\Sigma}_{xy}$ is identifiable, where $\bm{\Sigma}_{xy} = {\mathbb E}[\mathbf{x} \, y] - \sum_{j=1}^\infty {\mathbb E}[{\bf x} \, \psi_j({\bf x})] \, {\mathbb E}[y \, \psi_j({\bf x})]$. Furthermore, Lemma 2.2 shows that $\bm{\beta}_0$ is invariant to the choice of the basis $\{\psi_j(\cdot): j\geq 1\}$.

Probability density function of $\mathbf{x}$ unknown

Since the space ${\mathcal L}^2$ is defined as the set of all $g({\bf x})$ such that ${\mathbb E}[g^2({\bf x})]<\infty$, when the density of $\mathbf{x}$ is unknown, it is not easy to determine which function of ${\bf x}$ belongs to the space. In this case, we may relax the restriction ${\mathbb E}[g^2({\bf x})]<\infty$.

Let $\pi(x)$ be a density satisfying $c\pi(x)\le f_{\bf x}(x)\le C\pi(x)$ for some constants $C>c>0$, where $f_{\bf x}(x)$ is the density of ${\bf x}$. Define ${\mathcal L}^2_{\pi}=\{h({\bf x}): \int h^2(x)\pi(x)dx<\infty\}$. It is clear that $g({\bf x})\in {\mathcal L}^2$ iff $g({\bf x})\in {\mathcal L}^2_{\pi}$. Let ${\mathcal S}_{\pi}$ be a proper closed subspace of ${\mathcal L}^2_{\pi}$, and ${\mathcal S}_{\pi}={\rm span}\{\xi_j({\bf x}), j\ge 1\}$ where $\{\xi_j(x)\}$ is an orthonormal sequence w.r.t. $\pi(x)$, that is, $\int \xi_i(x)\xi_j(x)\pi(x)dx=\delta_{ij}$.

Because of Proposition 2.1 in BCCK2015, for any vector $U_k({\bf x})=(\xi_1({\bf x}), \cdots, \xi_k({\bf x}))$, the eigenvalues of the matrix ${\mathbb E}[U_k({\bf x})U_k({\bf x})^\top]$ are bounded above and away from zero, uniformly over $k$. This implies that any subset of $\{\xi_j({\bf x}), j\ge 1\} $ is linearly independent, so we use Gram-Schmidt procedure to orthogonalize $\{\xi_j({\bf x}), j\ge 1\} $ to be $\{{\psi}_j({\bf x}), j\ge 1\}$, such that (1) ${{\mathbb E}}[{\psi}_i({\bf x}){\psi}_j({\bf x})]=\delta_{ij}$, and ${\mathbb E}[{\psi}_j({\bf x})]=0$ for all $i, j$; (2) ${\rm span}\{{\psi}_j({\bf x}), j\ge 1\}={\rm span}\{{\xi}_j({\bf x}), j\ge 1\}$.

Therefore, $\{{\psi}_j({\bf x}), j\ge 1\}$ is an orthonormal sequence under inner product $\langle\xi, \eta\rangle={\mathbb E}(\xi \eta)$, and ${\mathcal S}={\rm span}\{{\psi}_j({\bf x}), j\ge 1\}$ is a proper subspace of ${\mathcal L}^2$; moreover, $m({\bf x})\in {\mathcal S}$ iff $m({\bf x})\in {\mathcal S}_\pi$.

We shall show these assertions using induction. First, let ${\psi}_1({\bf x})=\frac{\displaystyle1}{\displaystyle{\rm s.d.}(\xi_1({\bf x})) }\{\xi_1({\bf x})-{\mathbb E}[\xi_1({\bf x})]\}.$ Clearly, the assertion holds for $k=1$, that is, ${\rm span}\{\psi_1({\bf x})\}={\rm span}\{{\xi}_1({\bf x})\}$ and ${\mathbb E}[{\psi}_1({\bf x})]=0$.

Suppose the assertion holds for $k-1$, i.e. $\{{\psi}_j({\bf x}), j=1,\cdots, k-1\}$ is a set of orthonormal sequence, ${\rm span}\{\psi_j({\bf x}), j=1,\cdots, k-1\}={\rm span}\{{\xi}_j({\bf x}), j=1,\cdots, k-1\}$, and ${\mathbb E}[{\psi}_j({\bf x})]=0$, $\forall \, j\le k-1$. Define $h_k({\bf x})=\xi_k({\bf x})-{\mathbb E}[\xi_k({\bf x})]-\sum_{j=1}^{k-1}{\mathbb E}[\xi_k({\bf x}) \psi_j({\bf x})] \, {\psi}_j({\bf x}).$

$h_k({\bf x})$ is orthogonal with each of $\{{\psi}_j({\bf x}), j=1,\cdots, k-1\}$, ${\mathbb E}[h_k({\bf x})]=0$, and ${\mathbb E}[h_k^2({\bf x})]={\rm Var}[\xi_k({\bf x})]-\sum_{j=1}^{k-1}\{{\mathbb E}[\xi_k({\bf x}) \, {\psi}_j({\bf x})]\}^2>0;$ Otherwise $\xi_k({\bf x})$ is a combination of $\{{\psi}_j({\bf x}), j=1,\cdots, k-1\}$, a contradiction. Let ${\psi}_k({\bf x})=\frac{\displaystyle1}{\displaystyle{\rm s.d.}(h_k({\bf x}))} \, h_k({\bf x}).$ Hence, the assertion holds for $k$. By the method of induction, the assertion holds for any given $k$ and therefore for the sequence.

Noting that ${\rm span}\{\psi_j({\bf x}), j=1,\cdots\}={\rm span}\{{\xi}_j({\bf x}), j=1,\cdots\}$, $m({\bf x})$ can be expanded as an orthogonal infinite series in terms of $\{{\psi}_j({\bf x}), j\geq 1\}$ iff it can be expanded as an orthogonal infinite series in terms of $\{{\xi}_j({\bf x}), j\geq 1\}$. So, $m({\bf x})\in {\mathcal S}$ iff $m({\bf x})\in {\mathcal S}_\pi$. It is worth mentioning that in the Gram-Schmidt procedure, all mathematical moments in population can be replaced by their sampling versions, such as replacing ${\mathbb E}\left[\psi_j(\mathbf{x}) \, \psi_l(\mathbf{x})\right]$ by $\frac{1}{n} \sum_{i=1}^n \psi_j(\mathbf{x}_i) \, \psi_l(\mathbf{x}_i)$ for all $(j,l)$, in practice, when there is no prior knowledge about the distributional structure of the data under study.

Nonlinear regression

Estimation of nonlinear regression

Consider a general nonlinear regression model of the form:

where the functional form of $g(\cdot, \bm{\theta}_0)$ is assumed to be parametrically known, but is indexed by ${\bm \theta}_0$ as a vector of unknown parameters, and $(\varepsilon, \mathbf{x})$ is the same as defined in Section (ref), in which the $p$--dimensional ${\bm \theta}_0$ is an interior point of a compact set $\Theta\subset \mathbb{R}^p$.

The relevant literature (see, amemiya1977, amemiya1985, BT1985, newey1990, for example) discusses endogeneity issues for model ((ref)). Before we discuss about how to identify and then estimate $\bm{\theta}_0$, we provide the following examples.

{\bf Example A.1}: Consider the following log(wage) and education model:

where $y$ denotes the log(wage), $x_{1}$, $x_{2}$ and $x_2$ denote the education, experience, and experience$^2$, respectively, and $\varepsilon$ is uncorrelated with $x_2$ and $x_2^2$, but is correlated with $x_1$ (see, for example, wooldridge2016), and $\bm{\theta}_0 = \left(\theta_{00}, \theta_{01}, \theta_{02}, \theta_{03}\right)^{\top}$ is a vector of unknown parameters.

Appendix E shows that the proposed SP estimation method works well numerically when $m(\mathbf{x})$ contains second--order polynomial terms of $\mathbf{x}$.

{\bf Example A.2}: Consider the following additive case where

where $\bm{\theta}_0 = \left(\bm{\alpha}_{01}^{\top}, \cdots, \bm{\alpha}_{0q}^{\top}; \gamma_{01}, \cdots, \gamma_{0q}\right)^{\top}$ is a vector of unknown parameters of interest, and each $g_j(\cdot)$ is a known function commonly used in empirical applications.

Model ((ref)) covers a class of important models often used in the neural network literature for the case where $g_j(\cdot) = \sigma(\cdot)$ is chosen as the so--called activation function (see, for example, BK2019). While permitting possible endogeneity, Example E3 of Appendix E shows that the SP method works well numerically for a special case of model ((ref)).

We now extend our SP method to identify and estimate $\bm{\theta}_0$ in model ((ref)). Recall $m(\mathbf{x}) = {\mathbb E}[\varepsilon|\mathbf{x}]$ and then rewrite model ((ref)) as follows:

Note that the introduction of model ((ref)) not only covers the linear model in ((ref)), but also facilitates the discussion of a wide range of nonlinear and non--separable econometric models as outlined in Appendix A.3 of the supplementary document.

As in model ((ref)), we rewrite model ((ref)) as follows:

where $r(\mathbf{x}) = \sum_{j=k+1}^{\infty} \psi_j(\mathbf{x}) \, \gamma_j$, and the true version: $(\bm{\theta}_0, \bm{\gamma}_0)$ is chosen such that to minimize

Letting $\frac{\partial S(\bm{\theta}, \bm{\gamma})}{\partial \bm{\theta}}|_{(\bm{\theta} = \bm{\theta}_0, \bm{\gamma} = \bm{\gamma}_0)}=0$ and $\frac{\partial S(\bm{\theta}, \bm{\gamma})}{\partial \bm{\gamma}}|_{(\bm{\theta} = \bm{\theta}_0, \bm{\gamma} = \bm{\gamma}_0)}=0$, we then have

which imply

where $\mathbf{g}_1(\mathbf{x}, \bm{\theta}_0) =\frac{\partial g(\mathbf{x}, \bm{\theta})}{\partial \bm{\theta}}|_{\bm{\theta} = \bm{\theta}_0}$.

Recall ${\mathcal L}^2$, ${\mathcal S}$, ${\mathcal P}_{\mathcal S}$ and ${\mathcal M}_{\mathcal S}$ as defined in the same way as in Section 2.2 above. Recall also that the operator ${\mathcal P}_{\mathcal S}$ projects any element of ${\mathcal L}^2$ into ${\mathcal S}$, while the operator ${\mathcal M}_{\mathcal S}$ projects any element of ${\mathcal L}^2$ into ${\mathcal S}^\bot$. In particular, we have ${\mathcal M}_{\mathcal S}(\varepsilon)=\varepsilon-{\mathcal M}_{\mathcal P}(\varepsilon)=e$.

If we use the projection mapping operator, equation ((ref)) reduces to

which is the first order condition of $\min_{{\bm \theta}}{\mathbb E}[{\mathcal M_S}(y-\mathbf{g}(\mathbf{x}, \bm{\theta}))]^2$ at $\bm{\theta} = \bm{\theta}_0$. It follows from model ((ref)), ${\mathcal M_S}(y-\mathbf{g}(\mathbf{x}, \bm{\theta}_0))=e$, which is the corresponding nonlinear parametric model without endogeneity. Therefore, the conventional NLS applies. In a similar way to the proof of Lemma (ref) above, it can be shown that the identifiability of $\bm{\theta}_0$ by ((ref)) is invariant to the choice of $\{\psi_j(\cdot): j\geq 1\}$.

In the case of $g(\mathbf{x}, \bm{\theta}_0) = \mathbf{x}^{\top} \bm{\beta}_0$ in Section (ref), Equation ((ref)) reduces to Lemma (ref)(ii) by requiring

to ensure that $\bm{\beta}_0 = {\bm\Sigma}_x^{-1} \, \bm{\Sigma}_{xy}$ is identifiable, where $\bm{\Sigma}_{xy} = {\mathbb E}[\mathbf{x} \, y] - \sum_{j=1}^\infty {\mathbb E}[{\bf x} \, \psi_j({\bf x})] \, {\mathbb E}[y \, \psi_j({\bf x})]$.

We now introduce the following assumption.

Therefore, operating ${\mathcal M_S}$ on both sides of model (ref) yields ${\mathcal M_S}(y)={\mathcal M_S}(g({\bf x},{\bm \theta}_0))+e$. This motivates the following observation:

where the equality holds if and only if ${\mathbb E}\{[{\mathcal M_S}(g({\bf x},{\bm \theta}_0)-g({\bf x},{\bm \theta}))]^2\}=0$. Thus, under Equation ((ref)) and Assumption (ref), ${\bm \theta}_0$ is the unique minimum point of ${\mathbb E}\{[{\mathcal M_S}(y-g({\bf x},{\bm \theta}))]^2\}$, so that ${\bm \theta}_0$ is uniquely identifiable and NLS is applicable to estimate it.

Given $\{(y_i, {\bf x}_i), i=1, \cdots,n\}$, a sampling version of model (ref) is as follows:

Note that, under Assumption (ref), we have $m(\mathbf{x})=\sum_{j=1}^\infty \psi_j(\mathbf{x}) \, \gamma_j$ with $\gamma_j={\mathbb E}[m(\mathbf{x}) \, \psi_j(\mathbf{x})]$, and for a given truncation parameter $k>1$, define the partial sum $m_k(\mathbf{x})=\mathbf{V}_k(\mathbf{x})^\top {\bm\gamma}$, where $\mathbf{V}_k(\mathbf{x})=(\psi_1(\mathbf{x}), \cdots, \psi_k(\mathbf{x}))^\top$ and ${\bm\gamma}=(\gamma_1, \cdots, \gamma_k)^\top$.

We also define $\delta_k(\mathbf{x})=\sum_{j=k+1}^\infty \gamma_j \psi_j(\mathbf{x})$ for us to rewrite ((ref)) as

which can be written in matrix form:

where $\mathbf{G}({\bm \theta}_0)=(g({\bf x}_1, {\bm \theta}_0), \cdots, g({\bf x}_n, {\bm \theta}_0))^\top$.

Defining ${\bf P}_v={\bf V}({\bf V}^\top {\bf V})^{-1}{\bf V}^\top$ and ${\bf M}_v={\bf I}_n-{\bf P}_v$, we then have from (ref) that ${\bf M}_v({\bf y}-\mathbf{G}({\bm \theta}_0))={\bf M}_v({\bm\delta}+{\bf e})$. Letting $L_n( {\bm \theta})=\frac{\displaystyle1}{\Dn}\|{\bf M}_v({\bf y}- \mathbf{G}({\bm \theta}))\|^2$, the estimator $\widehat{\bm\theta}$ is then defined by

The additional condition $\lambda(\{x: \; g(x, {\bm \theta})\ne g(x, {\bm \theta}_0)\})>0$ is necessary for the asymptotic consistency of the NLS estimator, as it helps identify ${\bm \theta}$ from ${\bm \theta}_0$ when $g(x, {\bm \theta})\ne g(x, {\bm \theta}_0)$ for a set of points, $x$'s, whose measure is greater than zero.

Let $\mathbf{S}_n({\bm \theta})= \frac{\displaystyle\partial}{\displaystyle\partial {\bm \theta}}L_n( {\bm \theta})=-\frac{\displaystyle2}{\Dn}\frac{\displaystyle\partial}{\displaystyle\partial {\bm \theta}} \mathbf{G}({\bm \theta})^\top {\bf M}_v({\bf y}- \mathbf{G}({\bm \theta}))$ and

be the score function and Hessian matrix of $L_n( {\bm \theta})$, respectively. Before we establish Theorem (ref) below, we introduce the following assumption.

The unknown quantities involved in Theorem (ref) can be consistently estimated respectively by their sampling versions. It can be shown that $\bm{\Sigma}_g$ is invariant to the choice of $\{\psi_j(\cdot): j\geq 1\}$ in an analogous way to Lemma 2.2. We also have the following remarks.

In view of the above discussion in Section (ref).2 on weak endogeneity, we also propose a test for a full level of exogeneity versus a wide range of nonlinear weak endogeneity.

Testing for nonlinear weak endogeneity

Consider $y =g(\mathbf{x}, \bm{\theta}_0) + \varepsilon = g(\mathbf{x}, \bm{\theta}_0) + m_n(\mathbf{x}) + e$ under the following null hypothesis:

with $m_n(\mathbf{x})={\mathbb E}[\varepsilon|\mathbf{x}]$ satisfies ${\mathbb E}[m_n^2({\bf x})]\rightarrow 0$ as $n\rightarrow \infty$.

To test $H_0: \, {\mathbb P}(m_n({\bf x}) =0)=1$, in a similar fashion to that discussed in Section (ref), we can establish an asymptotic normality of $T_n$ that is defined in the same way as in ((ref)) in Section (ref) with $\widetilde{e}_i = y_i - g({\bf x}_i, \widehat{\bm\theta})$.

Equation ((ref)) shows that $T_n$ is asymptotically normally distributed under $H_0$. The proof of Theorem (ref) is given in Appendix C, where it can also be shown $\sigma_n^2= \frac{2 \,\widehat{\sigma}_e^4}{3} \cdot n^2 \, k_{\max}^3\, (1+ o(1))$.

To show the consistency of our testing statistic under a sequence of local alternatives, we test

where positive sequence $a_n\to 0$ with certain rate while $m(\mathbf{x})\in {\mathcal L}^2$ and ${\mathbb E}[m^2({\bf x})]>0$.

It is pointed out that a sequence of local alternatives covers a wide range of weak endogeneity. The following theorem establishes the consistency of the proposed test under $H_1$.

The proofs of Theorems (ref)--(ref) are given in Appendix D. Discussions of Theorems (ref) and (ref) are similar to those for Theorems (ref) and (ref) in Section (ref) above.

Nonlinear and non--separable models

We now show that the proposed SP method can also be extended in a unified way to a number of non-- and semi--parametric regression models associated endogeneity issues, such as

as a semiparametric regression model, along with the following nonparametric regression:

Appendices (ref).1 and (ref).2 below outline the main ideas and steps about how to estimate models ((ref)) and ((ref)) consistently and unbiasedly.

In Appendix (ref).3, we will also discuss one class of binary models of the form:

where the quantities are the same as before.

Semiparametric regression

Observe that model ((ref)) may be motivated as follows:

where a semiparametric projection gives us the following model:

where $\mathbf{x}_v$ is either an observed variable or an IV available to the econometrician satisfying

under which we may identify and estimate $\bm{\alpha}_0$ consistently in the same way as has been done in the relevant literature. Model ((ref)) with ((ref)) has been employed in Example 5.1.

When condition ((ref)) is not satisfied, we consider model ((ref)) for the case where

Note that the proposed approach below is valid regardless of whether ${\mathbb E}[\varepsilon|\mathbf{x}_u]\neq 0$ or ${\mathbb E}[\varepsilon|\mathbf{x}_u]=0$. Let $q_1(\mathbf{x}_v) = {\mathbb E}[y|\mathbf{x}_v]$, $\mathbf{q}_2(\mathbf{x}_v) ={\mathbb E}[\mathbf{x}_u|\mathbf{x}_v]$, $\widetilde{\mathbf{x}}_u=\mathbf{x}_u - \mathbf{q}_2(\mathbf{x}_v)$, $\widetilde{y}=y - q_1(\mathbf{x}_v)$, $q_3(\mathbf{x}_v) = {\mathbb E}[\varepsilon|\mathbf{x}_v]$, $\eta = \varepsilon - q_3(\mathbf{x}_v)$, $m(\widetilde{\mathbf{x}}_u) = {\mathbb E}[\eta|\widetilde{\mathbf{x}}_u]$ and $e=\eta - m(\widetilde{\mathbf{x}}_u)$. Then model ((ref)) can be written as

with ${\mathbb E}[e|\widetilde{\mathbf{x}}_u]=0$. Notationally, model ((ref)) is the same as model (2.3) with $(y, \mathbf{x})$ being replaced by $(\widetilde{y}, \widetilde{\mathbf{x}}_u)$, respectively. The main parameter--of-interest, $\bm{\alpha}_0$, can then be estimated in the same way as in Section 2.1 based on a sampling version of the form:

where $\widetilde{y} = y - \widehat{q}_1(\mathbf{x}_v)$ and $\widetilde{\mathbf{x}}_u = \mathbf{x}_u - \widehat{\mathbf{q}}_2(\mathbf{x}_v)$, in which $\widehat{q}_1(\cdot)$ and $\widehat{\mathbf{q}}_2(\cdot)$ can be constructed as in Section 2.1, before $m(\cdot)$ can be estimated in the same way as in Section 2.1.

It is noted that our discussion covers the case where $\mathbf{x}_u$ is a binary variable, and $\mathbf{x}_v$ reduces to a univariate fixed--design variable of the form $\mathbf{x}_v = \tau \in [0,1]$.

Nonparametric regression

Consider model ((ref)), where $g(\mathbf{x})$ is an unknown function of interest, and the error term $\varepsilon$ satisfies ${\mathbb E}[\varepsilon]=0$, but ${\mathbb E}[\varepsilon|\mathbf{x}] \neq 0$. In such cases, there are several different approaches proposed in dealing with potential endogeneity issues. One of the existing methods is the so--called “Nonparametric IV" approach, as discussed and reviewed in a recent paper by cff2025, for example, which can be summarized as follows.

Let $\mathbf{z}$ be an IV such that ${\mathbb E}[\varepsilon|\mathbf{z}]=0$. One may then use ((ref)) to derive

for which $g_1(\mathbf{z})={\mathbb E}[y|\mathbf{z}]$ can be estimated by an existing nonparametric method before one might be able to recover $g(\mathbf{x})$ by solving an inverse problem ${\mathbb E}[g(\mathbf{x})|\mathbf{z}]={g}_1({\bf z})$.

As stressed before, the main drawback is that such an approach depends on the availability and validity of $\mathbf{z}$ as the IV, in addition to addressing ill--posed inverse issues.

Our discussion is to offer an alternative by decomposing $\varepsilon$ and rewriting model ((ref)) as

where ${\mathbb E}[e|\mathbf{x}]=0$ is automatically satisfied.

We first consider the case where $m(\mathbf{x}) = {\mathbb E}[\varepsilon|\mathbf{x}] = m(\mathbf{x}, \bm{\gamma}_0)$. Model ((ref)) then becomes

which is the same notation as discussed in Appendix A.2.1, where $\widetilde{g}_1(\mathbf{x}) = g(\mathbf{x}) - \mu_0$ with $\mu_0 = {\mathbb E}[g(\mathbf{x})]$. Both the identifiability and estimation of $(\mu_0, \bm{\gamma}_0, \widetilde{g}_1(\cdot))$ follows from that for $(\bm{\theta}_0, m(\cdot))$ discussed in Appendix A.2.1.

It is pointed out that model ((ref)) allows the case where $g(\mathbf{x})$ is nonlinear in $\mathbf{x}$ and $m(\mathbf{x}, \bm{\gamma}_0) = \mathbf{x}^{\top}\bm{\gamma}_0$. Consequently, the case where $(\varepsilon, \mathbf{x})$ follows a joint Gaussian distribution can be covered in ((ref)).

In the case where $m(\mathbf{x})={\mathbb E}[\varepsilon|\mathbf{x}]$ is specified nonparametrically, substantially new developments are required. We therefore wish to write up such details into a different paper.

Binary regression

Regression models with binary dependent variables have many theoretical investigations and empirical studies, and the relevant literature is comprehensive, such as lewbel1998, and blundell2004, for example. To show that the proposed SP method is also useful to address certain types of endogeneity involved in binary models, we start with the following linear probability model:

where the exogeneity case of ${\mathbb E}[\varepsilon|\mathbf{x}]=0$ has been discussed in the relevant literature, such as Chapter 11 of jsmw2018. \, We then rewrite model ((ref)) as

where $m(\mathbf{x}) = {\mathbb E}[\varepsilon|\mathbf{x}]\neq 0$, $e = \varepsilon - m(\mathbf{x})$, and ${\mathbb E}[e|\mathbf{x}]=0$ by definition. Model ((ref)) can then be estimated in the same way as in Section 3, although $y$ is now a binary response variable.

Meanwhile, we consider a non--separable binary model of the form:

which can be rewritten as ${\mathbb E}[y|\mathbf{x}] = {\rm P}_{e}(e\leq \mathbf{x}^{\top} \bm{\theta}_0 - m(\mathbf{x})) = F_{e}(\mathbf{x}^{\top} \bm{\theta}_0 - \mathbf{v}(\mathbf{x})^{\top} \bm{\gamma}) = F_{e}(\mathbf{w}_{-}(\mathbf{x})^{\top} \bm{\theta}_{-})$ when ignoring the approximation error term and $m(\mathbf{x})$ is replaced by $\mathbf{v}(\mathbf{x})^{\top} \, \bm{\gamma}$ as in Section 2, where $F_{e}(\cdot)$ denotes the cumulative distribution function (CDF) of $e$, $(\mathbf{v}(\mathbf{x}), \bm{\gamma})$ is the same as in Section 2, $\mathbf{w}_{-}(\mathbf{x}) = \left(\mathbf{x}^{\top}, -\mathbf{v}^{\top}(\mathbf{x})\right)^{\top}$ and $\bm{\theta}_{-} = \left(\bm{\theta}_0^{\top}, \bm{\gamma}^{\top}\right)^{\top}$.

When $F_{e}(\cdot)$ is parametrically known, the parameter--of--interest, $\bm{\theta}_0$, can then be estimated consistently and unbiasedly by MLE. When $F_e(\cdot)$ is nonparametrically unknown, it can be expanded by an infinite sum of known CDFs before a semiparametric MLE estimation method may be developed.

Discussion of other estimation methods

SP estimation efficiency

Consider the case of $d=1$ and $\sigma_i^2 \equiv \sigma_e^2$ for notational simplicity in the following derivations. Recall the standard OLS estimator by

where $s_n \equiv: \sum_{i=1}^n \mathbf{x}_i \, \mathbf{x}_i^{\top}$.

Let $\bm{\Sigma}_{xm} = {\mathbb E}[e^2] \, {\mathbb E}\left(\left[\mathbf{x} \, \mathbf{x}^{\top}\right]\right) + {\mathbb E}\left(\left[\mathbf{x} \, m(\mathbf{x}) - {\ell}_{\rm endo}\right] \, \left[\mathbf{x} \, m(\mathbf{x}) - {\ell}_{\rm endo}\right]^{\top}\right)$, $\bm{\Sigma}_{\rm LS} = {\mathbb E}^{-2} \left[\mathbf{x}_1^2\right]\, \bm{\Sigma}_{xm}$ and ${\ell}_{\rm endo} = {\mathbb E}[\mathbf{x} \, m(\mathbf{x})]$, and $\bm{\Sigma}_{\rm SP} = \sigma_e^2 \, \bm{\Sigma}_x^{-1}$ with $\bm{\Sigma}_x = {\mathbb E}[\mathbf{x}_1^2] - \sum_{j=1}^{\infty} {\mathbb E}^2\left[\mathbf{x}_1 \, \psi_j(\mathbf{x}_1)\right]$. We now show that $\bm{\Sigma}_{\rm SP} \leq \bm{\Sigma}_{\rm LS}$ as follows:

when ${\mathbb E}\left(\left[(\mathbf{x}_1 \, m(\mathbf{x}_1) - {\ell}_{\rm endo})\right]^2\right) \geq \frac{\bm{\sigma}_{11} \, \bm{\sigma}_{12}}{\bm{\sigma}_{11} - \bm{\sigma}_{12}} \, \sigma_e^2$, where $\bm{\sigma}_{11} = {\mathbb E}[\mathbf{x}_1^2]$ and $\bm{\sigma}_{12} = \sum_{j=1}^{\infty} {\mathbb E}^2\left[\mathbf{x}_1 \, \psi_j(\mathbf{x}_1)\right]$.

Equation ((ref)) remains true even in the case of ${\ell}_{\rm endo}=0$ and $\bm{\sigma}_{12} = 0$. In other words, $\widehat{\bm{\beta}}_{\rm SP} \equiv: \widehat{\bm{\beta}}$ is more efficient than $\widehat{\bm{\beta}}_{\rm LS}$ even when ${\ell}_{\rm endo}=0$.

GMM estimation method

Recall the notation and symbols introduced in Section 3. Letting $\mathbf{Q}_{d\times d}$ be a known positive definite weight matrix to be chosen by the user, we estimate $\bm{\beta}_0$ by

which offers a closed--form expression as follows:

It is known from the discussion in Section 3 that

We then have

The second term is $o_P(1)$ as $k\to \infty$ under Assumption 3.1(iii)(iv) as shown in the proof of Theorem 3.1. For the first term, the conditional covariance matrix is

as $(n,k)\to(\infty, \infty)$ under the conditions of Theorem 3.1(ii) since ${\bf Q}$ and $\bm{\Sigma}_x(k)$ are invertible. The result is the same as Theorem 3.1.

Estimation of weakly identified linear models

We start with the linear case of $m(\mathbf{x}) = \left(\mathbf{x} - {\mathbb E}[\mathbf{x}]\right)^{\top} \bm{\gamma}_0$. Let us define a truncated version of $\bm{\Sigma}_x$ of the form:

If we assume that we can expand $\mathbf{x} - {\mathbb E}[\mathbf{x}]= \sum_{j=1}^{\infty} \psi_j(\mathbf{x}) \, \bm{\gamma}_j$ by the same orthonormal series: $\{\psi_j(\cdot): j\geq 1\}$ as used in Section 2, we have each given $k\geq 1$

which has a reduced rank when ${\mathbb E}[\mathbf{x}]\neq 0$ and $d>1$. In this case, the proposed LASSO selection estimation method in Section 5.1, which allows for $m(\mathbf{x})$ to be linear in $\mathbf{x}$, addresses such reduced--rank issues. The finite--sample evaluation results in Section 5 and Appendix (ref) support the LASSO selection method.

In the case of $d=1$ and ${\mathbb E}[x]\neq 0$, it can be seen that the SP method itself is directly applicable, and Theorem 3.1 remains true. Example B.2.2 of Appendix (ref) shows that without using the LASSO selection method, the SP selection method works well with commonly used functions, including $m({x}) = \left({x} - {\mathbb E}[{x}]\right) {\gamma}_0$ and second--order polynomial functions.

When ${\mathbb E}[\mathbf{x}] = 0$, model (2.3) reduces to

which means that one may only be able to correctly identify $\left(\bm{\beta}_0 + \bm{\gamma}_0\right) = {\mathbb E}^{-1}\left[\mathbf{x} \, \mathbf{x}^{\top}\right] \, {\mathbb E}\left[\mathbf{x} \, y\right]$ collectively, rather than $\bm{\beta}_0$ individually.

For model ((ref)) in the case of ${\mathbb E}[\mathbf{x}] = 0$, meanwhile, equation ((ref)) implies that as $k\rightarrow \infty$

Due to ((ref)), for model ((ref)), we cannot assume Assumption 2.1(iii). Instead we replace it by {\bf Assumption 2.1(iii)*}: Assume that there are a positive definite matrix of real numbers, $\mathbf{Q}_n$, and another positive definite matrix, $\bm{\Sigma}_{x, \ast}$, such that $\lambda_{\rm \min}(\mathbf{Q}_n) \rightarrow \infty$, $ n \, \lambda_{\rm \min}(\mathbf{Q}_n^{-1}) \rightarrow \infty$ and $\mathbf{Q}_n \, \bm{\Sigma}_x(k)\rightarrow_P \bm{\Sigma}_{x, \ast}$ as $n\rightarrow \infty$, where $\lambda_{\min}(A)$ denotes the smallest eigenvalue of matrix $A$.

In view of the proof of Theorem 3.1, replacing Assumption 2.1(iii) by Assumption 2.1(iii)*, it can then be established that the following asymptotic normality holds:

which still offers a consistent estimator for $\bm{\beta}_0$, with a reduced rate of convergence of an order of $\left(\sqrt{n} \, \mathbf{Q}_n^{-1/2}\right)$, slower than $\sqrt{n}$, in which

Finite--Sample Evaluations

Implementational Issues

We discuss several important issues about how to implement the proposed SP estimation procedure in practice. We offer our recommendations on the choice of orthonormal series functions and the truncation parameters involved in the proposed SP estimation.

Dimension reduction

The main model and estimation method proposed in Section 3 remains valid for the case where the dimensionality, $d$, is large but fixed in theory, although we assume in Section 2 that we focus on the case where $d$ is small. To explain this in a bit more detail, we recall from Section 2 that we expand $m({\bf x})$ as $m(\mathbf{x}) = \sum_{j=1}^{\infty} \psi_j(\mathbf{x}) \, \gamma_j$, in which the dimensionality of $\mathbf{x}$ is only involved in the chosen series $\{\psi_j(\mathbf{x}): j\geq 1\}$. For ease of implementation, we suggest using $\psi_j(\mathbf{x}) = \prod_{k=1}^d \psi_{jk}(x_k)$ when $d\geq 2$, where $\{\psi_{jk}(\cdot): k\geq 1\}$ is an array of univariate series functions, and $\mathbf{x} = (x_1, \cdots, x_d)$.

Meanwhile, some other dimension reductions might be employed, such as an additivity structure of the form: $m(\mathbf{x}) = \sum_{j=1}^d m_j(x_j)$ directly as in ln1995, in which it is expected that the construction of $\mathcal{S}$ should be a direct sum of $\mathcal{S}_j$ for $1\leq j\leq d$. There are also dimension reduction methods proposed by involving certain types of single--index or multi--index modelling methods (see, dg2025, for proposing an additive single--index structure form for each $\psi_j(\cdot)$). Our experience with the finite--sample studies in Section 5 and Appendix (ref) shows that the multiplicative form: $\psi_j(\mathbf{x}) = \prod_{k=1}^d \psi_{jk}(x_k)$ works well and better than some other competing forms available in the relevant literature.

Choice of series functions

In both theory and practice, we need not require orthogonality on $\{\psi_l(\cdot): l\geq 1\}$, although the orthogonality assumption simplifies the notation involved in the theoretical derivations. In practice, our experience suggests using either one of the following series, or a mixture of both.

1. \, The probabilist's Hermite polynomials $\{H_j(x),j\ge 1\}$ are a set of orthonormal basis functions defined on $L^2(\mathbb{R}, \exp(-x^2/2))$ for the univariate setting;

2. \, Let the trigonometric polynomials be $p_j(x) = \sqrt{2}\cos(\pi j x)$ with $j\ge 1$. Lemma E.1 of Appendix E in the online supplement shows that $\{p_j(x): j\ge 1\}$ is an orthogonal set of basis functions on $L^2([a, b])$ as long as $(a,b)$ are different integers. For the case $a=0$ and $b=1$, $\{p_j(x): j \ge 1\}$ constitutes an orthonormal family.

3. \, In the multivariate setting, we propose using the following form of either $p_j(\mathbf{x}) = \prod_{k=1}^d p_{jk}(x_k)$ or $H_j(\mathbf{x}) = \prod_{k=1}^d h_{jk}(x_k)$ in simulations and empirical applications, where $p_{jk}(\cdot)$ and $h_{jk}(\cdot)$ are the corresponding univariate functions that may be chosen as in Steps 1 and 2 above.

As discussed in Section 4.2.3 and Appendix C.4 of dg2025, we need not know the distributional structure of the data under analysis in practice as long as the support of $\{\mathbf{x}_i: i\geq 1\}$ becomes available to the practitioner.

In theory, the choice of $\{\psi_j(\cdot): j\geq 1\}$ can be flexible as long as $\bm{\Sigma}_x>{\bf 0}$. As discussed in Section 5.1 of the main submission, moreover, the proposed LASSO selection not only facilitates the choice of an optimal orthonormal series for $\{\psi_j(\cdot): j\geq 1\}$ in practice, but also helps address possible reduced--rank issues as alluded in Appendix (ref) above.

Main simulations

Finite--sample properties of the estimation theory

In this section, we consider a number of scenarios to demonstrate the finite-sample performance of the proposed SP estimation method under orthonormality. To show that the orthonormality and even the orthogonality on $\{\psi_j(\cdot): j\geq 1\}$ may all be relaxed, we present extensive numerical evaluations in Appendix E of the online supplementary document to demonstrate that the SP method still works well numerically. We provide the code at \url{https://github.com/pengbin430/SIV_2026/tree/main}.

{\bf Example B.2.1}: We consider the following data generating processes:

In Appendix (ref).3 below, we show that $E[x_im(x_i)]\ne 0$ and $E[m(x_i)]=0$ for all cases, so the endogeneity exits. Case C can also be considered as a special case of the NLS framework discussed in Appendix A.2. Both Cases B and C have truncation residuals. As shown in the justification, Case B has more sparsity in $\{\gamma_j \}$, so it is reasonable to expect our approach has better finite sample performance for Case B.

After identifying $\phi_j(x)$'s with nonzero coefficients from the LASSO method, we then construct $\{\psi_j(\cdot): j\geq 1\}$ to run SP to obtain $\widehat{\beta}_{\rm SP}$. We alway enforce $\mathbf{V}_k(\cdot)$ to include the constant term 1 in order to remove the intercept.

{

}

We repeat the above procedure $R=1000$ times and then report the following measures in Table (ref):

where $\overline{\beta} = \frac{1}{R}\sum_{r=1}^R\widehat{\beta}_r$ with $\widehat{\beta}_r \in\{\widehat{\beta}_{\rm LS},\widehat{\beta}_{\rm SP}\}$ in each replication. Additionally, we calculate $\left|\frac{\Delta \beta_{SP}}{\Delta \beta_{LS}}\right|$ in Table (ref) in order to show the relative magnitude of the biases associated with both OLS and SP methods. As indicated in Table (ref), the OLS method is apparently biased and has large standard deviation. The SP method works well for all cases. Although the magnitudes of biases vary in different cases, the ratio $\left|\frac{\Delta \beta_{SP}}{\Delta \beta_{LS}}\right|$ is reasonably stable across all cases. The only exception is the ratio associated with $\beta_1$ of Case D, which should be expected due to the fact that $w_i$ is an exogenous variable.

{

}

We then examine the LASSO selection method to further reiterate our discussion in the above sections. Specifically, in Table (ref), we report the following measure: $P_j =\frac{1}{R}\sum_{r=1}^R I(\phi_j\in S_r)$, where $S_r$ stands for the set $S$ selected by the LASSO method in the $r^{th}$ replication. Therefore, $P_j$ measures the probability of $\phi_j$ being selected over $R$ replications.

It is clear that that the LASSO method works well. A few findings emerge. Take Case C as an example. Justification of Example B.2.1 given in Appendix (ref) shows that the odd indexed $\phi_j$'s have nonzero coefficients, while the even indexed $\phi_j$'s have coefficients 0. Consequently, Table (ref) shows that the odd indexed $P_j$'s are often selected with high probability, while the even indexed $P_j$'s are barely selected.

We then add some simulation results for the case where $m(\mathbf{x})$ itself contains a linear component in each case under examination in Example B.2.2 below.

{\bf Example B.2.2}: We now consider model $y_i = x_i \, \beta_0 + \varepsilon_i$, $\varepsilon_i = m(x_i) + e_i$, $e_i\sim N(0,1)$ and $x_i \sim U(0, \pi)$ under the following scenarios:

Apparently, the restrictions: ${\mathbb E}[m(x_i)]=0$ and ${\mathbb E}[x_i \, m(x_i)]\neq 0$ can be verified in a similar way to those for Example B.2.1. We summarize the simulation results in Table (ref) below, wherein the biases and standard deviations (in parentheses) are calculated in exactly the same way as before.

Table (ref) obviously shows that the SP method works well numerically. Moreover, the numerical results are better than those in Example B.2.1, probably because the elementary functional forms of $m(\cdot)$ considered in Cases E--G were more accurately approximated by the trigonometric series $\left\{\psi_j(x) = \sqrt{\frac{2}{\pi}} \, \cos(j\, x): j\geq 1\right\}$ for the case of $x\in U[0, \pi]$.

Finite--sample properties of the testing theory

In this section, we evaluate the finite--sample property of the test proposed in Section 3.2.

{\bf Example B.2.3}: Without loss of generality, we consider Cases A and B of Example B.2.1, and slightly make the following modification in order to evaluate size and power respectively.

where $x_i \sim U(0, \pi)$, and $a_{n,j} =c_j \, \sqrt{\frac{k_{\max}}{n}}$ for $1\leq j\leq 4$. Without of loss generality, we chose $c_1\equiv 0$ to evaluate the sizes, and $c_2 \equiv 0.75$, $c_3 \equiv 1$, and $c_4 \equiv 1.25$ to evaluate the local powers. Accordingly, we have $\{\psi_{j}(x)=\sqrt{2/\pi}\cos(jx)\mid j \ge 1 \}$, which are defined precisely in multiple places. When calculating the test statistic, we chose $k_{\min}=\lceil \log \log n\rceil$ and $k_{\max}=\lceil 3 \log \log n\rceil$. By construction, we have

so the requirement to ensure local power properties is fulfilled. We generate $m(x)$ in the following two forms, which are almost identical to those in Example B.2.1.

To improve the finite sample performance of the size function in each case, we propose using the following wild bootstrap procedure (see, for example, gg2008) to obtain bootstrapping critical values.

After $R=1000$ simulation replications, we report the following rejection rate for $n=200,300, 400$ respectively: $\text{RJ} =\frac{1}{R} \sum_{r=1}^R I(L_{n,r}\not \in \mathscr{C}_{n,r}),$ where $L_{n,r}$ and $\mathscr{C}_{n,r}$ respectively stand for the values of $L_n$ and $\mathscr{C}_n$ obtained in the $r^{th}$ replication. The results are summarized in Table (ref).

In Table (ref), the rejection rates with $c_1=0$ are around the nominal rejection rate (i.e., 5%). When using $c_2$, $c_3$ and $c_4$, the rejection rates reflect the local power of the proposed test. It is not surprising that the local power varies across two cases. In Case B, the rejection rates are slightly lower than those in Case A, which might be due to the impact of the truncation residual. In summary, Table (ref) shows that the proposed test can detect the weakest possible endogeneity at an optimal order of $n^{-1/2} \, \delta_n$ for such $\delta_n$ that diverges to $\infty$ at one of the slowest possible rates of an order of $\log(\log(n))$.

Verification of the simulation designs in Examples B.2.1 and B.2.2

We consider Case C of Example B.2.1 only here. As the derivations for Cases A, B and D of Example B.2.1 are similar and simpler, and for Example B.2.2 are also similar, we omit the details. To show that $E[m(x_i)]=0$ and $E[x_im(x_i)]\ne 0$, it is sufficient to consider

so endogeneity exists.

We then show that

and

Therefore, we have shown that $m(x)$ can be fully expanded by the odd terms of $\{\cos(j \, x)\}$, while fully expanding the regressor $x^2$ requires all of $\{\cos(j \, x)\}$.

Extra empirical analysis

{\bf Example 5.1} (continued)

As alluded before, we change the basis functions to $\phi_j(x) =(x-\pi/2)^j$ for $j\geq 1$, and keep everything else unchanged. Obviously, the new set of basis functions are not even orthogonal to each other in the space $L^2([0,\pi])$. The LASSO method selects $\phi_2(x) =(x-\pi/2)^2$ only, so $\widehat{m}(x)= -0.0643(x-\pi/2)^2$.

We then obtain the following table, wherein the last two rows are corresponding to the results obtained from the 2sLS method by using $(z_1,\phi_2)$ and $\phi_2(x)=(x-\pi/2)^2$ as IVs, respectively.

It is observed that the new estimation results are very similar to those reported in Section 5. This may indicate the robustness of our method, and the choice of $\{\phi_j(\cdot): j\geq 1\}$ and whether it is orthonormal or not may not affect the corresponding estimation results very much.

Possible sources of endogeneity

As an additional evidence to show that there is some substantial nonlinear endogeneity involved in Example 5.1, we plot the estimated $m(x)$ with its 95% confidence interval (obtained via wild bootstrap) in Figure (ref). It is clear the SP method supports a clearly nonlinear functional form for the specification of $m(x)$.

{\samepage {

}}

As briefly pointed out in the previous sections of this paper, an additional novelty of the SP approach is that we are able to estimate $\varepsilon$ of model (2.1) by $\widehat{\varepsilon} = y - \mathbf{x}^{\top} \widehat{\boldsymbol{\beta}}$. As a consequence, we may be able to make use of $\widehat{\varepsilon}$ to establish an approximate regression model of the form:

for the econometrician to try to reveal possible sources of endogeneity, where each $\mathbf{w}_i$ is a vector of variables possibly containing omitted variables and/or errors in variables, which may be highly correlated with $\mathbf{x}$, and ${\eta}$ is the error term involved.

Due to the fact that $\varepsilon$ is usually unobservable and latent, it is reasonable to assume the non--separable functional form in ((ref)). If each $\mathbf{w}_i$ is available to the econometrician, the discussion in Appendices A.2 and (ref) suggests that one may be able to reveal possible sources of endogeneity by estimating the functional form of $G(\mathbf{w}; \cdot)$ with respect to $\mathbf{w}$.

Let us now come back to Example 5.1. We would like to investigate how $\varepsilon$ might be uncorrelated with $\{z_j: 1\leq j\leq 5\}$ used above. We compute the sample averages of the following quantities: $E[\widetilde{x}\cdot m(\widetilde{x})]$ and $E[z_j\cdot \varepsilon]$, where $\widetilde{x}$, $m(\widetilde{x})$ are defined in (5.3), $z_j$'s are the different IV variables adopted above, and those unobservable (i.e., $m(\cdot)$ and $\varepsilon = y- \widetilde{\alpha}_0- \widetilde{x}\, \widetilde{\beta}_0$ are replaced with their estimates via the SP method).

As shown in Table (ref), $E[\widetilde{x}\cdot m(\widetilde{x})]$ is clearly non-zero, which justifies the existence of the endogeneity. As an valid IV, we would expect $E[z_j\cdot \varepsilon]$ all to be zero. However, only $z_1$ and $z_5$ generate averages sufficiently closed to 0. Numerically, it offers a reason to explain why the results associated with $z_1$ and $(z_1, z_5)$ in Table 2 are close to that from the SP approach. However, the other quantities: $E[z_j\cdot \varepsilon]$, for $2\leq j\leq 4$, don't justify the validity of those IVs.

Finally, we estimate the following three nonparametric regressions models:

where $\widetilde{x}$, $z_1$ and $\varepsilon$ have been defined previously, $g_1(\cdot)$, $g_2(\cdot)$ and $g_3(\cdot)$ measure possible relationships among them, and $\xi_1$, $\xi_2$ and $\xi_3$ are the error terms.

{\samepage

}

The nonparametrically estimated relationships, along with their corresponding $95\%$ confidence intervals, are plotted in Figure (ref). The second plot clearly shows that the relationship between $\varepsilon$ and $z_1$ is clearly nonlinear, although being symmetrically fluctuating around zero. For simplicity, we use the GCV method of gao2002 to select the truncation parameter here. In the second case, the GCV suggests the number of truncation parameters being $20$. The second plot infers that the typical IV method might be modelled beyond linear regression.

The first and third plots together infer the relationship between $x$ and $z_1$ has a linear trend but not exactly linear, as the truncation parameters for these two regressions are 13 and 7 respectively according to the GCV approach. This further explains why the results associated the SP method and the IV method (using $z_1$) are close to each other in Table 2.

In summary, due to the unobservable nature of $\varepsilon$, it is difficult to establish and reveal a relationship between $\varepsilon$ and potential instrumental variables. Without assuming the knowledge of such relationships, the existence and the validity of potential IVs, the proposed SP approach offers a feasible and robust way to deal with certain types of endogeneity issues.

As discussed above, the proposed SP estimation method is applicable to time--series data. So we now have a look at the following application.

Endogeneity in stock return

We apply the SP method to revisit the classic market model (M1997, campbell1998econometrics). See Eq. (4.3.2) of campbell1998econometrics for example.

Consider the following linear model:

where $\varepsilon_t=m(x_t) + e_t$, $y_t$ represents daily stock return for a chosen stock, and $x_t$ stands for the market portfolio. In this study, we use daily return of the SP500 index as $x_t$. We aim to argue that endogeneity should also be taken into account for the market model.

For the response $y_t$, we collect data for Apple, Google, Meta and Microsoft from Yahoo finance covering period from the first trading day of 2012 to the last trading day of 2020, which gives $n=2169$ observations for us to conduct regressions. As in our study about the return to schooling, we normalize $x_t$ to get $\widetilde{x}_t$, so the latter belongs to $[0,\pi]$. Thus, the model used for regression is again written as

For each stock, we consider (i) The SP method; \ (ii) The LS method; and (iii) The 2sLS method using the basis functions selected by the SP method as IVs.

{

}

The estimated $m(x)$ for different companies are given as follows:

For the purpose of demonstration, we also plot $\widehat{m}(x)$ from Google with its 95% confidence interval in Figure (ref) below, which clearly shows some nonlinear fluctuation of the estimated endogeneity component represented by $\widehat{m}(\cdot)$.

Tables (ref) and (ref) summarize the results, revealing several key findings. First, Table (ref) indicates similar features about the SP estimates to those discussed in Section 5 with relatively smaller RMSEs in comparison with those of LS and 2sLS estimates. Second, Table (ref) confirms the presence of endogeneity, as the estimated values of $E[\widetilde{x} \, m(\widetilde{x})]$ are negative across all companies, which should be expected given that these four companies account for approximately more than 19% of the entire S&P 500 index.

Third, the estimated $\beta_0$ for all companies is greater than 1 indicating that these four stocks are more volatile than the market after accounting for the endogeneity.

Bias correction for weak endogeneity

This section discusses about how to eliminate the bias term by a simple jackknife method.

To show the main idea, we focus on a weak endogeneity case of the form:

where $\mu=4$, $x_i\sim U(0,\pi)$, and $e_i\sim N(0,1)$.

We partition our sample in two groups: $\mathscr{N}_1$ and $\mathscr{N}_2$, where $\sharp \mathscr{N}_j =n_j$ for $j=1,2$, and $\mathscr{N}_1\cap \mathscr{N}_2=\emptyset$ before we define the OLS estimators: $\widehat{\bm{\beta}}_{\rm LS, 1} = \left(\sum_{\mathscr{N}_1}\mathbf{x}_i\mathbf{x}_i^\top\right)^{-1} \sum_{\mathscr{N}_1}\mathbf{x}_i y_i$ and $\widehat{\bm{\beta}}_{\rm LS, 2} = \left(\sum_{\mathscr{N}_2}\mathbf{x}_i\mathbf{x}_i^\top\right)^{-1} \sum_{\mathscr{N}_2}\mathbf{x}_i y_i$.

Define $\widehat{\bm{\beta}}_{\rm BC} = c \, \widehat{\bm{\beta}}_{\rm LS, 1} + (1-c) \, \widehat{\bm{\beta}}_{\rm LS, 2}$. Let $n=n_1 + n_2$. Simple algebra shows that

Note that ${\rm Var}\left(\frac{c}{\sqrt{n_1}}\sum_{\mathscr{N}_1}\mathbf{x}_i e_i+ \frac{(1-c)\sqrt{n_1}}{n_2}\sum_{\mathscr{N}_2}\mathbf{x}_i e_i\right) = 2 \, c^2\, {\mathbb E}\left[\mathbf{x} \, \mathbf{x}^{\top}\right] \sigma_e^2$ due to requiring $c = (c-1)\frac{\sqrt{n_1}}{\sqrt{n_2}}$. Note also that while the choice of $c$ affects the variance component, it doesn't have any impact on the bias evaluation as long as $c>1$ satisfies $c = (c-1)\frac{\sqrt{n_1}}{\sqrt{n_2}}$.

We therefore choose $c=1.5$ and then $n_2=\left[\frac{n_1}{9}\right]$ in the following numerical evaluation. We consider the case of $n_1\in \{270, 405, 540\}$, so the value of $n_2$ is defined accordingly as follows. We report the absolute bias and the standard derivation based $R=1000$ replications as follows:

where $\overline{\beta} =\frac{1}{R}\sum_{r=1}^R \widehat{\beta}_r $, and $\widehat{\beta}_r$ stands for the value of $\widehat{\beta}_{\rm LS}$ or $\widehat{\beta}_{\rm BC}$ in the $r^{th}$ replication.

The results are summarized in Table (ref). It is clear that $\widehat{\beta}_{\rm LS}$ yields a much larger bias in each individual case, although $\widehat{\beta}_{\rm LS}$ has a smaller standard derivation correspondingly. The unbiased estimate $\widehat{\beta}_{\rm BC}$ has been obtained by partitioning our full sample size $n=n_1+n_2$, as a consequence, it sacrifices the standard deviations slightly.

Proofs for Sections 2, 3 and 5

{

}

\setcounter{section}{3}