# Another Look at the Linear Probability Model and Nonlinear Index Models

Kaicheng Chen Department of Economics, Michigan State University. Email:
<EMAIL_ADDRESS>Robert S. Martin Division of Price and Index Number Research,
Bureau of Labor Statistics. Email<EMAIL_ADDRESS>Jeffrey M. Wooldridge
Department of Economics, Michigan State University<EMAIL_ADDRESS>
###### Abstract

We reassess the use of linear models to approximate response probabilities of
binary outcomes, focusing on average partial effects (APE). We confirm that
linear projection parameters coincide with APEs in certain scenarios. Through
simulations, we identify other cases where OLS does or does not approximate
APEs and find that having large fraction of fitted values in $[0,1]$ is
neither necessary nor sufficient. We also show nonlinear least squares
estimation of the ramp model is consistent and asymptotically normal and is
equivalent to using OLS on an iteratively trimmed sample to reduce bias. Our
findings offer practical guidance for empirical research.

Keywords: Binary response; linear probability model; average partial effect;
nonlinear least square; probit model.

JEL Classification Code: C25

##  1 Introduction

When an outcome variable, $y$, is binary, empirical researchers usually choose
between two general strategies given a vector of (exogenous) explanatory
variables, $\mathbf{x}$: (i) approximate the response probability,
$P\left(y=1|\mathbf{x}\right)$, using a model linear in parameters or (ii) use
a nonlinear model, such as logit or probit. The first strategy is commonly
known as using a linear probability model (LPM). The benefits of the LPM are
well-known and include ease of interpretation, simple estimation, and
straightforward extension to situations with endogenous explanatory variables
(so that instrumental variables are used) and panel data settings with
unobserved heterogeneity. The shortcomings of the LPM are also well known, and
discussed in most introductory econometrics texts; see, for example,
Wooldridge,, 2019, Section 7.5. More advanced discussions of the LPM recognize
that one should not take the linear model for $P\left(y=1|\mathbf{x}\right)$
literally but only as an approximation. The approximation can be exact in
special cases—such as when $\mathbf{x}$ consists of binary indicators that are
exhaustive and mutually exclusive—and it may be poor in other cases. However,
for the most part, prediction is not the primary use of LPMs specifically or
binary response models generally. Rather, researchers are largely interested
in using binary response models to measure ceteris paribus or causal effects,
and it is from this perspective that the LPM approximation should be
evaluated. Angrist and Pischke,, 2009, Section 3.4.1 and Wooldridge,, 2010,
Section 15.2 take this perspective. Wooldridge,, 2010, Section 15.6, p. 579
shows how the results of Stoker, (1986) can be applied to OLS estimation of
the parameters in a linear probability model. Remarkably, there are situations
where the linear projection exactly recovers the APEs across a broad range of
binary response models. As is well known—see, for example, Wooldridge,, 2010,
Chapter 4—under standard sampling assumptions—OLS consistently estimates the
parameters of the linear projection (LP).

Even though it is natural to study the LPM from the linear projection
perspective, this opinion is not universally held. In an influential paper,
Horrace and Oaxaca, (2006) study both the bias and inconsistency of the OLS
estimator for the parameters of an underlying piecewise linear model for the
response probability that ensures the probabilities are in the unit
interval.111Horrace and Oaxaca, (2006) defines the LPM as the piecewise linear
ramp model. However, in this paper we differentiate between the “ramp model”
and the “LPM” (which is linear everywhere). The Horrace-Oaxaca (H-O) paper is
regularly cited in empirical research, sometimes as a cautionary tale in using
the LPM and sometimes as support for using the LPM when relatively few fitted
values lie outside the unit interval. (In the previous two years, H-O has
almost 200 Google Scholar citations.) While H-O take the piecewise linear
model seriously, much if not most of the citing literature seeks to use their
results to choose between the LPM and an alternative like probit or
logit.222See, for example, Footnote 20 of van den Berg and Siflinger, (2022).

In the current paper, we revisit the H-O framework but, rather than focus on
parameters, we focus on APEs. We show that H-O set up the problem so that, in
general, the response probability is nonlinear in the underlying linear index,
$\mathbf{x\beta}=\beta_{1}+\beta_{2}x_{2}+\cdots+\beta_{K}x_{K}$. The
nonlinear function of $\mathbf{x\beta}$—sometimes called the ramp function—is
piecewise linear and continuous, but it is not strictly increasing, and it is
nondifferentiable at two inflection points. Nevertheless, under fairly weak
assumptions, one can define the average partial effects, and these are
necessarily smaller in magnitude than the corresponding parameter in the
underlying nonlinear model. Consequently, H-O’s focus on parameters rather
than APEs is essentially the same as focusing on parameters in smooth response
probabilities such as the logit and probit functions. Therefore, any
conclusions about the usefulness of the LPM should be reexamined from the
perspective of identifying APEs rather than coefficients.

It is important to understand that we are not necessarily advocating the H-O
ramp function as an especially sensible model of the response probability.
Rather, we primarily study that specification from the perspective of average
partial effects to determine how the OLS estimator holds up. Briefly, in some
cases OLS does a very good job of approximating the APEs even when a larger
percentage of the fitted values are outside the unit interval. Conversely, in
other cases, OLS does a very poor job of approximating the APEs even when a
high percentage of the fitted values are within the unit interval. A practical
implication is that there is little justification for how the H-O study is
cited in empirical research.

We compare OLS to a few nonlinear competitors, including probit and logit
quasi-maximum likelihood estimation (QMLE), as natural benchmarks. H-O cite a
few theoretical rationalizations for the ramp model, so it also makes sense to
see if a consistent estimator exists which takes it seriously. H-O suggest
trimming the sample of fitted values outside the unit interval and re-
estimating using OLS, but do not present any theoretical or simulation
results.333In unreported simulations, we found that trimming the sample once
did not necessarily improve performance over OLS for estimating the APEs. In
Section 4, we show that nonlinear least squares estimation (NLS) using the
ramp function is consistent and asymptotically normal under mild assumptions.
We also give a variance estimator for the asymptotic variance of the NLS
estimator and provide a consistency result for this. We are unaware of any
other studies using NLS to estimate the ramp model. We also examine an
iterative trimming OLS (ITO) procedure similar in spirit to H-O’s suggestion
and show it is equivalent to numerically minimizing the NLS objective function
using the well-known Newton-Raphson algorithm. In simulations, for estimating
the APEs, we find that NLS estimation of the ramp function performs comparably
to quasi-MLE estimation of the logit and probit models and has good finite
sample properties even when OLS estimation of the LPM does not. For
completeness, we also consider a local linear estimation of a nonparametric
model of the conditional mean, but it seems for our data generating processes
there is not much gain in using a nonprarametric model compared to other
nonlinear models.

In Section 2 we present the population model equivalent of the Horrace-Oaxaca
model and show that it is equivalent to a latent variable model with uniformly
distributed errors. We also derive the average partial effects of the so-
called “ramp function.” In Section 3 we extend the discussion in Wooldridge,,
2010, Section 15.6 and show that, when the covariates have a multivariate
normal distribution, the linear projection identifies the APEs. Section 5
contains several simulations comparing sample APEs produced by OLS, NLS,
Probit QMLE, Logit QMLE, and Local Linear estimations when the true response
probability follows the ramp model. We provide some results not covered by the
existing theories and show that a large fraction of fitted values in $[0,1]$
is neither sufficient or necessary condition for OLS to well-approximate the
APEs.

In Section 6, we apply LPM, ramp, probit, and logit models to an empirical
study of discrimination in mortgage lending decisions. The LPM estimated by
OLS estimation, with a full set of interactions between the race indicator and
the control variables, delivers a notably smaller and marginally statistically
significant estimate of the race effect on the approval probability. The NLS,
probit QMLE, and logit QMLE are very similar and all statistically significant
at the 0.2% level—both because the estimated effects are larger but also
because the (robust) standard errors are notably smaller. In Section 7, we
conclude with some implications for empirical research.

##  2 The Ramp Model and the Linear Projection

One of the key features of the Horrace-Oaxaca setting is that it imposes the
logical bounds on the response probability in a setting where, over some of
its range, the response probability is linear in an index. Almost all of the
important arguments are in terms of the underlying population model, so that
is our focus. Let $y$ be the binary outcome variable and $\mathbf{x}$ the
$1\times K$ vector of explanatory vairables, where $x_{1}\equiv 1$ allows for
an intercept in the index. Defining $p(\mathbf{x})=P(y=1|\mathbf{x})$, then
H-O’s specification can be written as
$p\left(\mathbf{x}\right)=R(\mathbf{x\beta})$ where where
$\mathbf{x\beta}=\beta_{1}+\beta_{2}x_{2}+\cdots+\beta_{K}x_{K}$ is the linear
index and

$\displaystyle R\left(z\right)=\begin{cases}0,&\text{$\mathbf{z}\leq 0$}\\\
\mathbf{z},&\text{$\mathbf{z}\in(0,1)$}\\\ 1,&\text{$\mathbf{z}\geq
1$}\end{cases}\vspace{-3mm}$ (2.1)

is the “ramp function”. This response probability was also suggested by
Horowitz and Savin, (2001) as being suitable when one starts with a linear
model for $p\left(\mathbf{x}\right)$ but wants to ensure that the
probabilities are within the unit interval. The ramp function, which is
piecewise linear, is plotted in Figure 1.

Figure 1: The Ramp Function

In defining partial effects we must recognize that the function
$R\left(z\right)$ is nondifferentiable at $z=0$, $z=1$. A very common
assumption imposed in the semiparametric literature on binary response models
is that at least one element of $\mathbf{x}$ is continuous, and that element
has a nonzero coefficient. Then $\mathbf{x\beta}$ is continuous, and so

$P\left(\mathbf{x\beta}=0\right)=P\left(\mathbf{x\beta}=1\right)=0.\vspace{-3mm}$

In what follows, we maintain that $\mathbf{x\beta}$ is continuous so that
partial effects are well-defined with probability one.

Now let $x_{j}$ be a continuously distributed explanatory variable. For
simplicity, the discussion here assumes that $x_{j}$ appears only by itself.
If the model includes quadratics, interactions, and so on then the details
become more complicated but the conclusions do not change substantively.

We can define a partial effect function as the derivative of
$R\left(\mathbf{x\beta}\right)$ and ignore points where the derivative does
not exist:

$PE_{j}\left(\mathbf{x}\right)=\frac{\partial p}{\partial
x_{j}}\left(\mathbf{x}\right)=\beta_{j}1\left[0\leq\mathbf{x\beta}\leq
1\right]\vspace{-3mm},$

where $1\left[\cdot\right]$ is the indicator function. This expression for
$PE_{j}\left(\mathbf{x}\right)$ uses the fact that

$\displaystyle\frac{dR}{dz}\left(z\right)$ $\displaystyle=$ $\displaystyle
1\text{, }z\in\left(0,1\right)$ $\displaystyle=$ $\displaystyle 0\text{ if
}z<0\text{ or }z>1,\vspace{-3mm}$

and is undefined at $z=0$ or $z=1$. In obtaining the average partial effect of
$x_{j}$ we need not define the derivative of $R\left(\cdot\right)$ at the
points $z=0$ and $z=1$ because with a continuously distributed variable in
$\mathbf{x}$ and non-zero coefficient it takes on these two values with
probability zero. Therefore, the APE is

