Endogeneity Corrections in Binary Outcome Models with Nonlinear Transformations: Identification and Inference

-2cmEndogeneity Corrections in Binary Outcome Models with Nonlinear Transformations: Identification and Inference

\thispagestyle{empty}

Introduction

Estimating regression models in the presence of endogeneity without external instruments has become popular in econometrics (see, e.g. lewbel:1997, rigobon:2003, ebbes:2005, kleinvel:09, dong:2010, esc:16, tran:22, lewetal:23, kiviet:23, or gaowang:2023) and in many applications in business and economics. Examples include, among others, empirical marketing science (see, e.g. burmeister:15, bach:21, or zhang:24), productivity analysis (see, e.g. tran:2015 or haschka:24a) or energy economics (see, e.g. aloui:16).

The starting point of this paper is the seminal work by guptapark:2012, who propose such an estimator based on a specific copula assumption concerning the dependence between the endogenous regressor and the error term. yang:2022 and haschka:2022, haschka:24 provide extensions for the case that there is dependence between endogenous and exogenous regressors. bmw:24 consider an approach based on control functions and derive asymptotic properties for their estimator. A recent literature review can be found in parkgupta:2024, see also becker:22, papa:2022, lien:2024, or yietal:24.

We fill a gap in the literature by considering binary outcome models in the presence of endogeneity without external instruments and with possibly nonlinear dependence between endogenous and exogenous regressors. Similar to the seminal work of rivers:88, we identify the structural parameters using a control function approach. The crucial distinction is, however, that we do not require outside instruments. In particular, as a direct extension to bmw:24, we propose to estimate these models with a nonparametrically generated control function to take potential endogeneity into account. This control function is derived from a rank-based transformation of the reduced-form residuals. In contrast to bmw:24, however, we explicitly allow for nonlinear dependence between endogenous regressor and exogenous regressor in the first step, which is similar in nature to the approach taken by dong:2010 or esc:16. In doing so, we allow for parametric, semi-parametric, and nonparametric estimation of the first stage. The so-obtained residuals are transformed to ranks, the ranks are transformed with the standard normal quantile function, and the resulting term is added as a control term to the regression equation.

This way, we obtain identification of the parameters, if the endogenous regressor is a nonlinear function of one component of the error term, conditional on the exogenous regressors. In this case, the transformation in the first step may be linear or nonlinear. Moreover, there is identification if the endogenous regressor is a linear function of the error term, but the dependence in the first step is nonlinear.

The possibility of linearity in the first step is an improvement over the existing literature. For example, dong:2010 considers a special case of our model, in which the transformation in the first step regression has to be nonlinear to achieve identification. In this case, no restriction on the dependence between the latent error terms is required. This is also similar to the identification strategy in esc:16, who require nonlinearity of the first step. Moreover, our estimation strategy, based on rank-based transformations, is inspired by the seminal work of guptapark:2012 and differs from the kernel-based estimators used, for example, in dong:2010 and esc:16.

The difficulty of deriving an asymptotic theory stems from the nonparametrically generated control function. Using recent results from residual empirical processes theory (zhao:2020 and zhao:2022) for (nonparametrically) estimated normal scores, we are able to show that the estimator is consistent and asymptotically normal. We do so by establishing sufficient high-level conditions that allow for parametric, semi-, and nonparametric estimation of the reduced form regression function. Similar to pagan:1984, the sampling uncertainty of having to estimate the control function affects the sampling distribution of the estimator. We thus propose a bootstrap procedure to take it into account.

A simulation study yields numerical evidence and an empirical application on German insolvency data illustrates the usefulness of our estimator. Here, we model the probability of the start of an insolvency case as a function of the recent company growth. We use a unique administrative dataset from the German Forschungsdatenzentren der Statistischen Ämter des Bundes und der Länder \citepalias{rdc:2019} which contains over one million companies in the year 2018 and 2019. Correcting for potential endogeneity yields substantively different results, which are robust to controls and are in line with other results in the literature.

The remainder of this paper is organised as follows. Section (ref) introduces the model, the (limited) maximum likelihood estimator, and discusses identification. We lay out assumptions, discuss consistency and asymptotic normality of the estimator as well as inferential methods in Section (ref). Section (ref) contains a Monte Carlo study, while the empirical application is presented in Section (ref).

Model, Identification, and Estimator

Following rivers:88, we consider the following structural model

where \(Z\) is a \(k \times 1\) vector of exogenous regressors (including a constant) and \(D\) is a scalar endogenous regressor correlated with the error term \(U\). Let us further assume that the endogenous variate and the error can be decomposed as