$\alpha_{j}\equiv
E\left[PE_{j}\left(\mathbf{x}\right)\right]=\beta_{j}P\left(0\leq\mathbf{x\beta}\leq
1\right)=\beta_{j}P\left(0<\mathbf{x\beta}\leq
1\right)=\beta_{j}\left[F_{\mathbf{x\beta}}\left(1\right)-F_{\mathbf{x\beta}}\left(0\right)\right],\vspace{-3mm}$
(2.2)

where $F_{\mathbf{x\beta}}\left(\cdot\right)$ is the CDF of $\mathbf{x\beta}$.

There are some simple but useful observations about (2.2). First, $\alpha_{j}$
always has the same sign as $\beta_{j}$. Second, because
$F_{\mathbf{x\beta}}\left(1\right)-F_{\mathbf{x\beta}}\left(0\right)\leq 1$,
$\left|\alpha_{j}\right|\leq\left|\beta_{j}\right|$; with wide support for
$\mathbf{x\beta}$, $\alpha_{j}$ can be much smaller in magnitude than
$\beta_{j}$. Moreover, $\alpha_{j}=\beta_{j}$ if and only if

$P\left(\mathbf{x\beta}\in\left[0,1\right]\right)=1,\vspace{-3mm}$ (2.3)

which means the support of $\mathbf{x\beta}$ is inside the unit interval. This
is essentially the condition used by Horrace and Oaxaca (2006) to conclude
that the OLS estimator in a linear regression is unbiased and consistent for
$\mathbf{\beta}$. Our goal here is to compare the OLS estimators with the APEs
in the general case where
$P\left(\mathbf{x\beta}\in\left[0,1\right]\right)<1$; the H-O condition is
then a special case where the index coefficient, $\beta_{j}$, is identical to
the APE, $\alpha_{j}$.

There is a third set of parameters important for the discussion, and those are
the linear projection parameters. Assume that the $x_{j}$ have finite second
moments and that the $K\times K$ matrix
$E\left(\mathbf{x}^{\prime}\mathbf{x}\right)$ is nonsingular; this simply
rules out perfect collinearity in the population. Then we can always define
the $K\times 1$ vector $\mathbf{\gamma}$ as

$\mathbf{\gamma}=\left[E\left(\mathbf{x}^{\prime}\mathbf{x}\right)\right]^{-1}E\left(\mathbf{x}^{\prime}y\right).\vspace{-3mm}$

We then write the linear projection of $y$ on $\left(1,x_{2},...,x_{K}\right)$
as.

$L\left(y|\mathbf{x}\right)=L\left(y|1,x_{2},...,x_{K}\right)=\gamma_{1}+\gamma_{2}x_{2}+\cdots+\gamma_{K}x_{K}=\mathbf{x\gamma}\vspace{-3mm}.$

In understanding the H-O findings, and their limitations, it is important to
know that $\alpha_{j}$, $\beta_{j}$, and $\gamma_{j}$ are all well-defined
parameters and, in general, they are all different. Defining $\mathbf{\beta}$
and $\mathbf{\alpha}$ requires an underlying model for the response
probability whereas defining $\mathbf{\gamma}$ does not.

As is well known, under random sampling the OLS estimator consistently
estimates the parameters of the linear projection; see, for example,
Wooldridge (2010, Chapter 4.2). In other words, if we run the OLS regression
underlying LPM estimation,

$y_{i}\text{ on }1\text{, }x_{i2}\text{, ..., }x_{iK}\text{,
}i=1,...,N,\vspace{-3mm}$

and obtain the $\hat{\gamma}_{j}$, then
$\hat{\gamma}_{j}\overset{p}{\rightarrow}\gamma_{j}$. Again, this result holds
free of any kind of underlying model.

H-O study the consistency of the $\hat{\gamma}_{j}$ when considered as
estimators of $\beta_{j}$—the coefficients in the index. In other words, their
asymptotic analysis is the same as comparing the linear projection parameters
$\gamma_{j}$ to the index parameters $\beta_{j}$. Our view is that this does
usually not make much sense—for the same reason we do not study consistency of
the OLS estimator for the index parameters in, say, probit or logit. If one
explicitly models the response probability as a nonlinear function of
$\mathbf{x\beta}$ then one must recognize that nonlinearity when defining the
parameters of interest. When interest is in the effects of the explanatory
variables on the response probability—which describes almost all modern usages
of the LPM—it only makes sense to compare the linear projection parameters to
the APEs. In other words, we should ask: When is $\gamma_{j}$ “close” to
$\alpha_{j}$? This is not the same as studying when $\gamma_{j}$ is “close” to
$\beta_{j}$ (except in the special case where (2.3) holds).

Under the H-O ramp model we can write

$E\left(y|\mathbf{x}\right)=p\left(\mathbf{x}\right)=1\left[0\leq\mathbf{x\beta}\leq
1\right]\mathbf{x\beta}+1\left[\mathbf{x\beta}>1\right]\vspace{-3mm}.$

If (2.3) holds then, with probability one,

$E\left(y|\mathbf{x}\right)=\mathbf{x\beta}=L\left(y|\mathbf{x}\right),\vspace{-3mm}$

in which case $\alpha_{j}=\beta_{j}$ and the OLS estimators,
$\hat{\gamma}_{j}$ are consistent for $\beta_{j}$ (which is the APE of
$x_{j}$). If for a random sample of size $N$,
$\mathbf{x}_{i}\mathbf{\beta}\in\left[0,1\right]$ for all $i$, then

$E\left(y_{i}|\mathbf{x}_{1},\mathbf{x}_{2},...,\mathbf{x}_{n}\right)=\mathbf{x}_{i}\mathbf{\beta},\vspace{-3mm}$

and it follows that the OLS estimators are conditionally unbiased for the
$\beta_{j}$ – the conclusion reached in H-O.

If (2.3) fails then $\beta_{j}$ measure the partial effect when
$0\leq\mathbf{x\beta}\leq 1$, but this restriction depends on the unknown
vector $\mathbf{\beta}$. If
$P\left(\mathbf{x\beta}\in\left[0,1\right]\right)<1$ then the $\beta_{j}$ need
not be very useful as summary measures of the partial effects. In the next
section we discuss when the LP parameters identify APEs, with (2.3) being a
special case.

In the next section, it is useful to observe that the response probability in
(2.1) can be derived from a latent variable formulation. Suppose that

$\displaystyle y^{\ast}$ $\displaystyle=\mathbf{x\beta}+u,$ (2.4)
$\displaystyle u|\mathbf{x}$
$\displaystyle\sim\mathrm{Uniform}\left(0,1\right),$ (2.5) $\displaystyle y$
$\displaystyle=1\left[y^{\ast}>0\right]\vspace{-3mm}.$ (2.6)

Because the CDF of $u$ is identical to the ramp function
$R\left(\cdot\right)$, it follows immediately that (2.4), (2.5), and (2.6)
lead to the response probability in (2.1).

##  3 When are the Linear Projection Parameters
Identical to the APEs?

In addition to being easy to interpret and readily extending to cases with
endogenous explanatory variables and unobserved heterogeneity, empirically the
OLS estimates of the LPM are often similar to the corresponding APEs from
nonlinear index models—particularly logit or probit. Wooldridge,, 2010,
Section 15.6 provides a discussion based on a results of Stoker (1986) that
helps one understand these empirical findings. Here we expand that discussion
to allow for an extension of the H-O framework.

It is useful to start with a general setting. Consider an index model

$P\left(y=1|\mathbf{x}\right)=G\left(\mathbf{x\beta}\right)=G\left(\beta_{1}+\beta_{2}x_{2}+\cdots+\beta_{K}x_{K}\right)\vspace{-3mm}$

where $G:\mathbb{R}\rightarrow\left[0,1\right]$. As argued in Wooldridge,,
2010, Section 15.6, the results of Stoker, (1986) imply that, if
$\left(x_{2},...,x_{K}\right)$ has a multivariate normal distribution and
$G\left(\cdot\right)$ is differentiable almost everywhere on $\mathbb{R}$
(with respect to Lebesgue measure), then

$\gamma_{j}=\beta_{j}E\left[g\left(\mathbf{x\beta}\right)\right]\equiv\alpha_{j}\text{,
}j=2,...,K,\vspace{-3mm}$

where $\gamma_{j}$ is the slope coefficients on $x_{j}$ in
$L\left(y|\mathbf{x}\right)=\mathbf{x\gamma}$, $g\left(\cdot\right)$ is the
almost everywhere derivative of $G\left(\cdot\right)$, and $\alpha_{j}$ is the
APE. The ramp function $R\left(\cdot\right)$ is differentiable everywhere
except at zero and one, and so it satisfies Stoker’s (1986) assumptions. The
result is that OLS consistently estimates the APEs, $\alpha_{j}$, even though
the $\alpha_{j}$ are attenuated versions of the $\beta_{j}$:

$\alpha_{j}=\beta_{j}\cdot P\left(0\leq\mathbf{x\beta}\leq 1\right)$

This equality holds even when $P\left(0\leq\mathbf{x\beta}\leq 1\right)$ can
be very close to zero. Horrace and Oaxaca, (2006), and many papers citing
their findings, focus on the inconsistency of OLS for $\beta_{j}$, failing to
recognize that the OLS estimators from the linear model could be consistent
for the more interesting quantities, the $\alpha_{j}$. This point is key to
our argument: If the model of the response probability is nonlinear so that
$0\leq p\left(\mathbf{x}\right)\leq 1$ is ensured, one should study estimation
of APEs, not underlying index parameters.

Clearly the assumption of multivariate normality of $\mathbf{x}$ is too
restrictive to be widely applicable. Nevertheless, the results of Stoker,
(1986) are suggestive, especially when combined with Ruud, (1983). Ruud
studies smooth nonlinear function forms that never hit the endpoints of the
unit interval, like probit and logit. In these cases, quasi-MLE identifies the
index coefficients up to scale. If $\mathbf{x}$ has a centrally symmetric
distribution—of which the multivariate normal is a special case—then Ruud’s
(1983) conditions hold. In the next section, we will find that when
$\mathbf{x}$ are symmetrically distributed with small variance (not too spread
out), the average partial effects are still approximated well by the LP
parameters; as the variance increases, however, the approximation breaks down.

To facilitate further discussion, including the simulations in the next
section, it is helpful to modify and extend the H-O setup. In particular,
write

$\displaystyle y^{\ast}$ $\displaystyle=\mathbf{x\beta}+u$ $\displaystyle
u|\mathbf{x}$ $\displaystyle\sim\mathrm{Uniform}\left(-a,a\right)$ (3.1)
$\displaystyle y$ $\displaystyle=1\left[y^{\ast}>0\right]\vspace{-3mm}$

for some $a>0$. Compared with H-O, we have shifted the intercept so that $u$
has a symmetric distribution about its mean of zero. Also, we allow $u$ to
have narrow or wide support, depending on $a$. While the latent error support
is not identified, it is a convenient device for generating data where the
unit interval for probabilities is binding to varying degrees. The CDF for the
$\mathrm{Uniform}\left(-\sqrt{3},\sqrt{3}\right)$ distribution, which has unit
variance, is graphed in Figure 2.

Figure 2: The CDF of $y$ with
$u|x\sim\mathrm{U}\left(-\sqrt{3},\sqrt{3}\right)$.

Given the latent variable model in (3.1), we can derive the response
probability:

$\displaystyle p\left(\mathbf{x}\right)$ $\displaystyle\equiv$ $\displaystyle
P\left(y=1|\mathbf{x}\right)=P\left(y^{\ast}\geq
0|\mathbf{x}\right)=P\left(u\geq-\mathbf{x\beta}|\mathbf{x}\right)=P\left(u\leq\mathbf{x\beta}|\mathbf{x}\right)=F_{u}\left(\mathbf{x\beta}\right)$
$\displaystyle=$ $\displaystyle 0\text{ if }\mathbf{x\beta}<-a$
$\displaystyle=$ $\displaystyle\frac{\mathbf{x\beta}+a}{2a}\text{ if
}-a\leq\mathbf{x\beta}\leq a$ $\displaystyle=$ $\displaystyle 1\text{ if
}\mathbf{x\beta}>a\vspace{-3mm}$

We write this function as $F_{u}\left(\mathbf{x\beta}\right)\equiv
R_{a}\left(\mathbf{x\beta}\right)$, which is a ramp function that is
nondifferentiable at $-a$ and $a$. For an $x_{j}$ with a positive coefficient,
the response probability has the same shape as in Figure 1. As $a$ increases
relative to $\mathbf{\beta}$, the response probability is linear over more of
the support of $\mathbf{x}$. If

$P\left(-a\leq\mathbf{x\beta}\leq a\right)=1\vspace{-3mm}$ (3.2)

then, with probability one,
$R_{a}\left(\mathbf{x\beta}\right)=\left(\mathbf{x\beta}+a\right)/2a$, a
linear function of $\mathbf{x}$. In this case, the partial effects are
constant and equal to $\beta_{j}/2a$, $j=2,...,K$. These are also the linear
projection parameters $\gamma_{j}$ and so OLS consistently estimates the APEs
under (3.2).

If $x_{j}$ is a continuous variable, we are interested in the APE defined as a
derivative, which exists with probability one when $\mathbf{x\beta}$ is
continuous. At $\mathbf{x\beta}\in\left\\{-a,a\right\\}$ the definition of the
partial effect is immaterial. To be concrete, take

$PE_{j}\left(\mathbf{x}\right)=\frac{\beta_{j}}{2a}\cdot
1\left[-a\leq\mathbf{x\beta}\leq a\right].\vspace{-3mm}$

Notice that $PE_{j}\left(\mathbf{x}\right)=0$ if $\mathbf{x\beta}<-a$ or
$\mathbf{x\beta}>a$ because we are on one of the flat parts of the ramp. This
feature of $PE_{j}\left(\mathbf{x}\right)$ is taken into account in computing
the APE:

$\alpha_{j}=E\left[PE_{j}\left(\mathbf{x}\right)\right]=\frac{\beta_{j}}{2a}\cdot
P\left(-a\leq\mathbf{x\beta}\leq a\right)\vspace{-3mm}$

Furthermore, the previous results based on Stoker (1986) still hold: If
$\left(x_{2},...,x_{K}\right)$ is multivariate normal then, again letting
$\gamma_{j}$ denote the LP parameter,

$\gamma_{j}=\alpha_{j}\vspace{-3mm}$ (3.3)

Again, it is important to understand that (3.3) holds even if
$P\left(-a\leq\mathbf{x\beta}\leq a\right)$ is close to zero (with zero being
ruled out); consequently, H-O’s discussion about the amount of inconsistency
in the OLS estimators if $P\left(-a\leq\mathbf{x\beta}\leq a\right)<1$ is
incomplete because they focus on $\beta_{j}$, not $\alpha_{j}$. In the
extended model (3.1), depending on the values of $a$ and
$P\left(-a\leq\mathbf{x\beta}\leq a\right)$, $\left|\alpha_{j}\right|$ need
not be smaller than $\left|\beta_{j}\right|$. The case that aligns with H-O is
$a=1/2$—so that the $Uniform\left(0,1\right)$ distributed has just been
shifted to have zero mean—in which case
$\left|\alpha_{j}\right|\leq\left|\beta_{j}\right|$, and the difference
between $\alpha_{j}$ and $\beta_{j}$ can be large. It is easily seen that
$\left|\alpha_{j}\right|<\left|\beta_{j}\right|$ for any $a\geq 1/2$.

We can also define APEs for discrete changes in the explanatory variables. For
example, if $x_{K}$ is binary, its partial effect is

$PE_{K}\left(\mathbf{x}\right)=R_{a}\left(\beta_{1}+\beta_{2}x_{2}+\cdots+\beta_{K-1}x_{K-1}+\beta_{K}\right)-R_{a}\left(\beta_{1}+\beta_{2}x_{2}+\cdots+\beta_{K-1}x_{K-1}\right),\vspace{-3mm}$

which corresponds to setting $x_{K}$ at its two values and obtaining the
difference in probabilities. Averaging across the joint distribution of the
other explanatory variables gives the APE:

$\alpha_{K}\equiv
APE_{K}=E\left[PE_{K}\left(\mathbf{x}\right)\right]\vspace{-3mm}$

If $x_{K}$ is a binary intervention or treatment indicator, $\alpha_{K}$ is
the average treatment effect.

Other than the case of multivariate normality of
$\left(x_{2},...,x_{K}\right)$, there is another case where the LP parameters,
$\gamma_{j}$, $j=2,...,K$, equal the APEs: $x_{2}$, …, $x_{K}$ are mutually
exclusive binary indicators that, along with a base group given by
$x_{2}=x_{3}=\cdots=x_{K}=0$, are exhaustive. See Angrist and Pischke,, 2009,
Section 3.1.4 and Wooldridge,, 2010, Section 15.2. If $x_{1}=1$ denotes the
base group then the APEs are simply

$\alpha_{j}=E\left(y|x_{j}=1\right)-E\left(y|x_{1}=1\right)\text{,
}j=2,...,K,\vspace{-3mm}$

and these are identical to the corresponding LPM coefficients.

Beyond the extreme cases described here, there appears to be no general theory
to determine when the LP coefficients will be the same or “close” to the APEs.
Many empirical applications include a combination of continuous, discrete, and
even mixed explanatory variables. Rarely do these have marginal symmetric
distributions, let alone a symmetric joint distribution. Plus, such
explanatory variables often appear appear as quadratics, interactions, and
other functional forms—which also do not have symmetric distributions. In
Section 5, we use simulations to shed light on when the the LPM coefficients
closely approximate the APEs—and when they do not. First, however, we describe
an estimator which is consistent under the ramp model.

##  4 The NLS Estimator of the Ramp Model

We have already seen how if
$P(\mathbf{x}_{i}\mathbf{\beta}\in\left[0,1\right])=1$, then OLS is consistent
for the $\beta_{j}$, which are equal to the APEs $\alpha_{j}$ in the case of a
continuous covariate $x_{j}$ under model (2.1). If the probability that
$\mathbf{x}_{i}\mathbf{\beta}$ lies outside the unit interval is nonzero, then
OLS is no longer consistent for the $\beta_{j}$, and it may or may not
approximate the $\alpha_{j}$ depending on the distribution of $\mathbf{x}$. In
addition to probit and logit quasi-MLE, it makes sense to consider an
estimator which is consistent if the ramp model is true. Of course, Bernoulli
MLE using the ramp model as the conditional response probability is not
feasible because the log-likelihood is not defined for
$\mathbf{x\beta}\notin(0,1)$. Instead, we consider nonlinear least squares
(NLS) using the piecewise ramp function $R(\mathbf{x}_{i}\mathbf{\beta})$ from
(2.1) as the conditional mean. In addition, since there may not be much
justification to think the ramp function is the true response probability, we
allow for general misspecification. Therefore, we define $\mathbf{\beta}_{o}$
as the pseudo-true value in the sense that $\mathbf{\beta}_{o}$ is the unique
solution to

$\displaystyle\underset{\mathbf{\beta}}{\text{min
}}E\left[\left(y_{i}-R(\mathbf{x}_{i}\mathbf{\beta})\right)^{2}\right]\equiv\underset{\mathbf{\beta}}{\text{min
}}Q(\mathbf{\beta}).\vspace{-3mm}$ (4.1)

We say that the model is misspecified if there is no such $\mathbf{\beta}$
such that $E[y|\mathbf{x}]=R(\mathbf{x}\mathbf{\beta})$. By construction,
$\mathbf{\beta}_{o}$ is the true coefficient when the model is correctly
specified and otherwise we view $R(\mathbf{x}\mathbf{\beta}_{o})$ as the best
mean squared error approximation to $E[y|x]$ over all ramp functions
$R(\mathbf{x}\mathbf{\beta})$.

As a sample analogue of (4.1), we define the objective function
$Q_{N}(\mathbf{\beta})$ as

$\displaystyle Q_{N}(\mathbf{\beta})\equiv$
$\displaystyle\frac{1}{N}\sum_{i=1}^{N}\left(y_{i}-R(\mathbf{x}_{i}\mathbf{\beta})\right)^{2}$
$\displaystyle=$
$\displaystyle\frac{1}{N}\sum_{i=1}^{N}\left(y_{i}^{2}1\left\\{\mathbf{x}_{i}\mathbf{\beta}\leq
0\right\\}+\left(y_{i}-\mathbf{x}_{i}\mathbf{\beta}\right)^{2}1\left\\{\mathbf{x}_{i}\mathbf{\beta}\in(0,1)\right\\}+(y_{i}-1)^{2}1\left\\{\mathbf{x}_{i}\mathbf{\beta}\geq
1\right\\}\right),\vspace{-3mm}$

where $N$ is the sample size. We define the NLS estimator
$\hat{\mathbf{\beta}}$ as

$\hat{\mathbf{\beta}}\equiv\text{arg}\underset{\mathbf{\beta}}{\text{min
}}\,Q_{N}(\mathbf{\beta}).\vspace{-3mm}$

The following theorem gives the consistency of the NLS estimator for the
pseudo-true value, allowing for misspecification of the conditional mean
model.

######  Theorem 1.

Let $\left\\{y_{i},\mathbf{x}_{i}\right\\}_{i=1}^{\infty}$ be an i.i.d.
sequence with $y$ only taking on values zero and one, and let
$R:\mathbb{R}\to[0,1]$ be the ramp function defined in (2.1). Suppose
$\mathbf{\beta}\in\mathbf{\mathcal{B}}$ such that
$\mathbf{\mathcal{B}}\subset\mathbb{R}^{K}$ is compact, and
$\mathbf{\beta}_{o}$ is identified in the sense that
$\forall\,\mathbf{\beta}\in\mathbf{\mathcal{B}},\mathbf{\beta}\neq\mathbf{\beta}_{o}$,

$E\left[\left(y_{i}-R(\mathbf{x}_{i}\mathbf{\beta}_{o})\right)^{2}\right]<E\left[\left(y_{i}-R(\mathbf{x}_{i}\mathbf{\beta})\right)^{2}\right]$

Then, $\hat{\mathbf{\beta}}\overset{p}{\rightarrow}\mathbf{\beta}_{o}$ as
$N\to\infty$.

The consistency result of Theorem 1 follows directly from Theorem 12.2 of
Wooldridge, (2010).

If $\mathbf{x}$ contains a continuously distributed $x_{j}$ and $\beta_{jo}$
is nonzero, then the probability of $\mathbf{x}_{i}\mathbf{\beta}_{o}$ being
equal to 0 or 1 is zero. Then, with suitable moment conditions on $\mathbf{x}$
(so the Leibniz integral rule applies), the FOC of the (4.1) is well defined
with probability 1 as follows:

$\displaystyle
E\left[\mathbf{x}_{i}^{\prime}u_{i}1\\{\mathbf{x}_{i}\mathbf{\beta}_{o}\in(0,1)\\}\right]=0,\vspace{-3mm}$
(4.2)

where $u_{i}(\mathbf{\beta})=y_{i}-R(\mathbf{x}_{i}\mathbf{\beta})$ and
$u_{i}\equiv u_{i}(\mathbf{\beta}_{0})$. Define the score function for random
draw $i$:

$\displaystyle\mathbf{s}_{i}(\mathbf{\beta})=-\mathbf{x}_{i}^{\prime}u_{i}(\mathbf{\beta})1\\{\mathbf{x}_{i}\mathbf{\beta}\in(0,1)\\}.\vspace{-3mm}$

Then, $\mathbf{\beta}_{o}$ solves $E[\mathbf{s}_{i}(\mathbf{\beta}_{o})]=0$.
The variance-covariance matrix of $\mathbf{s}_{i}(\mathbf{\beta})$ is

$\displaystyle\mathbf{\Omega}(\mathbf{\beta})=E\left[\mathbf{x}_{i}^{\prime}\mathbf{x}_{i}u_{i}(\mathbf{\beta})^{2}1\left\\{\mathbf{x}_{i}\mathbf{\beta}\in(0,1)\right\\}\right].\vspace{-3mm}$
(4.3)

The natural definition of the Jacobian of $\mathbf{s}_{i}(\mathbf{\beta})$ is

$\displaystyle\mathbf{A}_{i}(\mathbf{\beta})=\mathbf{x}_{i}^{\prime}\mathbf{x}_{i}1\left\\{\mathbf{x}_{i}\mathbf{\beta}\in(0,1)\right\\}.\vspace{-3mm}$

For the similar reason as (4.2), the Hessian of $Q(\mathbf{\beta})$ is well-
defined with probability 1 at $\mathbf{\beta}_{o}$ as follows

$\displaystyle\mathbf{A}(\mathbf{\beta}_{o})=E\left[\mathbf{x}_{i}^{\prime}\mathbf{x}_{i}1\left\\{\mathbf{x}_{i}\mathbf{\beta}_{o}\in(0,1)\right\\}\right].\vspace{-3mm}$
(4.4)

Note that (4.3) and (4.4) are the same whether the conditional mean model is
correctly specified or not. Therefore, the following asymptotic distribution
result allows for misspecification of the model.

######  Theorem 2.

Suppose that the assumptions from Theorem 1 hold, and (i) $\mathbf{\beta}_{o}$
is an interior point of $\mathbf{\mathcal{B}}$; (ii) $\mathbf{x}_{i}$ contains
a continuously distributed random variable with a nonzero coefficient; (iii)
$E\left\|\mathbf{x}_{i}\right\|^{2}<\infty$ and
$E\left[\mathbf{x}_{i}^{\prime}\mathbf{x}_{i}1\left\\{\mathbf{x}_{i}\mathbf{\beta}_{o}\in(0,1)\right\\}\right]>0$,
where $\|.\|$ denotes the $l^{2}-norm$. Then, as $N\to\infty$,

$\displaystyle\sqrt{N}\left(\hat{\mathbf{\beta}}-\mathbf{\beta}_{o}\right)\Rightarrow\mathbb{N}(0,\mathbf{A}(\mathbf{\beta}_{o})^{-1}\mathbf{\Omega}(\mathbf{\beta}_{o})\mathbf{A}(\mathbf{\beta}_{o})^{-1}).\vspace{-3mm}$

The proof of Theorem 2 is given in the Appendix. The asymptotic normality
results does not follow directly from the M-estimator due to the non-
smoothness of the objective function. We therefore leverage an asymptotic
normality result for estimators with non-smooth objective function from Newey
and McFadden, (1994).

To estimate $\mathbf{\beta}_{o}$, H-O suggest running OLS on a trimmed sample
(i.e., those observations for which initial OLS fitted values are inside the
unit interval) to reduce bias. We find in practice that a single round of
trimming may not reduce the bias for the APEs in the cases where OLS is not
consistent for them. However, we find an iterative trimming OLS procedure
(ITO) does reduce the bias for estimating APEs, as well as
$\mathbf{\beta}_{o}$.444The procedure goes: 1) estimate the LPM by OLS. 2)
Compute fitted values. 3) Drop observations with fitted values outside the
unit interval, and 4) Repeat starting at 1) until no further observations are
dropped. In fact, we find in simulations that the NLS estimates are
numerically the same as the ITO estimates up to machine precision.555With some
DGP, it was occasionally necessary to specify OLS starting values for the NLS
function evaluator for this result to hold. It turns out that ITO is
implicitly minimizing the NLS sample objective function using the OLS
estimates as starting values and following the Newton-Raphson numerical
method, which is iterative (see Wooldridge,, 2010, Section 12.7.1). Given an
estimate $\mathbf{\beta}^{\left\\{g\right\\}}$, the next iteration is given
(using our notation) by

$\displaystyle\mathbf{\beta}^{\left\\{g+1\right\\}}$
$\displaystyle=\mathbf{\beta}^{\left\\{g\right\\}}-\left[N^{-1}\sum_{i=1}^{N}\mathbf{A}_{i}(\mathbf{\beta}^{\left\\{g\right\\}})\right]^{-1}N^{-1}\sum_{i=1}^{N}\mathbf{s}_{i}(\mathbf{\beta}^{\left\\{g\right\\}})$
$\displaystyle=\mathbf{\beta}^{\left\\{g\right\\}}+\left[N^{-1}\sum_{i=1}^{N}\mathbf{x}_{i}^{\prime}\mathbf{x}_{i}1\left\\{\mathbf{x}_{i}\mathbf{\beta}^{\left\\{g\right\\}}\in(0,1)\right\\}\right]^{-1}N^{-1}\sum_{i=1}^{N}\mathbf{x}_{i}^{\prime}\left(y_{i}-\mathbf{x}_{i}\mathbf{\beta}^{\left\\{g\right\\}}\right)1\left\\{\mathbf{x}_{i}\mathbf{\beta}^{\left\\{g\right\\}}\in(0,1)\right\\}$
$\displaystyle=\left[N^{-1}\sum_{i=1}^{N}\mathbf{x}_{i}^{\prime}\mathbf{x}_{i}1\left\\{\mathbf{x}_{i}\mathbf{\beta}^{\left\\{g\right\\}}\in(0,1)\right\\}\right]^{-1}N^{-1}\sum_{i=1}^{N}\mathbf{x}_{i}^{\prime}y_{i}1\left\\{\mathbf{x}_{i}\mathbf{\beta}^{\left\\{g\right\\}}\in(0,1)\right\\}.\vspace{-3mm}$

The second equality above substitutes our expressions for
$\mathbf{s}_{i}(\mathbf{\beta})$ and $\mathbf{A}_{i}(\mathbf{\beta})$ and uses
the fact that $R(\mathbf{x}_{i}\mathbf{\beta})=\mathbf{x}_{i}\mathbf{\beta}$
for $\mathbf{x}_{i}\mathbf{\beta}\in(0,1)$. This shows that
$\mathbf{\beta}^{\left\\{g+1\right\\}}$ is simply the OLS estimator on the
sample with $\mathbf{x}_{i}\mathbf{\beta}^{\left\\{g\right\\}}\in(0,1)$.

As a consequence, the preceding consistency and asymptotic normality results
for the NLS estimator justify using the ITO procedure to reduce the OLS bias.
However, it is worth mentioning that, at least in Stata, the pre-loaded NLS
solver (the“nl” command) may have a performance advantage over ITO in
practice. We find in simulations that ITO can result in a dead loop when only
a very small portion of observations are left for estimation after iterative
trimming. The pre-loaded NLS algorithm continues to work well in those cases.

Taking the sample analogue of the asymptotic variance from Theorem 2, we
define a variance estimator of
$\sqrt{N}(\hat{\mathbf{\beta}}-\mathbf{\beta}_{o})$ as

$\displaystyle\mathbf{\hat{V}}=\mathbf{A}_{N}(\hat{\mathbf{\beta}})^{-1}\mathbf{\Omega}_{N}(\hat{\mathbf{\beta}})\mathbf{A}_{N}(\hat{\mathbf{\beta}})^{-1},\vspace{-3mm}$

where
$\mathbf{A}_{N}(\hat{\mathbf{\beta}})=N^{-1}\sum_{i=1}^{N}\mathbf{x}_{i}^{\prime}\mathbf{x}_{i}1\\{\mathbf{x}_{i}\hat{\mathbf{\beta}}\in(0,1)\\}$,
$\mathbf{\Omega}_{N}(\hat{\mathbf{\beta}})=N^{-1}\sum_{i=1}^{N}\mathbf{x}_{i}^{\prime}\mathbf{x}_{i}\hat{u}^{2}_{i}1\\{\mathbf{x}_{i}\mathbf{\hat{\beta}}\in(0,1)\\}$,
and $\hat{u}_{i}=y_{i}-R(\mathbf{x}_{i}\hat{\mathbf{\beta}})$. Standard errors
are obtained the usual way from $\mathbf{\hat{V}}/N$. The next theorem gives
the consistency result of the variance estimator.

######  Theorem 3.

Under the same assumption of Theorem 2 and $E\|x\|^{4}<\infty$, as
$N\to\infty$, $\mathbf{\hat{V}}\stackrel{{\scriptstyle
p}}{{\to}}\mathbf{A}(\mathbf{\beta}_{o})^{-1}\mathbf{\Omega}(\mathbf{\beta}_{o})\mathbf{A}(\mathbf{\beta}_{o})^{-1}$.

The proof of Theorem 3 is given in the Appendix. As before, we are interested
in the APE. Consider the best ramp approximation in (4.1), the APE of a
continuous random variable $x_{k}$ is defined as

$\displaystyle APE_{k}=E\left[\frac{\partial
R(\mathbf{x}_{i}\mathbf{\beta}_{o})}{\partial
x_{k}}\right]=\beta_{ko}P\left(\mathbf{x}_{i}\mathbf{\beta}_{o}\in(0,1)\right).\vspace{-3mm}$

A sample-analogue estimator of the APE is then given by

$\displaystyle
A\widehat{P}E_{k}=\hat{\mathbb{\beta}}_{k}\frac{1}{N}\sum_{i=1}^{N}1\left\\{\mathbf{x}_{i}\hat{\mathbf{\beta}}\in(0,1)\right\\}\vspace{-3mm}$

Define
$g(\mathbf{x}_{i},\mathbf{\beta})=\beta_{k}1\\{\mathbf{x}_{i},\mathbf{\beta}_{o}\\}$,
$\delta_{o}=E[g(\mathbf{x}_{i},\mathbf{\beta}_{o})]$, and
$\mathbf{G}_{o}=\nabla_{\mathbf{\beta}}g(\mathbf{x}_{i},\mathbf{\beta}_{o})$.
Following problem 12.17 of Wooldridge (2010), the asymptotic variance of the
estimated APE is given by

$\displaystyle
AVar\left(\sqrt{N}\left(A\widehat{P}E_{k}-APE_{k}\right)\right)=Var\left(g(\mathbf{x}_{i},\mathbf{\beta}_{o})-\delta_{o}-\mathbf{G}_{o}\mathbf{A}(\mathbf{\beta}_{o})^{-1}\mathbf{s}_{i}(\mathbf{\beta}_{o})\right),\vspace{-3mm}$