where \(V\) (continuous) and \(E\) are mean-zero error terms independent of $Z$ and $Z$ and $D$, respectively. The functions $\pi(\cdot)$ and $m(\cdot)$ are in general unknown. Thus, unless $\rho=0$, $D$ is endogenous due to the presence of the term $m(V)$.

One object of interest could be for example the average structural function (ASF). Assuming $E/\sigma \sim F$, for some symmetric cumulative distribution function (cdf) $F$, we get \( \textnormal{\textsf{E}}[Y \mid Z,D,V]= F(\sigma^{-1}(\alpha^{\textnormal{\textsf{T}}} Z+\beta D+\rho m(V))), \) because $E$ is independent of $Z$, $D$, and thus also of $V = D-\pi(Z)$. For example, if $F$ is the standard normal cdf $\Phi$, say, and $m(V) \sim \Phi$, then the ASF given by

To make these objects operational, an estimator of the unknown parameters is needed.

Our estimator is based on the following fundamental identification assumption.

If the marginal distribution $G$ of the reduced-form error term $V$ were known and Assumption (ref) ($a$) would hold, then $m(V) = H^{-1}(G(V)) \eqqcolon \eta$ almost surely (see the proof of bmw:24), and we could (up to $\sigma$) identify $\theta = (\gamma^{\textnormal{\textsf{T}}},\rho)^{\textnormal{\textsf{T}}}$, $\gamma = (\alpha^{\textnormal{\textsf{T}}},\beta)^{\textnormal{\textsf{T}}}$, by augmenting the model using the so-called {\it `control function'} $\eta$, i.e.\ $Y = 1\{\alpha^{\textnormal{\textsf{T}}} Z + \beta D + \rho\eta + E\geq0\}$. This is essentially the identification strategy of the popular copula-based endogeneity correction of guptapark:2012 and related approaches; see parkgupta:2024 for a recent review. \footnote{ Monotonicity of \( m \), together with that of either \( H \) or \( G \), implies strict monotonicity of the triple \( (m, H, G) \), which is essential to pointwise identify the control function via \( m(V) = H^{-1}(G(V)) \), and, in conjunction with Assumption (ref) (a), to ensure \( G \neq H \). If any of the functions \( m \), \( G \), or \( H \) fail to be strictly monotonic, then the identity \( m(V) = H^{-1}(G(V)) \) may hold only in distribution, not almost surely. Moreover, if \( G = H \), then strict monotonicity implies \( G(v) = H(v)={\sf P}(V \leq m^{-1}(v))=G(m^{-1}(v)) \), and by injectivity of \( G \), it follows that \( m = {\sf id} \), contradicting Assumption (ref) (a).} If, however, $G=H$, then part ($a$) is violated and identification breaks down as $(Z,D,\eta) = (Z, \pi(Z)+V, V)$ are perfectly collinear unless $\pi(\cdot)$ is nonlinear. In other words, non-linearity of the first stage provides an additional source of identification and yields a `{\it robustification}' against violations of part ($a$) (i.e. $G=H$). This aspect of our identification strategy is novel relative to the earlier work cited above and applies also to the linear models considered in bmw:24 and thus also to several specifications derived from the seminal approach of guptapark:2012.

In practice, $G$ is typically unknown and has to be estimated from a sample $\mathcal{S}_n \coloneqq \{X_i,Y_i\}_{i=1}^n$, say, which is assumed to be an {\sf IID} sample independently drawn from $X\coloneqq (Z^{\textnormal{\textsf{T}}},D)^{\textnormal{\textsf{T}}}$ and $Y$. More specifically, we could estimate $G$ in a first step nonparametrically and construct

Note that $\tilde G_n(V_i)$ is just the relative empirical rank, i.e.\ the rank of $V_i$ among $\{V_1,\dots,V_n\}$ divided by $n+1$. Since also $\pi(\cdot)$ is typically unknown, we could estimate the regression function in a preliminary step to obtain the residuals $V_{i,n} = D_i - \pi_n(Z_i)$ for some estimator $\pi_n(\cdot)$ constructed from $\mathcal{S}_n$ so that

Our instrument-free approach comes at the cost of having to specify $H$, the cdf of the source of endogeneity $m(V)$. In principle, different choices are possible. Here, we follow the literature (see parkgupta:2024 and the references therein), and make the following normality assumption as this allows us to leverage corresponding theoretical results for normal scores (zhao:2020,zhao:2022).

Finally, in order to derive our estimator a link function, i.e.\ the cdf of the innovation $E$ has to be fixed. The following conditions restrict the link function and impose restrictions on $V$ and $E$.

Part ($i$) is a common assumption in the literature and identical to amemiya:1985. Popular choices for $F$, that satisfy part ($i$), include the logistic (i.e.\ $F(z) = \Lambda(z) \coloneqq {\sf exp}(z)/(1+{\sf exp}(z))$) or the standard normal distribution (i.e.\ $F = \Phi$), giving rise to a logit and a probit specification, respectively. Part ($ii$) implies that $E$ is independent of $Z$, $D$, and $V$, while part ($iii$) requires independence between $V$ and $Z$.

Now, set $X \coloneqq (Z^{\textnormal{\textsf{T}}},D)^{\textnormal{\textsf{T}}}$ and define the log-likelihood contribution

where we use the normalization $\sigma = 1$. Our limited information maximum likelihood estimator $\theta_n$, say, of the true parameter vector $\theta_0$ is an extremum estimator defined as a solution (if it exists) of $$ \frac{\partial}{\partial \theta}{\cal L}_n(\theta) = 0, \quad {\cal L}_n(\theta) \coloneqq \frac1{n}\sum_{i=1}^n\ell(\theta;Y_i,X_i,\eta_{i,n}). $$ This is similar to rivers:88 in that we use a control function to cope with endogeneity in a binary response model, with the distinction, however, that our approach does not require outside instruments.

Asymptotic Theory and Inference

Consistency

A crucial tool to derive the asymptotic properties of $\theta_n$ is the following high-level assumption on the difference between the infeasible control function $\tilde\eta_{i,n}$ and its feasible counterpart $\eta_{i,n}$, both defined in Eqs. (ref) and (ref), respectively.

If the regression function is linear $\pi(z) = z^{\textnormal{\textsf{T}}}\delta$, and the $k \times 1$ vector $\delta$ is estimated by the OLS estimator $\delta_n$, say, then as discussed in bmw:24, it follows from zhao:2020 that Assumption (ref) is satisfied for \[ \eta_{i,n}-\tilde\eta_{i,n} = -\kappa(U_i)(\delta_n-\delta)^{\textnormal{\textsf{T}}}(Z_i-\textnormal{\textsf{E}}[Z]) + R_{i,n}, \quad \kappa(u) = \frac{g(G^{-1}(u))}{\phi(\Phi^{-1}(u))}, \] with $\operatorname*{\textnormal{\textsf{max}}}_i |R_{i,n}| = o_p(n^{-1/2})$. If, in addition, $\textnormal{\textsf{E}}[\Vert \kappa(U)Z\Vert^2] < \infty$, then also Assumption (ref) ($ii$) will be satisfied. If $\pi(\cdot)$ is a smooth function estimated nonparametrically using common kernel-based estimation techniques, then Assumption (ref) holds for $\kappa(\cdot)$ defined above under standard regularity conditions as demonstrated in zhao:2022; for part ($ii$) see, e.g. hardle:88.

To derive consistency, we show that ${\cal L}_n(\theta)$ is, uniformly in $\theta \in \Theta \subseteq \mathbb{R}^{k+2}$, close to the infeasible objective function ${\cal L}_{0,n}(\theta) = \frac{1}{n}\sum_{i=1}^{n} \ell(\theta;Y_i, X_i,\eta_i)$. To do so, we have to impose a mild regularity condition to control small perturbations in a neighbourhood around ${\cal L}_{0,n}(\theta)$. More specifically, define \[ \psi(\theta;Y,X,\eta) \coloneqq \frac{Y-F(X^{\textnormal{\textsf{T}}}\gamma + \rho \eta)}{F(X^{\textnormal{\textsf{T}}}\gamma + \rho \eta)(1-F(X^{\textnormal{\textsf{T}}}\gamma + \rho \eta))}f(X^{\textnormal{\textsf{T}}}\gamma + \rho \eta), \] and note that $\rho\psi(\cdot;\cdot,\cdot,t) = \partial \ell(\cdot; \cdot,\cdot,t)/\partial t$. For example, in the probit case ($F=\Phi$), we obtain \[ \psi(\theta;Y,X,\eta) = \frac{Y-\Phi(X^{\textnormal{\textsf{T}}}\gamma + \rho \eta)\phi(X^{\textnormal{\textsf{T}}}\gamma + \rho \eta)}{(1-\Phi(X^{\textnormal{\textsf{T}}}\gamma + \rho \eta))\Phi(X^{\textnormal{\textsf{T}}}\gamma + \rho \eta)}, \] while for a logit specification ($F = \Lambda$), we get $\psi(\theta;Y,X,\eta) = (Y-\Lambda(X^{\textnormal{\textsf{T}}}\gamma + \rho \eta))$.

Following the discussion in baing:2008, Assumption (ref) can be shown to be satisfied for logit or probit specifications of the link function.

Limiting distribution

Having established weak consistency of $\theta_n$, we now turn to the question of its asymptotic distribution. In a first step, the asymptotic normality of the score vector evaluated at the true value $\theta_0$ is established. To achieve this, we assume, similar to rothe:09, that the feasible score function satisfies a linear representation. In particular, define the $(k + 2) \times 1$ score vector $s(t;\cdot,\cdot,\cdot) \coloneqq \partial \ell(t;\cdot, \cdot,\cdot)/\partial t$, where

Let us also define the first derivative of $\psi$ with respect to the last argument $\rho\dot{\psi}(\cdot;\cdot,\cdot,t) =\partial \psi(\cdot; \cdot,\cdot,t)/\partial t$ given by

and set $S(\theta;Y,X,\eta) \coloneqq W\dot\psi(\theta;Y,X,\eta)$, $S_0(Y,X,\eta) \coloneqq W \dot\psi_0(Y,X,\eta),$ with $\dot\psi_0(Y,X,\eta) \coloneqq \dot\psi(\theta_0;Y,X,\eta).$\footnote{To illustrate, note that for a probit specification we get

where $\lambda(x) \coloneqq \phi(x)/\Phi^{-1}(x)$,while for a logit model \[ \dot{\psi}(\theta;Y,X,\eta) = -\Lambda(X^{\textnormal{\textsf{T}}}\gamma + \rho \eta)(1- \Lambda(X^{\textnormal{\textsf{T}}}\gamma + \rho \eta)). \]} We then use the following linear representation:

This assumption warrants some discussion. First, suppose that $\psi(Z) = Z^{\textnormal{\textsf{T}}}\delta$, and $\delta$ is estimated by the OLS estimator $\delta_n$, then $\pi_n(Z_i)-\mathbb{E}_n[\pi_n(Z)]-\pi(Z_i)+\textnormal{\textsf{E}}[\pi(Z)] = (Z_i-\textnormal{\textsf{E}}[Z])^{\textnormal{\textsf{T}}}(\delta_n-\delta_0)$ and Assumption (ref) is satisfied with \( h(Z_i) \coloneqq \textnormal{\textsf{cov}}[S_0(Y,X,\eta)\kappa(U),Z](\textnormal{\textsf{E}}[ZZ^{\textnormal{\textsf{T}}}])^{-1}(Z_i-\textnormal{\textsf{E}}[Z]). \) If, on the other hand, $\pi_n(Z)$ is e.g. the local linear kernel estimator, then Assumption (ref) holds with \( h(Z_i) \coloneqq \textnormal{\textsf{E}}[S_0(Y_i,X_i,\eta_i)\kappa(U_i) \mid Z_i]-\textnormal{\textsf{E}}[S_0(Y,X,\eta)\kappa(U)]. \) To see this, note that the left-hand side of the expression in Assumption (ref) can be written as

where \( \mathbb{G}_n(r_n(\cdot)) \coloneqq \frac1{\sqrt{n}}\sum_{i=1}^n(r_n(X_i)-\mathbb{E}_n[r_n(X)]),\) for some function $r_n$ that might depend on ${\cal S}_n$. By stochastic equicontinuity (for primitive sufficient conditions see escanciano2014uniform), $\mathbb{G}_n(S_0(\cdot)\kappa(\cdot) \pi_n(\cdot))-\mathbb{G}_n(S_0(\cdot) \kappa(\cdot)\pi(\cdot)) = o_p(1)$, while \[ \frac1{\sqrt{n}} \sum_{i=1}^n \mathbb{E}_n[S_0(Y_i,X_i,U_i)\kappa(U_i)(\pi_n(Z)-\pi(Z))] = \frac1{\sqrt{n}} \sum_{i=1}^n \textnormal{\textsf{E}}[S_0(Y,X,U)\kappa(U)]V_i + o_p(1), \] because, by the LLN, $n^{-1}\sum_{i=1}^nS_0(Y_i,X_i,U_i)\kappa(U_i) = \textnormal{\textsf{E}}[S_0(Y,X,U)\kappa(U)]+o_p(1)$ and \( \mathbb{E}_n[\pi_n(Z)-\pi(Z)] = \int (\pi_n-\pi)(z)f_Z(z)dz = \frac1{n}\sum_{j=1}^n V_j + o_p(1). \) Thus, Assumption (ref) holds.

Note that $\Omega_2 = 0$ if $G(\cdot)$ is known, $\Omega_3 = 0$ if $\pi(\cdot)$ is known, and $A = \Omega_1$ if $D$ is exogenous.

In order to derive the limiting distribution of $\theta_n$, we have to make sure that the Hessian \[ H(\theta;Y,X,\eta) \coloneqq \frac{\partial}{\partial \theta} s(\theta;Y,X,\eta) = -WW'\dot\psi(\theta;Y,X,\eta) \] obeys a uniform LLN. The following assumption is similar to baing:2008 and, using their arguments, can be shown to hold for logit and probit specifications.

Inference

As Proposition (ref) reveals, the limiting distribution depends on unknown nuisance parameters. An important special case concerns hypotheses that contain the restriction of no endogeneity (i.e.\ $\rho = 0$). In this case we can use common textbook standard errors, as $\sqrt{n}(\theta_n-\theta_0) \rightarrow_d \mathcal{N}(0,\Omega_1^{-1})$, where $\Omega_1$ is consistently estimated using standard approaches (amemiya:1985).

For more general hypotheses, we propose the following pairs bootstrap: Draw $(Y_{b,1},X^{\textnormal{\textsf{T}}}_{b,1})^{\textnormal{\textsf{T}}}$,$\dots$,$(Y_{b,n},X^{\textnormal{\textsf{T}}}_{b,n})^{\textnormal{\textsf{T}}}$ with replacement from the empirical distribution of the original data ${\cal S}_n$ and define, analogously to \(\theta_n\), \(\theta_{n,b}\) based on the bootstrap data. We can then construct bootstrap standard errors via $\Sigma_{n,B} \coloneqq \frac{n}{B}\sum_{b=1}^B (\theta_{n,b}-\theta_n)(\theta_{n,b}-\theta_n)^{\textnormal{\textsf{T}}}$. Following the discussion surrounding bmw:24, consistency of $\Sigma_{n,B}$ conditionally on the original data ${\mathcal S}_n$ (as $n$ and $B$ diverge) follows if we assume that $\sqrt{n}(\theta_n-\theta_0)$ possesses uniformly integrable second moments.

If interest lies in other functionals, for example, the ASF introduced in Eq. (ref), one could apply the delta-method in conjunction with $\Sigma_{n,B}$ (woold:2010). To fix ideas, consider the probit-specification (i.e. $F = \Phi$), then we can estimate the ASF via \[ {\sf ASF}_{n}(x) \coloneqq \Phi\left(\frac{\theta_{n,\alpha}^{\textnormal{\textsf{T}}} z+ \theta_{n,\beta} d}{\sqrt {1+\theta_{n,\rho}^2}}\right) \] for some $x = (z^{\textnormal{\textsf{T}}}, d)^{\textnormal{\textsf{T}}} \in \mathbb{R}^{k+1}$ and the partition $\theta_n = (\theta_{n,\alpha}^{\textnormal{\textsf{T}}},\theta_{n,\beta},\theta_{n,\rho})^{\textnormal{\textsf{T}}}$ corresponds to $W$. An estimator of the asymptotic variance is then given the sandwich form $(\nabla_\theta {\sf ASF}(x)\vert_{\theta = \theta_n})^{\textnormal{\textsf{T}}} \Sigma_{n,B} \nabla_\theta {\sf ASF}(x)\vert_{\theta = \theta_n}$, where $\nabla_\theta {\sf ASF}(x)$ is the $(k+2) \times 1$ gradient vector.

Relaxing Assumptions 2 and 3

{\it Unknown distribution of the endogeneity.} Instead of requiring that $m(V) \sim H=H_0$ is known, one could assume that the unknown $H_0$ belongs to a class of parametric distributions ${\mathcal H}$, say, for which $H_0 \coloneqq H(\lambda)$ is known up to a finite dimensional parameter $\lambda = \lambda_0 \in {\sf int}(\Lambda)$, $\Lambda \in \mathbb{R}^l$. Our estimator of $(\theta_0^{\textnormal{\textsf{T}}},\lambda_0^{\textnormal{\textsf{T}}})^{\textnormal{\textsf{T}}}$ then maximizes the following objective function:

where $\eta_{i,n}(\lambda) \coloneqq H^{-1}(G_n(V_{i,n});\lambda)$. Given that $H^{-1}(\lambda)$ enters the objective function, it would be desirable to specify $H(\lambda)$ such that the inverses have closed-form representations. One such computationally appealing yet flexible choice of ${\cal H}$ could be the class of asymmetric distributions considered by gijbels:2019. A special case is the two-piece skew-normal distribution of mudhut:00 for which

Adopting this specification, Assumption (ref) can be empirically tested via $H_0$: $\lambda = 0$. Although, in principle, the properties of the estimator could be investigated by leveraging the results developed here and the likelihood theory obtained by gijbels:2019, a theoretical treatment is well beyond the scope of the current paper.

{\it Unknown distribution of the innovation.} In case the true link function $F=F_0$ is unknown, one could replace $F_0$ in Eq. (ref) with a nonparametric estimate so that our estimator $\theta^\ddagger_n$, say, maximizes

where $\tau_i \coloneqq 1\{ (X_i,\eta_{i,n}) \in {\cal X} \}$ is a trimming function for a compact set ${\cal X}$ and, for a given $\theta$, $F_{n}$ is the Nadaraya-Watson estimator of $F_0$, i.e. a nonparametric kernel regression of $Y_i$ on $X_i^{\textnormal{\textsf{T}}}\gamma + \rho\eta_{i,n}$, $i \in \{1,\dots,n\}$. The dependence of $F_n$ on $\theta$ is implicitly understood for ease of notation. This is idea is similar to the approach proposed initially by klein:93 and then further extended by, among others, blundell:2004 and rothe:09. Adapting the arguments of rothe:09, it might be possible to show that $\theta^\ddagger_n$ is $\sqrt{n}$-consistent with asymptomatic Gaussian limiting distribution; a conjecture supported by the finite sample evidence of the following section. Finally, a nonparametric estimator of the ASF in Eq. (ref) can be obtained via \[ \widetilde{\sf ASF}_n(x) \coloneqq\frac1{n}\sum_{i=1}^n F_n(x^{\textnormal{\textsf{T}}} \theta^\dagger_{n,\gamma} + \theta^\dagger_{n,\rho}\eta_{i,n}), \quad x \in \mathbb{R}^{k+1}. \]

Monte Carlo Simulation

The Monte Carlo design is based on Eq. (ref), i.e.\ $Y = 1\{\alpha_0 + \alpha_1 Z + \beta D+ U>0\}$, where $Z \sim \Phi$ and $D=\pi(Z)+V$. We consider a probit specification, where $E \sim \Phi$ in $U = \rho m(V)+E$. The assumption that $m(V) \sim \Phi$ is maintained throughout, while the distribution $G$ of the reduced form error $V$ is $G = \Phi$ or $G = {\sf Gamma}(2,2)$, respectively. We distinguish between a linear ($\pi(z) = z$) and a non-linear ($\pi(z) = z^2$) specification of the reduced form. It is apparent from the discussion of Assumption (ref) that identification breaks down if $G$ is normal and $\pi(\cdot)$ is linear.

{1.4pt}

{1.4pt}

{1.4pt}

{1.4pt}

In each of the 1,000 Monte Carlo repetitions, we take {\sf IID} draws $\{Y_i, Z_i, D_i\}_{i=1}^n$, where $n \in \{500,$ 1,000$\}$, based on the specifications discussed earlier for estimation and inference. We consider eight different estimators:

Finally, we consider also the three last estimators (3)-(5) with the link function estimated using the Nadaraya-Watson estimator. These three additional estimators are denoted by (6) {\it np}{\sf MW}1, (7) {\it np}{\sf MW}2, and (8) {\it np}{\sf DONG}. The kernel-based estimation of the link function and the ASF (see Section 3.4.) is implemented using the {\sf R} function {\sf kreg} with default settings.

For estimators (1)-(5), we use the normalization $\sigma = 1$, while, following the literature, for (6)-(8), we use $\alpha_1 = 1$. Note that, due to the local level specification of the nonparametric estimator of the link function, estimators (6)-(8) do not include a constant.

We report the following metrics for each estimator: the mean, the standard deviation, the root-mean-squared error, the empirical size of a two-sided $t$-test at the nominal significance level of 5% for the estimators of $\theta$ and the ASF evaluated at the mean of $X$. For estimators (3)-(5), we compute test statistics using bootstrap standard errors with 499 repetitions. For computational reasons, estimators (6)-(8) use 99 bootstrap repetitions.

Table (ref) shows that in case of endogeneity ($\rho = 0.5$), the proposed endogeneity correction does its job as long as Assumption (ref) is satisfied. The naïve probit estimator ({\sf ML}) displays severe bias and size distortions unless endogeneity is absent (i.e.\ $\rho = 0.0$, see Table (ref)), in which case it coincides with {\sf CF}0. As expected, the case of a linear reduced form in conjunction with $H=G$ (i.e.\ violation of Assumption (ref)) leads to non-identification due to collinearity. Evidently, in this case {\sf CF}0 is not defined as collinearity becomes perfect. We note that {\sf MW}2 has severe problems in case the nonlinear first-stage is misspecified, while the efficiency loss of {\sf MW}1 that estimates the first stage nonparametrically relative to its infeasible counterpart {\sf CF}0 seems acceptable. When the link function is estimated nonparametrically, estimation precision--though still satisfactory--is generally lower compared to the probit specifications. As shown in Tables (ref) and (ref), both estimation performance and size control improve as the sample size increases from $n = 500$ to $n = $1,000.

Finally, the estimator proposed by dong:2010, which uses the nonparametrically estimated innovation $V$ as a control function, is more efficient than {\sf MW}1 in the setting where $\pi(z) = z^2$ and $V \sim \Phi$, so that $m(V) = V$. This is because the rank-based estimation of $m(\cdot)$, as employed by {\sf MW}1, is superfluous here. In all other scenarios, in particular if the dependence in the first step is linear, {\sf MW}1 outperforms {\sf DONG}.

Application to Insolvency Risk

We consider German administrative data from the Forschungsdatenzentren der statistischen Ämter des Bundes und der Länder, which contains all German companies (Rechtseinheiten) in the year 2018 and 2019. The general task is to model insolvency risk, which is a currently relevant topic, see e.g. weissbachwied:2022. In particular, we are interested in measuring the influence of company growth on insolvency risk: Is strong growth an indicator for a healthy company or does strong growth imply substantial risk? Dependencies between company growth and insolvency risk have been of interest in corporate development for a long time, see bensoussan:1981, santanna:2017 or xuezhouetal:2022.

With this question, potential endogeneity issues arise: There might be reverse causality (if insolvency lies on the table, employees might leave the company) or the existence of a latent relevant variable which measures the current quality of the management and related aspects.

The dependent variable is the indicator variable if an insolvency case starts in 2019. While such insolvency cases can take several years, typically a five-digit number of companies actually becomes insolvent in Germany per year (weissbachwied:2022). The base model is given by

Potential endogenous variables are the sales growth from 2018 to 2019 and the employee growth from 2018 to 2019. Exogenous controls are the sales in 2018, the employees in 2018 the German state and the legal status (Rechtsform). Companies whose sales in 2018 are below and above the 5% and 95% quantiles of these sales are excluded. This leads to a sample size of $n=1$,$131$,$230$. In the sample, for 3,412 companies, an insolvency case starts in 2019.

We believe that it is reasonable to assume the existence of normally distributed latent terms, which determine the sales and employee growths. These terms refer to “management intelligence” and there is evidence for the fact that such intelligence-related terms are normally distributed (bmw:24). On the other hand, growth values are typically non-normally distributed, so that our nonlinearity condition should be fulfilled. In our dataset, there is no explicit information about such terms. Moreover, no instruments such as external or internal firm growth as considered in xuezhouetal:2022 are available.

We show in Table (ref) the results for a probit regression without control terms ({\sf ML}), for control terms (one for both endogenous variables) with a nonparametric first step ({\sf MW}1) and a linear first step ({\sf MW}2). Moreover, we provide a comparison with the Probit estimator from dong:2010. The estimates for the state and legal status are omitted for brevity. The standard errors for the regressions with control terms are obtained via bootstrap (99 replications), the other ones via the standard Fisher information from the likelihood.

The results indicate that the control terms are very relevant. Without including them, the estimates for $\beta_1$ and $\beta_2$ are statistically significantly negative. With the terms, they become positive in most cases, whereas the $t$-statistics decrease in absolute values. Notably, the estimates for the control terms are statistically significantly negative.

This supports the interpretation that the control terms capture the current quality of management and related factors. A higher management quality reduces the likelihood of insolvency proceedings. When this factor is accounted for, an increase in sales and employment raises the probability of insolvency, likely because firms take on greater risks in pursuit of growth.

Interestingly, the standard errors of the dong:2010 estimates are substantially larger than those from our approach. This suggests that the dependence structure in the first step of our model is more linear than nonlinear. This conclusion is further supported by the fact that our estimator's results remain consistent between the linear and nonparametric specifications in the first step. More evidence is provided by additional plots of the fitted first stage values, which are not reported.

Similar results were obtained in xuezhouetal:2022. These authors consider a partly similar model, but use a different estimation approach. In the same spirit as our analysis, their estimates for firm growth increase, once “mediation variables” (with negative coefficient estimates) are included into the model.

Summary and Outlook

This paper addresses a gap in the literature by proposing a rank-based endogeneity correction for binary outcome models in the presence of endogeneity, without relying on external instruments. The approach allows for both linear and nonlinear dependence between endogenous and exogenous regressors. While we focus on the case of a known link function, we also discuss potential extensions, including methods for handling unknown link functions, drawing inspiration from klein:93, blundell:2004, and rothe2010nonparametric. Another promising direction for future research is the extension to distribution regression models, which build on binary outcome models as discussed recently by wied:2024.

Declarations

Data Availability Statement

We use a unique dataset from the Research Data Centres of the Federal Statistical Office and Statistical Offices of the Federal States of Germany, which is not publicly available. Fee-based access can be granted by signing a contract.

Funding and Competing Interests

All authors certify that they have no affiliations with or involvement in any organization or entity with any financial interest or non-financial interest in the subject matter or materials discussed in this manuscript. The authors have no funding to report.

Proofs

\setcounter{equation}{0} {\bf Proof of Proposition (ref).} The claim follows if we can show that the objective function ${\cal L}_n(\theta)$ is uniformly close to the infeasible objective ${\cal L}_{0,n}(\theta)$ that uses the unknown control function $\eta_i$ in place of $\eta_{i,n}$. To that end, let $\bar \eta_{i,n}$ (random) be on the line segment connecting $\eta_{i,n}$ and $\eta_i$. Then, by the mean-value theorem and Cauchy-Schwarz, we get

As argued in the treatment of their term `$B$' in the proof of zhao:2020 (which builds upon zhao:19), we get $\sum_{i=1}^{n}(\tilde\eta_{i,n}-\eta_i)^2 = o_p(n)$. Hence, as by Assumption (ref) $\sum_{i=1}^{n}(\tilde\eta_{i,n}-\eta_{i,n})^2 = o_p(n)$, it follows from the triangle inequality for the last term on the right-hand side of Eq. (ref), $\sum_{i=1}^{n}(\eta_{i,n}-\eta_i)^2 = o_p(n).$ Moreover, we obtain from Assumption (ref) \[ \operatorname*{\textnormal{\textsf{sup}}}\limits_{\theta \in \Theta}\frac{1}{n}\sum_{i=1}^{n} \psi^2(\theta;Y_i,X_i,\bar\eta_{i,n}) \leq \operatorname*{\textnormal{\textsf{sup}}}\limits_{\theta \in \Theta}\frac{1}{n}\sum_{i=1}^{n}\operatorname*{\textnormal{\textsf{sup}}}_{t_i: |t_i-\eta_i|\leq b_n}\psi^2(\theta;Y_i,X_i,t_i)= O_p(1). \] This shows $\operatorname*{\textnormal{\textsf{sup}}}\limits_{\theta \in \Theta}|{\cal L}_n(\theta) - {\cal L}_{0,n}(\theta)| = o_p(1)$ and the claim follows. $\square$\\

{\bf Proof of Proposition (ref)}. Define $A_1 \coloneqq n^{-1/2}\sum_{i=1}^ns_0(Y_i,X_i,\eta_i)$ and note that $A_1 \rightarrow_d {\cal A} _1 =_d \mathcal{N}(0,\Omega_1)$. Next, consider

say. Begin with $A$ and note that, by a first-order Taylor expansion, we obtain

with $A_2$ and $A_3$ being implicitly defined. Begin with $A_2$ and note that $\eta_i = \Phi^{-1}(U_i)$, with $U_i = G(V_i)$ being an {\sf IID} sequence of ${\sf Unif}[0,1]$ variates. Let $K_n$ denote the empirical cdf of $U_i$. Then, there exists an ordering of the indices $\{1,\dots ,n\}$, such that, by a first-order Taylor expansion, we obtain for $-A_2$

almost surely, for a standard Brownian Bridge $B(\cdot)$ so that $\textnormal{\textsf{var}}[{\cal A}_2] = \Omega_2$. Here we used a similar argument to the proof of bmw:24. Finally, $A_3 \rightarrow_d {\cal A}_3$, ${\cal A}_3 \sim \mathcal{N}(0,\Omega_3)$ follows by Assumptions (ref) and (ref). The claim thus follows because, by Cauchy-Schwarz and the previous result, $B$ and $C$ are $o_p(1)$ and $\operatorname*{\textnormal{\textsf{lim}}}\textnormal{\textsf{cov}}[A_j,A_i] = 0$, $i\neq j$. $\square$\\

{\bf Proof of Proposition (ref)}. Since $\theta_n$ is a solution of the maximisation problem ${\cal L}_n(\theta)$ it follows that $\sum_{i=1}^ns(\theta_n;Y_i,X_i,\eta_{i,n})=0_{k+2}$ and, by a first order Taylor expansion about $\theta_0$, we get $$\sqrt{n}(\theta_n-\theta_0) = \left[-\frac1{n}\sum_{i=1}^nH(\theta_0;Y_i,X_i,\eta_{i})+O_p(\vert \theta_0-\theta_n\vert^2)\right]^{-1}\frac1{\sqrt n}\sum_{i=1}^ns(\theta_0;Y_i,X_i, \eta_{i,n}),$$ where the remainder term is due to Assumption (ref). The claim then follows from standard arguments. $\square$\\

\singlespacing \addcontentsline{toc}{section}{References}