where $\mathbf{G}_{o}$ is a $1\times K$ vector with the $k^{th}$ element being
$p_{o}\equiv P\left(\mathbf{x}_{i}\mathbf{\beta}_{o}\in(0,1)\right)$ and all
else $0$. The asymptotic variance can be estimated by the sample variance of
$g(\mathbf{x}_{i},\hat{\mathbf{\beta}})-\hat{\delta}-\widehat{\mathbf{G}}\mathbf{A}_{N}(\hat{\mathbf{\beta}})^{-1}\mathbf{s}_{i}(\hat{\mathbf{\beta}})$,
where
$\hat{\delta}=\frac{1}{N}\sum_{i=1}^{N}g(\mathbf{x}_{i},\hat{\mathbf{\beta}})$,
$\widehat{\mathbf{G}}$ is a $1\times K$ vector with the $k^{th}$ element being
$\hat{p}=\frac{1}{N}\sum_{i=1}^{N}1\left\\{\mathbf{x}_{i}\hat{\mathbf{\beta}}\in(0,1)\right\\}$.

The APE for a discrete random variable $x_{k}$ can be defined as

$\displaystyle
APE_{k}=E\left[R(\mathbf{x}_{i,-k}\mathbf{\beta}_{-ko}+\beta_{ko})-R(\mathbf{x}_{i,-k}\mathbf{\beta}_{-ko})\right].\vspace{-3mm}$

A sample analogue estimator of $APE_{k}$ is given by

$\displaystyle
A\widehat{P}E_{k}=\frac{1}{N}\sum_{i=1}^{N}R(\mathbf{x}_{i,-k}\hat{\mathbf{\beta}}_{-k}+\hat{\beta}_{k})-R(\mathbf{x}_{i,-k}\hat{\mathbf{\beta}}_{-k}).\vspace{-3mm}$

The asymptotic variance can be found and estimated in a similar manner as the
continuous case.

##  5 Simulations

In this section we present several Monte Carlo simulations that provide
insights into the behavior of different modeling/estimation approaches. The
LPM is estimated by OLS and the ramp function is estimated by NLS. For NLS,
the average partial effects are estimated based on averages of derivatives and
differences of the ramp function. These resemble the familiar formulas for the
linear model, though the individual unit partial effects need to be scaled by
$1\left[0\leq\widehat{y}\leq 1\right]$ before averaging, where $\widehat{y}$
corresponds to the fitted values for each estimator. The logit and probit
parameters are estimated by the (quasi-) maximum likelihood estimator, and
then the average partial effects are estimated using the usual APE formulas.
We further consider a nonparametric model estimated by the local linear
estimator and the APEs are estimated by the sample average of the partial
effects. We used Stata®17 for simulation. The Stata code is available upon
request.

Initially, the true models take the form (we are dropping $o$ on beta here)

$y=1\left[\beta_{0}+\beta_{1}x_{1}+\beta_{2}x_{2}+u>0\right],\vspace{-3mm}$
(5.1)

where $u$ is independent of $\left(x_{1},x_{2}\right)$ with $u\sim
Uniform\left(-a,a\right)$ for $a>0$. The choice of $a$ is important because it
governs how close to linear is the response probability for a given
$\mathbf{\beta}$. The variable $x_{1}$ is continuous and $x_{2}$ is binary;
they are generated to be correlated. We indicate the intercept in the index
function by $\beta_{0}$.

When $u\sim Uniform(-a,a)$, the ramp model is correctly specified, but the LPM
is misspecified to varying degrees. For small $a$, the kinks in the ramp
function are binding and the LPM can provide a poor approximation to the
response probability. Naturally, the logit and probit models are always
misspecified in this case. As stated before, here we focus on the APEs rather
than the underlying parameters or how well the models approximate the true
response probability.

The sample size is $N=1,000$ and $1,000$ replications are used. The population
(or true) APEs are not available in closed form, and so we simulate these
along with the estimators. In the tables to follow, the columns labeled
“Simulated Truth” include the empirical means and standard deviations of the
sample average partial effects at the true parameter values. We also simulate
the probabilities

$P\left(-a\leq\mathbf{x\beta}\leq a\right)\vspace{-3mm}$

and

$P\left(0\leq\widehat{y}\leq 1\right),\vspace{-3mm}$

where $\hat{y}$ refers to predicted values for the linear index. The first of
these tells us how binding are the ramp function inflection points. The second
is practically relevant because researchers often check the fraction of fitted
values outside the unit interval as a way to determine the adequacy of the
LPM. The simulations show that having a large fraction of fitted values in
$\left[0,1\right]$ is neither necessary nor sufficient for OLS to produce
accurate estimates of the APEs.666H-O do note that consistency of OLS can no
longer be shown as soon as one observation has a true index outside the unit
interval.

We also consider the case where an interaction term, $x_{1}\cdot x_{2}$, is
included in the model, and the researcher includes in interaction in the
specification. In the Appendix, we choose
$u\sim\mathrm{Normal}\left(0,1\right)$ to compare the LPM with logit and
probit when the probit model is correct.

###  5.1 Symmetrically Distributed Explanatory Variables

In the first design, $\left(x_{1},x_{2}\right)$ are generated as

$\displaystyle x_{1}$ $\displaystyle=v/\sqrt{2}+e/\sqrt{2}$ $\displaystyle
x_{2}$ $\displaystyle=1\left[v/2+r>0\right],\vspace{-3mm}$

where $v$, $e$, and $r$ are independent standard normals. The index parameters
are set as

$\left(\beta_{0},\beta_{1},\beta_{2}\right)=\left(0.1,0.2,-0.3\right)\vspace{-3mm}$

The binary outcome $y$ is generated as in (5.1) with
$u\sim\mathrm{Uniform}\left(-a,a\right)$, where
$a\in\left\\{1/4,1/2,1\right\\}$. The case $a=1/2$ essentially corresponds to
H-O. When $a=1$, the response probability is essentially linear for these
index parameter values. When $a=1/4$, the kinks are binding and
$\mathbf{x\beta}$ is often outside the interval $\left[-a,a\right]$.

Table 1 reports the findings when $a=1/2$. There is a small probability that
$\mathbf{x\beta}\notin\left[-a,a\right]$ – roughly, about 0.013. Moreover,
across all simulations, about $1.3\%$ of the OLS fitted values are outside the
unit interval. The pattern is clear: All the estimators of the APEs show very
little bias and have the same precision. This is true for the continuous
variable, $x_{1}$, and the binary variable, $x_{2}$. Note that this is not
predicted by application of the Stoker results because $x_{2}$ is a discrete
variable. Nevertheless, this table illustrates what is often observed in
practice: the LPM coefficients estimated by OLS are often close to the probit
and logit APEs.

Table 1. No interaction, $x_{1}$ normal, $x_{2}$ sym. binary, $u\sim
U(-0.5,0.5)$
---
$N=1000$ | Simulated | LPM | Ramp | Probit | Logit | Nonpar.
|  | Truth* | (OLS) | (NLS) | (QMLE) | (QMLE) | (LL)
$APE_{1}$ | mean | 0.1975 | 0.1971 | 0.1972 | 0.1981 | 0.1969 | 0.1980
sd | 0.0007 | 0.0132 | 0.0134 | 0.0131 | 0.0132 | 0.0141
$APE_{2}$ | mean | -0.2925 | -0.2960 | -0.2919 | -0.2877 | -0.2851 | -0.2950
sd | 0.0010 | 0.0299 | 0.0291 | 0.0283 | 0.0283 | 0.0291
$P(y=1)$ | 0.4509 |  |  |  |  |
$P(0\leq\widehat{y}\leq 1)$ |  | 0.9875 | 0.9860 | 1.0000 | 1.0000 | 0.9867
$P(-a\leq\mathbf{x\beta}\leq a)$ | 0.9873 |  |  |  |  |

  * •

*This column contains the empirical means and standard deviations of the sample average partial effects at the true parameter values.

The story does not change when the flat parts of the ramp function are
strongly binding. In Table 2, ${\small P}\left(-a\leq\mathbf{x\beta}\leq
a\right)$ is only about $0.78$, and about 11.8 percent of the OLS fitted
values are outside $\left[0,1\right]$. And yet, for estimating the APEs, the
LPM does essentially as well as probit and logit, with logit having perhaps a
bit less bias. But, given the simulation error, these are not to be dwelt
upon. Table 3 shows the case where ${\small P}\left(-a\leq\mathbf{x\beta}\leq
a\right)$ is exactly one. We would expect the LPM to work very well in this
case, and it does—especially for the continuous variable $x_{1}$. What is,
perhaps, more surprising is that probit and logit work just as well, even
though the true response probability is linear over the support of
$\mathbf{x\beta}$. These findings are a good reminder of why statements such
as “the linear probability model is preferred to probit because the latter
assumes normality” are not just misleading: they are wrong. In the end, what
we care about is how well each approach approximates the partial effects on
$P\left(y=1|\mathbf{x}\right)$. When we consider the APEs, all methods do well
even when the response probability has the peculiar ramp shape.

Table 2. No interaction, $x_{1}$ normal, $x_{2}$ sym. binary, $u\sim
U(-0.25,0.25)$
---
$N=1000$ | Simulated | LPM | Ramp | Probit | Logit | Nonpar.
|  | Truth* | (OLS) | (NLS) | (QMLE) | (QMLE) | (LL)
$APE_{1}$ | mean | 0.3124 | 0.3137 | 0.3129 | 0.3152 | 0.3127 | 0.3270
sd | 0.0051 | 0.0097 | 0.0098 | 0.0092 | 0.0096 | 0.0112
$APE_{2}$ | mean | -0.4414 | -0.4774 | -0.4418 | -0.4424 | -0.4385 | -0.4567
sd | 0.0058 | 0.0229 | 0.0206 | 0.0198 | 0.0202 | 0.0240
$P(y=1)$ | 0.4215 |  |  |  |  |
$P(0\leq\widehat{y}\leq 1)$ |  | 0.8822 | 0.7768 | 1.0000 | 1.0000 | 0.8770
$P(-a\leq\mathbf{x\beta}\leq a)$ | 0.7809 |  |  |  |  |

  * •

*This column contains the empirical means and standard deviations of the sample average partial effects at the true parameter values.

Table 3. No interaction, $x_{1}$ normal, $x_{2}$ sym. binary, $u\sim U(-1,1)$
---
$N=1000$ | Simulated | LPM | Ramp | Probit | Logit | Nonpar.
|  | Truth* | (OLS) | (NLS) | (QMLE) | (QMLE) | (LL)
$APE_{1}$ | mean | 0.1000 | 0.1001 | 0.1001 | 0.1001 | 0.1001 | 0.1002
sd | 0.0000 | 0.0159 | 0.0159 | 0.0159 | 0.0159 | 0.0163
$APE_{2}$ | mean | -0.1500 | -0.1509 | -0.1509 | -0.1496 | -0.1492 | -0.1507
sd | 0.0000 | 0.0331 | 0.0331 | 0.0326 | 0.0325 | 0.0335
$P(y=1)$ | 0.4755 |  |  |  |  |
$P(0\leq\widehat{y}\leq 1)$ |  | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 0.9980
$P(-a\leq\mathbf{x\beta}\leq a)$ | 1.0000 |  |  |  |  |

  * •

*This column contains the empirical means and standard deviations of the sample average partial effects at the true parameter values.

We also generated the outcome $y$ using an interaction between $x_{1}$ and
$x_{2}$, with $u$ still having a uniform distribution. Specifically,

$y=1\left[\beta_{0}+\beta_{1}x_{1}+\beta_{2}x_{2}+\beta_{3}\left(x_{1}\cdot
x_{2}\right)+u>0\right]\vspace{-3mm}$

Remember, both $x_{1}$ and $x_{2}$ have symmetric distributions, but this
functional form falls outside Stoker’s results because $x_{2}$ is discrete and
so is $x_{1}\cdot x_{2}$: it has a mass point at zero and is otherwise
continuous. Across a few parameter settings, the three approaches—where the
interaction term is included in the estimation—delivered similar estimated
APEs that were close to the sample “true” APEs. (As previously, probit, logit,
and OLS approaches use a misspecified response probability.) The parameters in
that case are set at

$\left(\beta_{0},\beta_{1},\beta_{2},\beta_{3}\right)=\left(0.1,0.2,-0.3,-0.3\right)\vspace{-3mm}$

Table 4. With interaction, $x_{1}$ normal, $x_{2}$ sym. binary, $u\sim
U(-0.5,0.5)$
---
$N=1000$ | Simulated | LPM | Ramp | Probit | Logit | Nonpar.
|  | Truth* | (OLS) | (NLS) | (QMLE) | (QMLE) | (LL)
$APE_{1}$ | mean | 0.0495 | 0.0493 | 0.0494 | 0.0492 | 0.0491 | 0.0495
sd | 0.0049 | 0.0143 | 0.0146 | 0.0144 | 0.0144 | 0.0156
$APE_{2}$ | mean | -0.2989 | -0.3002 | -0.2990 | -0.2941 | -0.2929 | -0.2997
sd | 0.0099 | 0.0329 | 0.0326 | 0.0321 | 0.0322 | 0.0330
$P(y=1)$ | 0.3898 |  |  |  |  |
$P(0\leq\widehat{y}\leq 1)$ |  | 0.9938 | 0.9930 | 1.0000 | 1.0000 | 0.9946
$P(-a\leq\mathbf{x\beta}\leq a)$ | 0.9954 |  |  |  |  |

  * •

*This column contains the empirical means and standard deviations of the sample average partial effects at the true parameter values.

Tables 4 and 5 show the simulation findings. In Table 4, when the support of
$u$ is moderately wide, all methods perform about equally well. They have
little bias and their precisions are practically identical. In Table 5, when
the support of $u$ narrows to $\left(-1/4,1/4\right)$, the OLS estimation of
LPM still does as well as the other estimations. These findings would seem to
go against conventional wisdom because there is a non-trivial fraction of
fitted values outside the unit interval, about $0.15$. Finally, we note that
the OLS estimator for $P(0\leq\widehat{y}\leq 1)$ can be severley biased for
$P(-a\leq\mathbf{x\beta}\leq a)$.

Table 5. With interaction, $x_{1}$ normal, $x_{2}$ sym. binary, $u\sim
U(-0.25,0.25)$
---
$N=1000$ | Simulated | LPM | Ramp | Probit | Logit | Nonpar.
|  | Truth* | (OLS) | (NLS) | (QMLE) | (QMLE) | (LL)
$APE_{1}$ | mean | 0.1104 | 0.1115 | 0.1108 | 0.1105 | 0.1102 | 0.1187
sd | 0.0079 | 0.0125 | 0.0127 | 0.0126 | 0.0128 | 0.0138
$APE_{2}$ | mean | -0.5154 | -0.5401 | -0.5150 | -0.5074 | -0.5050 | -0.5256
sd | 0.0150 | 0.0275 | 0.0272 | 0.0268 | 0.0271 | 0.0280
$P(y=1)$ | 0.3094 |  |  |  |  |
$P(0\leq\widehat{y}\leq 1)$ |  | 0.8503 | 0.6805 | 1.0000 | 1.0000 | 0.9640
$P(-a\leq\mathbf{x\beta}\leq a)$ | 0.6850 |  |  |  |  |

  * •

*This column contains the empirical means and standard deviations of the sample average partial effects at the true parameter values.

###  5.2 Asymmetrically Distributed Explanatory Variables

The story changes markedly when the distributions of $x_{1}$ and $x_{2}$ are
asymmetric. With $v$, $e$, and $r$ generated as before, $x_{1}$ and $x_{2}$
are now generated as

$\displaystyle x_{1}$ $\displaystyle=$
$\displaystyle\exp\left(0.5+v/2+e/2\right)$ $\displaystyle x_{2}$
$\displaystyle=$ $\displaystyle 1\left[-0.5+v/2+r>0\right],\vspace{-3mm}$

so that $x_{1}$ has a lognormal distribution. The variable $x_{2}$ is still
binary but the response probability is well below $0.5$. Tables 6 and 7 repeat
the same experiments as Table 2 and Table 3, with $u\sim U(-0.25,0.25)$ and
$u\sim U(-1,1)$ respectively, but with the covariates generated as above. The
results for $u\sim U(-0.5,0.5)$ possess a similar pattern and thus omitted for
brevity. The parameter values are, again,

$\left(\beta_{0},\beta_{1},\beta_{2}\right)=\left(0.1,0.2,-0.3\right)\vspace{-3mm}$

Table 6. No interaction, $x_{1}$ asym., $x_{2}$ asym. binary, $u\sim
U(-0.25,0.25)$
---
$N=1000$ | Simulated | LPM | Ramp | Probit | Logit | Nonpar.
|  | Truth* | (OLS) | (NLS) | (QMLE) | (QMLE) | (LL)
$APE_{1}$ | mean | 0.2203 | 0.0478 | 0.2208 | 0.2131 | 0.2100 | 0.1573
sd | 0.0063 | 0.0100 | 0.0163 | 0.0149 | 0.0154 | 0.0289
$APE_{2}$ | mean | -0.3251 | -0.2772 | -0.3252 | -0.3283 | -0.3271 | -0.3420
sd | 0.0083 | 0.0263 | 0.0213 | 0.0203 | 0.0205 | 0.0233
$P(y=1)$ | 0.8149 |  |  |  |  |
$P(0\leq\widehat{y}\leq 1)$ |  | 0.9295 | 0.5483 | 1.0000 | 1.0000 | 0.8720
$P(-a\leq\mathbf{x\beta}\leq a)$ | 0.5507 |  |  |  |  |

  * •

*This column contains the empirical means and standard deviations of the sample average partial effects at the true parameter values.

Table 7. No interaction, $x_{1}$ asym., $x_{2}$ asym. binary, $u\sim U(-1,1)$
---
$N=1000$ | Simulated | LPM | Ramp | Probit | Logit | Nonpar.
|  | Truth* | (OLS) | (NLS) | (QMLE) | (QMLE) | (LL)
$APE_{1}$ | mean | 0.0932 | 0.0477 | 0.0939 | 0.1083 | 0.1102 | 0.0889
sd | 0.0008 | 0.0092 | 0.0095 | 0.0095 | 0.0100 | 0.0137
$APE_{2}$ | mean | -0.1379 | -0.1040 | -0.1391 | -0.1428 | -0.1426 | -0.1412
sd | 0.0012 | 0.0318 | 0.0289 | 0.0279 | 0.0280 | 0.0296
$P(y=1)$ | 0.6453 |  |  |  |  |
$P(0\leq\widehat{y}\leq 1)$ |  | 0.9802 | 0.9317 | 1.0000 | 1.0000 | 0.9803
$P(-a\leq\mathbf{x\beta}\leq a)$ | 0.9322 |  |  |  |  |

  * •

*This column contains the empirical means and standard deviations of the sample average partial effects at the true parameter values.

The findings in Table 7 are striking. Even though ${\small
P}\left(-a\leq\mathbf{x\beta}\leq a\right)$ is high—around $0.93$—and the OLS
fitted values are very rarely outside the unit interval (only about 2 percent
of the time), the OLS estimators of the LPM are badly biased for the APEs and
are notably worse than other methods. And among the other estimators, Ramp/NLS
has a smaller bias in terms of both APEs. In Table 6, it is the same case that
Ramp/NLS continues to work relatively better than any other methods in terms
of APEs, even though the fraction of fitted values of Ramp/NLS within the unit
interval is only about 0.55. The results with an interaction term is similar
and so are skipped for brevity.

###  5.3 Additional Simulations

Comparing Table 6 and Table 3, it seems like the symmetric distribution of
$\mathbf{x}$ is the key condition for OLS to consistently estimate $APEs$. In
Table 8, however, we give a counterexample where $\mathbf{x}$ is symmetrically
distributed, but $x_{1}$ has a $Uniform(-10,10)$ distribution. Compared to the
normal distribution of $x_{1}$ in Section 5.1, this distribution has higher
variance and lacks a mode. Unlike before, OLS poorly approximates the APEs in
this case. In addition to symmetry, therefore, the modality and higher moments
of the covariates may also be important in determining the performance of OLS.

Table 8. No interaction, $x_{1}\sim U(-10,10)$, $x_{2}$ sym. binary, $u\sim
U(-1,1)$
---
$N=1000$ | Simulated | LPM | Ramp | Probit | Logit | Nonpar.
|  | Truth* | (OLS) | (NLS) | (QMLE) | (QMLE) | (LL)
$APE_{1}$ | mean | 0.0501 | 0.0682 | 0.0501 | 0.0500 | 0.0498 | 0.0512
sd | 0.0016 | 0.0009 | 0.0016 | 0.0016 | 0.0017 | 0.0019
$APE_{2}$ | mean | -0.0750 | -0.0750 | -0.0754 | -0.0754 | -0.0754 | -0.0752
sd | 0.0022 | 0.0191 | 0.0184 | 0.0172 | 0.0180 | 0.0180
$P(y=1)$ | 0.4801 |  |  |  |  |
$P(0\leq\widehat{y}\leq 1)$ |  | 0.7325 | 0.4960 | 1.0000 | 1.0000 | 0.9440
$P(-a\leq\mathbf{x\beta}\leq a)$ | 0.5006 |  |  |  |  |

  * •

*This column contains the empirical means and standard deviations of the sample average partial effects at the true parameter values.

An additional set of simulations are included in the Appendix which suggest
our findings have more to do with the joint distribution of the explanatory
variables than with the choice of distribution for the latent model error.
Tables 11 corresponds to the DGPs of Tables 1-3, and Table 12 corresponds to
the DGP of Tables 4-5, but the Appendix tables are based on standard normal
error terms, corresponding to the probit model. As in main text, each of the
estimators (including OLS) has small bias for the APEs in Tables 11-12. Table
13 in the Appendix corresponds to Tables 6-7 of the main text, but also with a
normally distributed error. In this case, we find (as in the main text), OLS
has larger bias for the APEs than do the nonlinear or nonparametric
estimators.

##  6 Mortgage Approval Probabilities and Race

As an illustration of linear and nonlinear estimators for binary response
models, we revisit the analysis of discrimination in mortgage lending
decisions from Hunter and Walker, (1996).777We use a version of the loan
applications dataset provided by Mary Beth Walker for Wooldridge (2019). The
cultural affinity hypothesis posits that white loan officers may “rely more
heavily on basic objective loan application information in appraising the
creditworthiness of minorities” due to a lack of cultural familiarity. We
compare linear and nonlinear estimates of the average effect of being white on
the probability of loan approval, holding constant a number of loan, property,
and borrower characteristics. Table 9 presents basic summary statistics for
the dependent variable “approve” and 23 covariates.

Table 9: Loan Approval Summary Statistics ($N=1989$)
---
Variable | Description | Mean | SD | Skew. | Kurt.
approve | =1 if loan approved | 0.88 | 0.33 | -2.30 | 6.29
white | =1 if white | 0.85 | 0.36 | -1.91 | 4.64
loanamt | Loan amount $1000s | 143.25 | 80.52 | 3.13 | 20.36
suffolk | =1 if in Suffolk County | 0.15 | 0.36 | 1.91 | 4.66
appinc | Applicant income $1000s | 84.68 | 87.06 | 5.26 | 36.70
unit | Number of units in property | 1.12 | 0.44 | 4.01 | 19.89
married | =1 if applicant married | 0.66 | 0.47 | -0.67 | 1.45
dep | Number of dependents | 0.77 | 1.10 | 1.47 | 5.33
emp | Years employed in line of work | 0.21 | 1.00 | 6.69 | 50.57
yjob | Years at this job | 0.45 | 1.12 | 5.32 | 36.18
atotinc | Total monthly income | 5195.55 | 5269.06 | 6.36 | 65.34
self | =1 if self employed | 0.13 | 0.34 | 2.21 | 5.89
other | Other financing $1000s | 2.37 | 28.23 | 26.80 | 886.84
rep | Number of credit reports | 1.50 | 0.99 | 1.45 | 7.37
pubrec | =1 if filed bankruptcy | 0.07 | 0.25 | 3.40 | 12.59
hrat | Housing expense % of total inc. | 24.79 | 7.12 | 0.25 | 6.74
obrat | Other obligations % of total inc. | 32.39 | 8.26 | 0.44 | 7.40
cosign | =1 if there is a cosigner | 0.03 | 0.17 | 5.65 | 32.92
sch |  =1 if $>$ 12 years schooling | 0.77 | 0.42 | -1.29 | 2.68
mortno | =1 if no mortgage history | 0.33 | 0.47 | 0.71 | 1.51
mortlat1 | =1 if one or two late payments | 0.02 | 0.14 | 7.03 | 50.36
mortlat2 | =1 if more than two late payments | 0.01 | 0.10 | 9.58 | 92.72
chist |  =0 if accounts are delinq. $\geq$60 days | 0.84 | 0.37 | -1.83 | 4.35
loanprc | Loan amount / purchase price | 0.77 | 0.19 | 0.44 | 14.39

For our index model, we include interactions between “white” and all other
explanatory variables to allow for the factors like loan amount and credit
history to have a differential impact on approval probability by race. Let $w$
denote “white” and $\mathbf{z}$ be a vector including the 22 other covariates,
so that $\mathbf{x}=\left\\{1,\mathbf{z},w,w\mathbf{z}\right\\}$ and
$\mathbf{\beta}=\left\\{\mathbf{\beta}_{0},\mathbf{\beta}_{z},\beta_{w},\mathbf{\beta}_{wz}\right\\}$,
where $\mathbf{\beta}_{0}$ is the intercept, $\mathbf{\beta}_{z}$ and
$\beta_{w}$ are the coefficients on $\mathbf{z}$ and $w$, respectively, while
$\mathbf{\beta}_{wz}$ is the coefficient on $w\mathbf{z}$. Then the partial
effects we average are formed by evaluating the difference in the
probabilities evaluated at $w=1$ and $w=0$, respectively, as given below.

$APE_{w}=E\left[G(\mathbf{\beta}_{0}+\beta_{w}+\mathbf{z}(\mathbf{\beta}_{z}+\mathbf{\beta}_{wz}))-G(\mathbf{\beta}_{0}+\mathbf{z}\mathbf{\beta}_{z})\right],\vspace{-3mm}$

where $G()$ is either the identity function (for the LPM estimated by OLS),
the probit CDF, the logit CDF, or the ramp function. We also use the
nonparametric kernel estimator from our simulations, treating “white” as
discrete but not imposing a linear index or specific latent error
distribution.

Table 10 presents the results. Using the LPM estimated by OLS, about 18% of
observations have predicted probabilities outside the unit interval, so the
H-O results clearly imply OLS is inconsistent for the slope parameters if the
ramp model is correct. There is little reason to expect OLS will approximate
this APE, either based on the theoretical results of Stoker, (1986) or our
simulation study. Many of the explanatory variables are binary, and the
continuous variables (e.g., income) tend to be skewed. For each variable,
normality is strongly rejected by a Jarque-Bera test (a joint test of the
skewness and kurtosis) with p-values well below 1%. The model also includes
interactions between the continuous variables and a binary variable. Using the
LPM estimates, the APE for $white$ is $5.3$ percentage points and it is only
marginally significant. Using the nonlinear parametric estimators, the APE are
each a bit larger at about $7.0$ percentage points, and they are all
significant at the 1% level. Using the nonparametric estimator, the point
estimate is quite a bit larger at 17.2 percentage points, but it is only
marginally significant due to a very large standard error.

Table 10: Estimates of the APE of “White” on Loan Approval ($N=1976$)
---
| LPM | Ramp | Probit | Logit | Nonpar.
| (OLS) | (NLS) | (QMLE) | (QMLE) | (LL)
Estimate | 0.0532 | 0.0706 | 0.0695 | 0.0712 | 0.1721
Robust SE | 0.0278 | 0.0227 | 0.0220 | 0.0219 | 0.0955
$P(0\leq\hat{y}\leq 1)$ | 0.8173 | 0.6027 | 1.0000 | 1.0000 | 0.6587
Mean Squared Error | 0.1171 | 0.0839 | 0.0840 | 0.0837 | 0.1804

  * •

Note: There were only 1976 complete cases out of 1989 observations total. The
nonparametric regression used only 1975 observations. All robust standard
errors were computed using the sandwich forms and the delta method, save the
nonparametric regression for which we used a nonparametric bootstrap with 500
replications.

Interestingly, OLS predicts only 18% of observations with indexes outside the
unit interval, whereas NLS predicts nearly 40%, which follows the pattern of
many of our simulations from the previous section and suggests trimming the
sample once is not sufficient to consistently estimate the parameters or APEs
under the piecewise linear model. Model selection by the minimum mean squared
error favors logit, though the other nonlinear models are very similar.

##  7 Implications for Empirical Research

We have revisited the conclusions reached by Horrace and Oaxaca, (2006)
concerning the ability of the linear projection parameters—consistently
estimated by OLS—to recover interesting parameters. We argue that H-O’s focus
on the parameters in the underlying index model is misguided; instead, one
should focus on the APEs. Focusing on the APEs is hardly controversial, as
almost every study that employs any model nonlinear in the explanatory
variables reports estimated APEs.

Once the focus is on the APEs, a few useful conclusions emerge in an expanded
version of the H-O model, which allows for varying support in the underlying
uniform distribution. First, when the explanatory variables have a
multivariate normal distribution, the LP parameters are identical to the
population APEs under a general index model. Importantly, this is true even
when the flat parts of the ramp function occur with high probability. In this
case, the LP parameters, $\gamma_{j}$, will be greatly attenuated toward zero
compared with the index parameters, $\beta_{j}$. The logit and probit models,
estimated by quasi-MLE (because the response probabilities are misspecified),
also approximate the APEs very well. Nonlinear least squares estimation of the
ramp function is a new option, and we have shown the estimator is consistent
for the best MSE approximation and asymptotically normal.

When the explanatory variables have asymmetric distributions, the conclusions
for OLS are not as sanguine—unless the support of $\mathbf{x\beta}$ is
contained entirely in the support $\left[-a,a\right]$ of the uniform
distribution of the the latent error. Some simulations show that even if the
probability of $\mathbf{x\beta}$ being in the unit interval is high (e.g. 93%
in Table 7), the LP parameters are not very close to the true APEs. Especially
when the support $\left[-a,a\right]$ is narrow, the logit and probit
approximations to the APEs can be notably better than those for OLS.

To summarize, in evaluating different strategies, we need to make sure we have
carefully defined the population quantities of interest, and then we make
proper comparisons across different approaches. OLS estimation of the LPM has
good finite sample properties for the APE in many cases when the covariates
are symmetrically distributed. Probit, Logit, and the ramp model continue to
have good finite sample properties for estimating the APEs when the covariates
are asymmetrically distributed. Using the fraction of estimated response
probabilities in $\left[0,1\right]$ is neither necessary nor sufficient for
good performance of OLS. In the $a=1/4$ case with interactions, the estimated
probability is almost $0.95$ but OLS has a severe bias toward zero for the
APEs. By contrast, the three nonlinear models show very little bias. A
nonlinear model, of course, offers other advantages over the LPM, such as more
realistic response probabilities and nonconstant partial effects. However,
when the APEs are of interest, OLS is more widely applicable than a simple
reading of H-O might suggest.

The conclusions drawn here are easily extended to the case where $y$ is a
fractional response, where the limit values zero and one can occur with
positive probability. In particular, Stoker (1986) can be applied to
$E\left(y|\mathbf{x}\right)$. If this conditional mean follows the same ramp
function, the qualitative conclusions obtained in the binary case will remain.

## References

  * Angrist and Pischke, (2009) Angrist, J. D. and Pischke, J.-S. (2009).  Mostly harmless econometrics: An empiricist’s companion.  Princeton University Press.
  * Horowitz and Savin, (2001) Horowitz, J. L. and Savin, N. (2001).  Binary response models: Logits, probits and semiparametrics.  Journal of Economic Perspectives, 15(4):43–56.
  * Horrace and Oaxaca, (2006) Horrace, W. C. and Oaxaca, R. L. (2006).  Results on the bias and inconsistency of ordinary least squares for the linear probability model.  Economics Letters, 90(3):321–327.
  * Hunter and Walker, (1996) Hunter, W. C. and Walker, M. B. (1996).  The cultural affinity hypothesis and mortgage lending decisions.  The Journal of Real Estate Finance and Economics, 13:57–70.
  * Newey and McFadden, (1994) Newey, W. K. and McFadden, D. (1994).  Large Sample Estimation Testing.  Handbook of Econometrics, 4:2113–2245.
  * Ruud, (1983) Ruud, P. A. (1983).  Sufficient Conditions for the Consistency of Maximum Likelihood Estimation Despite Misspecification of Distribution in Multinomial Discrete Choice Models.  Econometrica, 51(1):225–228.
  * Stoker, (1986) Stoker, T. M. (1986).  Consistent Estimation of Scaled Coefficients.  Econometrica, 54(6):1461–1481.
  * van den Berg and Siflinger, (2022) van den Berg, G. J. and Siflinger, B. M. (2022).  The effects of a daycare reform on health in childhood–Evidence from Sweden.  Journal of Health Economics, 81:102577.
  * Wooldridge, (2010) Wooldridge, J. M. (2010).  Econometric analysis of cross section and panel data.  MIT Press.
  * Wooldridge, (2019) Wooldridge, J. M. (2019).  Introductory Econometrics: A Modern Approach.  Cengage Learning.

## Appendix

### Proof of Theorem 2

$Proof.$ We will obtain the asymptotic normality of the NLS estimator by
applying Theorem 7.1 of Newey and McFadden (1994). Condition (i) and (ii) of
Theorem 7.1 follows from our assumptions. As we discussed in the main context,
condition (iii) is satisfied as long as $\mathbf{x}$ contains a continuous
variables $x_{j}$ with nonzero $\beta_{jo}$ so that
$P(\mathbf{x}_{i}\mathbf{\beta}_{o}=0\text{ or
}\mathbf{x}_{i}\mathbf{\beta}_{o}=1)=0$.

For condition (iv), notice that the first derivative of the object function is
well defined at $\mathbf{\beta}_{o}$ with probability 1:

$\displaystyle
D_{N}(\mathbf{\beta}_{o})=\nabla_{\beta}Q_{N}(\mathbf{\beta}_{o})=\frac{1}{N}\sum_{i=1}^{N}\mathbf{x}_{i}^{\prime}(y_{i}-\mathbf{x}_{i}\mathbf{\beta}_{o})1\\{\mathbf{x}_{i}\mathbf{\beta}_{o}\in(0,1)\\}=\frac{1}{N}\sum_{i=1}^{N}\mathbf{x}_{i}^{\prime}u_{i}1\\{\mathbf{x}_{i}\mathbf{\beta}_{o}\in(0,1)\\},\vspace{-3mm}$

where $u_{i}=y_{i}-R(\mathbf{x}_{i}\mathbf{\beta}_{o})$. Since
$E\|\mathbf{x}_{i}^{\prime}u_{i}\|1\\{\mathbf{x}_{i}\mathbf{\beta}_{o}\in(0,1)\\}<\infty$
under the assumption $E\|x\|^{2}<\infty$, the vector Lindberg-Levy CLT
applies:

$\displaystyle\sqrt{N}D_{N}(\mathbf{\beta}_{o})\stackrel{{\scriptstyle
d}}{{\to}}\mathbb{N}\left(0,\mathbf{\Omega}(\mathbf{\beta}_{o})\right),\vspace{-3mm}$

giving condition (iv). Lastly, for condition (v), following Newey and McFadden
(1994), we can rewrite

$\displaystyle\sqrt{N}[Q_{N}(\mathbf{\beta})-Q_{N}(\mathbf{\beta}_{o})]$
$\displaystyle=$
$\displaystyle\sqrt{N}\left[D_{N}(\mathbf{\beta}_{o})(\mathbf{\beta}-\mathbf{\beta}_{o})+Q(\mathbf{\beta})-Q(\mathbf{\beta}_{o})\right]+\|\mathbf{\beta}-\mathbf{\beta}_{o}\|M_{N}(\mathbf{\beta}),\vspace{-3mm}$

where $M_{N}(\mathbf{\beta})$ is the remainder term, defined as:

$\displaystyle
M_{N}(\mathbf{\beta})=\frac{\sqrt{N}\left[Q_{N}(\mathbf{\beta})-Q_{N}(\mathbf{\beta}_{o})-D_{N}^{\prime}(\mathbf{\beta}_{o})(\mathbf{\beta}-\mathbf{\beta}_{o})-\left(Q(\mathbf{\beta})-Q(\mathbf{\beta}_{o})\right)\right]}{\|\mathbf{\beta}-\mathbf{\beta}_{o}\|}.\vspace{-3mm}$

Since $D_{N}(\mathbf{\beta}_{o})$ is the gradient of $Q_{N}(\mathbf{\beta})$
at $\mathbf{\beta}_{o}$,
$Q_{N}(\mathbf{\beta})-Q_{N}(\mathbf{\beta}_{o})-D_{N}(\mathbf{\beta}_{o})(\mathbf{\beta}-\mathbf{\beta}_{o})$
goes to zero faster than $\|\mathbf{\beta}-\mathbf{\beta}_{o}\|$ as
$\mathbf{\beta}$ goes to $\mathbf{\beta}_{o}$, by the definition of the
gradient. Similarly, due to
$\nabla_{\beta}Q(\mathbf{\beta}_{o})=E\left(s_{i}(\mathbf{\beta}_{o})\right)=0$,
$Q(\mathbf{\beta})-Q(\mathbf{\beta}_{o})$ goes to 0 faster than
$\|\mathbf{\beta}-\mathbf{\beta}_{o}\|$ as $\mathbf{\beta}$ goes to
$\mathbf{\beta}_{o}$. Under the moment conditions, we can easily show
$Q(\mathbf{\beta})-Q(\mathbf{\beta}_{o})\to 0$ in probability and so
$\sqrt{N}[Q_{N}(\mathbf{\beta})-Q(\mathbf{\beta})]$ is bounded in probability
for each $\beta$. Also note that $\sqrt{N}D_{N}(\mathbf{\beta}_{o})$ is
bounded in probability due to the asymptotic normality. Since the numerator is
bounded in probability and converges to zero faster than the denominator, we
conclude that for any $\varepsilon_{N}>0$,
$\lim_{N\to\infty}sup_{\|\mathbf{\beta}-\mathbf{\beta}_{o}\|<\varepsilon_{N}}M_{N}(\mathbf{\beta})\to
0$ in probability, which implies condition (v). ∎

### Proof of Theorem 3

$Proof.$ Consider $\mathbf{\Omega}(\mathbf{\hat{\beta}})$:

$\displaystyle\mathbf{\Omega}(\mathbf{\hat{\beta}})=$
$\displaystyle\frac{1}{N}\sum_{i=1}^{N}\mathbf{x}_{i}^{\prime}\mathbf{x}_{i}\left(y_{i}-R(\mathbf{x}_{i}\hat{\mathbf{\beta}})\right)^{2}1\\{\mathbf{x}_{i}\mathbf{\hat{\beta}}\in(0,1)\\}$
$\displaystyle\equiv$
$\displaystyle\frac{1}{N}\sum_{i=1}^{N}a(\mathbf{x_{i}},\hat{\mathbf{\beta}})\vspace{-3mm}$

Note that $E|y_{i}-R(\mathbf{x}_{i}\mathbf{\beta})|^{4}\leq 1$ for any
$\mathbf{\beta}\in\mathcal{B}$ since both $y_{i}$ and $R(.)$ are naturally
bounded in $[0,1]$ with probability 1. Then, we have

$\displaystyle
E\sup_{\mathbf{\beta}\in\mathcal{B}}\|a(x,\mathbf{\beta})\|\leq(E\left\|\mathbf{x}_{i}\right\|^{4}E|y_{i}-R(\mathbf{x}_{i}\mathbf{\beta})|^{4})^{1/2}<\infty,\vspace{-3mm}$

where the first inequality follows from Hölder’s inequality. Also note that
$a(\mathbf{x_{i}},\mathbf{\beta})$ is continuous at $\mathbf{\beta}_{o}$ with
probability one given that
$P(\mathbf{x}_{i}\mathbf{\beta}_{o}=0)=P(\mathbf{x}_{i}\mathbf{\beta}_{o}=1)=0$.
Then, we can apply Lemma 4.3 of Newey and McFadden, (1994):

$\displaystyle\mathbf{\Omega}(\mathbf{\hat{\beta}})=\frac{1}{N}\sum_{i=1}^{N}a(\mathbf{x_{i}},\hat{\mathbf{\beta}})\stackrel{{\scriptstyle
p}}{{\to}}E(a(\mathbf{x_{i}},\mathbf{\beta}_{o}))=\mathbf{\Omega}(\mathbf{\beta}_{o}).\vspace{-3mm}$

Similarly, Lemma 4.3 also applies to
$\mathbf{A}_{N}(\mathbf{\hat{\beta}})=\frac{1}{N}\sum_{i=1}^{N}\mathbf{x}_{i}^{\prime}\mathbf{x}_{i}1\\{\mathbf{x}_{i}\mathbf{\hat{\beta}}\in(0,1)\\}$:

$\displaystyle\operatorname*{plim}_{N\to\infty}\frac{1}{N}\sum_{i=1}^{N}1\\{\mathbf{x}_{i}\mathbf{\hat{\beta}}\in(0,1)\\}\mathbf{x}_{i}^{\prime}\mathbf{x}_{i}=\mathbf{A}(\mathbf{\beta}_{o}).\vspace{-3mm}$

So, we conclude that

$\displaystyle\mathbf{\hat{V}}=\mathbf{A}_{N}(\hat{\mathbf{\beta}})^{-1}\mathbf{\Omega}_{N}(\hat{\mathbf{\beta}})\mathbf{A}_{N}(\hat{\mathbf{\beta}})^{-1}\stackrel{{\scriptstyle
p}}{{\to}}\mathbf{A}(\mathbf{\beta}_{o})^{-1}\mathbf{\Omega}(\mathbf{\beta}_{o})\mathbf{A}(\mathbf{\beta}_{o})^{-1}.\,\,\qed{}\vspace{-3mm}$

Additional Simulations

Table 11. No interaction, $x_{1}$ sym, $x_{2}$ sym. binary, $u\sim N(0,1)$
---
$N=1000$ | Simulated | LPM | Ramp | Probit | Logit | Nonpar.
|  | Truth* | (OLS) | (NLS) | (QMLE) | (QMLE) | (LL)
$APE_{1}$ | mean | 0.0782 | 0.0774 | 0.0774 | 0.0774 | 0.0774 | 0.0776
sd | 0.0001 | 0.0171 | 0.0171 | 0.0171 | 0.0171 | 0.0175
$APE_{2}$ | mean | -0.1168 | -0.1160 | -0.1160 | -0.1154 | -0.1152 | -0.1160
sd | 0.0001 | 0.0339 | 0.0339 | 0.0336 | 0.0335 | 0.0344
$P(y=1)$ | 0.4805 |  |  |  |  |
$P(0\leq\widehat{y}\leq 1)$ |  | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 0.9983

  * •

*This column contains the empirical means and standard deviations of the sample average partial effects at the true parameter values.

Table 12. With interaction, $x_{1}$ sym, $x_{2}$ sym. binary, $u\sim N(0,1)$
---
$N=1000$ | Simulated | LPM | Ramp | Probit | Logit | Nonpar.
|  | Truth* | (OLS) | (NLS) | (QMLE) | (QMLE) | (LL)
$APE_{1}$ | mean | 0.0199 | 0.0191 | 0.0191 | 0.0191 | 0.0191 | 0.0193
sd | 0.0019 | 0.0171 | 0.0170 | 0.0171 | 0.0171 | 0.0176
$APE_{2}$ | mean | -0.1179 | -0.1177 | -0.1177 | -0.1173 | -0.1172 | -0.1178
sd | 0.0037 | 0.0337 | 0.0337 | 0.0336 | 0.0336 | 0.0347
$P(y=1)$ | 0.4563 |  |  |  |  |
$P(0\leq\widehat{y}\leq 1)$ |  | 0.9999 | 0.9999 | 1.0000 | 1.0000 | 0.9983

  * •

*This column contains the empirical means and standard deviations of the sample average partial effects at the true parameter values.

Table 13. No interaction, $x_{1}$ asym, $x_{2}$ asym. binary, $u\sim N(0,1)$
---
$N=1000$ | Simulated | LPM | Ramp | Probit | Logit | Nonpar.
|  | Truth* | (OLS) | (NLS) | (QMLE) | (QMLE) | (LL)
$APE_{1}$ | mean | 0.0722 | 0.0410 | 0.0667 | 0.0728 | 0.0753 | 0.0672
sd | 0.0005 | 0.0077 | 0.0101 | 0.0094 | 0.0099 | 0.0134
$APE_{2}$ | mean | -0.1070 | -0.0840 | -0.1064 | -0.1075 | -0.1083 | -0.1092
sd | 0.0007 | 0.0322 | 0.0305 | 0.0296 | 0.0294 | 0.0315
$P(y=1)$ | 0.6143 |  |  |  |  |
$P(0\leq\widehat{y}\leq 1)$ |  | 0.9873 | 0.9659 | 1.0000 | 1.0000 | 0.9844

  * •

*This column contains the empirical means and standard deviations of the sample average partial effects at the true parameter values.