A Kernelization-Based Approach to Nonparametric Binary Choice Models

A Kernelization-Based Approach to Nonparametric Binary Choice Models

\singlespacing

\singlespacing JEL Classification: C13, C14, C25, C81 \\ Keywords and phrases: nonparametric, binary choice models, reproducing kernel Hilbert space, sieve estimation, SNP

\setcounter{footnote}{0} \onehalfspacing

Introduction

Binary choice problems arise widely in economics. Examples include an individual's choice to work or not, a firm's decision to enter a market, and a household's intention to migrate. Binary choice models (BCMs) are extensively used to analyze these problems due to its underlying microeconomic interpretations of latent utility maximization (e.g., bhattacharya-21). Typically, the latent utility—net utility of one choice over another—comprises two components: a systematic component, which is a deterministic function $G$ of covariates, and a random component $\varepsilon$ representing idiosyncratic error.

To avoid possible misspecification and resulting inconsistency for the estimation, matzkin-92 first studied fully nonparametric BCMs, allowing for nonparametric $G$ and $\varepsilon$. However, the proposed estimator is not practical, as it relies on maximizing empirical likelihood under constraints, e.g., monotonicity of $\varepsilon$'s cumulative distribution function (CDF), whose computation becomes intractable as the number of regressors or sample size increases. Moreover, this approach cannot be used to estimate partial effects, which are important policy parameters.

One might consider using sieve approximations for nonparametric functions. However, allowing for a nonparametric $G$ can pose computational challenges when dealing with multiple covariates, even if the error distribution is known. This is because estimating nonlinear models requires numerical optimization, and sieve approximations of $G$ can require a large number of basis functions.\footnote{For example, the polynomial expansions of 50 variables up to the 2nd and 3rd orders produce 1,325 and 23,425 basis functions, respectively. } Although there is recent work addressing computational concerns in linear index models (e.g., ahn-ichimura-powell-ruud-18; khan-lan-tamer-21), an additional challenge here lies in handling a nonparametric $G$, which was assumed to be linear in these studies.

In light of these practical challenges, we propose a new estimation method for a broad class of nonparametric BCMs. We approximate the nonparametric component of covariates using functions in a reproducing kernel Hilbert space (RKHS), which can be viewed as a special sieve space, and couple it with further regularization through spectral cutoff for dimensional reduction. For the nonparametric error component, we follow gallant-nychka-87 and approximate its density by squared Hermite polynomials, resulting in simple closed-form approximate CDFs that can be easily evaluated without numerical integration.

We highlight the key computational differences between using classical sieve choices (e.g., polynomials or splines) and using RKHS as special sieves.\footnote{The notations used here are temporary for illustration, with formal results presented in Section (ref).} An estimator for a nonparametric function $G$ is often obtained by optimizing over a set of functions with certain basis functions, either by maximizing likelihood or minimizing least squares. Conventional sieve methods typically optimize over the coefficients of basis functions, which can become high-dimensional with multiple covariates. In contrast, when optimizing $G$ over an RKHS $\mathbb G_k$ with reproducing kernel $k(\cdot,\cdot)$, the estimator takes a different form. Here, $\mathbb G_k$ consists of functions spanned by $\{k(x,\cdot): x\in \mathbb R^d\}$, with the inner product induced by $\langle k(s,\cdot), k(t,\cdot) \rangle_{\mathbb{G}_k} = k(s,t)$; see Appendix (ref) for a brief introduction to RKHSs and further references. Specifically, for observed covariates $X_1,\dots,X_n$, the optimization over $\mathbb G_k$ effectively reduces to optimizing over the coefficients of $k(\cdot,X_i)$'s.\footnote{This follows from the representer theorem; see, e.g., Theorem 4.2 in scholkopf-smola-02, and the references therein for the history of the development of the representer theorem. } Notably, the number of these coefficients is independent of the covariate dimension.

To balance the bias-variance trade-off, we use RKHS balls with radii that increase to infinity at certain rates of the sample size. This radius constraint simplifies to a quadratic constraint in optimization. The optimization over the coefficients of $k(\cdot,X_i)$'s can be computationally challenging when the sample size $n$ is large. To address this, we employ spectral cutoff regularization on the $n\times n$ matrix given by $k(X_i,X_j)$ to further reduce the dimensionality for optimization, which is particularly convenient in our setting. We provide an upper bound on the difference between the objective function values at the optima with and without regularization. In our theory, this difference is assumed to vanish asymptotically, allowing the spectral cutoff regularized estimator to be considered as a near-optimal solution to the original problem.

A key theoretical contribution of this paper is a simple perspective of viewing RKHS balls as special sieve spaces—an idea that can be applied in many nonparametric problems beyond BCMs, and allows for RKHS-based methods to be seamlessly integrated into existing sieve estimation frameworks (e.g., chen-07). It is especially helpful when there are multiple covariates and classical sieves are practically challenging. Furthermore, our approach is more robust to misspecification compared to recent literature on RKHS-based methods in econometrics (e.g., singh-22; singh-xu-gretton-24), which typically assumes that the true nonparametric function belongs to a specific RKHS—an assumption that can be restrictive or misspecified.\footnote{ E.g., assuming that the true function lies in the RKHS with a Gaussian kernel—one of the most widely used kernels in practice—requires the function to be infinitely differentiable, which may be overly restrictive. } In our theoretical framework, we impose standard smoothness conditions on the differentiability and boundedness of $G$, as is common in the nonparametric literature, and appropriately choose an RKHS so that $G$ can be approximated by functions within RKHS balls. By explicitly accounting for the approximation error rate that arises when approximating a smooth function using elements from RKHS balls, our approach ensures robustness to misspecification when the true function does not belong to a prespecified RKHS.\footnote{ Our theoretical analysis builds on approximation results for RKHSs (e.g., steinwart-01, micchelli-xu-zhang-06) and on bounds for approximation error and entropy numbers of Gaussian RKHS balls (e.g., smale-zhou-03, kuhn-11). }

Our proposed method is not only computationally effective but also theoretically sound. We show the consistency of the proposed kernelized non-parametric (KNP) estimator for both the systematic component and the distribution function of the random component. The KNP estimation procedure provides a natural plug-in estimator for the conditional choice probability (CCP) function, for which we establish the convergence rate.

The KNP approach is useful for estimating important parameters of policy interest, including average partial effects (APEs) and, when accounting for heterogeneity, conditional APEs. Both APEs and conditional APEs are special cases of weighted average derivative functionals of the CCP. We establish the asymptotic normality of the estimators for weighted average derivatives. Moreover, these estimators are easy to compute, with a computational procedure that remains unchanged regardless of the covariate dimension.

The effectiveness of the KNP estimator is demonstrated using extensive simulation studies. We find that, compared to parametric estimation methods, the proposed method effectively improves the finite sample performance in case of misspecification and has rather mild efficiency loss if the model is correctly specified. To demonstrate the practical use of our proposed method, we revisit an empirical application of heyes-saberian-19,heyes-saberian-22, examining the effect of outdoor temperature on court judges' decisions.

\paragraph{Outline} The rest of the paper is organized as follows. Section (ref) describes the model. Section (ref) defines the proposed estimator and describes its implementation. Section (ref) presents the asymptotic properties. Simulation studies are in Section (ref). In Section (ref), the KNP estimation procedure is applied in a model on judges’ decisions and outdoor environments. Section (ref) concludes the paper. All of the proofs and other technical details are collected in the Appendix. Programs for implementation, along with replication packages for the simulation studies and the empirical application, are available online on the author’s webpage.

The Model

We consider the BCM that generates the binary outcome variable $Y \in \{0,1\}$ as follows.

where $G_0$ is an unknown function of covariates $X \in \mathbb{R}^{d_x}$, and $\varepsilon$ is the idiosyncratic error term. The conditional choice probability (CCP) is therefore

Let $F_0$ and $f_0$ denote the CDF and Lebesgue density of $\varepsilon$. Let $\mathcal{X} \subset \mathbb{R}^{d_x}$ be the support of $X$.\footnote{We follow the conventional definition that the support of a random vector $Z$ with distribution $P_Z$ is the smallest closed set $A$ which satisfies $P_Z(A)=1$. See, e.g., Page 181 in billingsley-95 for the existence and uniqueness of the support.}

For the identification of $G_0$ and $F_0$, we assume that one component of $X$, $V$, which has large support, enters $G_0$ linearly and is separable from the other components, $W$. This assumption is more general than that of many parametric (e.g., probit or logit) or semiparametric BCMs (e.g., manski-75,manski-85), which typically impose a fully linear form on $G_0$. We let $X = (V,W')'$, with $\mathcal V$ and $\mathcal W$ denoting the supports of $V$ and $W$, respectively. The assumptions for the identification of $G_0$ and $F_0$ are as follows.

Assumption (ref)((ref)) requires that $V$ has support $\mathbb R$ both marginally and conditional on $W=w_\ast$.\footnote{Note that $\mathcal{L}(V|W=w_\ast)$ having support $\mathbb R$ does not ensure that $\mathcal L(V)$ has support $\mathbb R$: If for any $w\neq w_\ast$ the support of $\mathcal L(V|W=w)$ is contained in a fixed bounded set, and $\mathbb P\{W=w_\ast\} = 0$, then the support of $\mathcal L(V)$ is bounded. } When $V$ is independent of $W$, Assumption (ref)((ref)) reduces to requiring that the support of $V$ is $\mathbb R$. Moreover, Assumption (ref)((ref)) requires that $\varepsilon$ also has support $\mathbb{R}$.

We denote the true parameter by $\theta_0 = (g_0,F_0)$, and let $\Theta=\mathcal G\times \mathcal F$. Under the independence of $\varepsilon$ from $X$, we define the CCP given by $\theta = (g,F) \in \Theta$ as $p_{\theta}(x) = F(v+g(w))$, where $x = (v,w')'$. We also write $p_{0} = p_{\theta_0}$ the true CCP, i.e. $p_0(x) = F_0(v + g_0(w) )$.

For a criterion function $\ell$ given by either

or

where $z=(y,x')'$, the following theorem gives the identification of $\theta_0$ in the sense that $\theta_0 \in \Theta$ is the unique minimum of the population objective function $Q(\theta) = \mathbb{E} \ell(Z,\theta)$.

Here, the uniqueness is in the sense that, for any $\theta\in \Theta$ which minimizes $\theta \mapsto \mathbb{E} \ell(Z,\theta)$, it must hold that $g(w) = g_0(w)$ for any $w\in \mathcal W$ and $F(u) = F_0(u)$ for any $u\in \mathbb R$.

The identification argument follows the general framework in matzkin-92,matzkin-93,matzkin-94 for nonparametric binary choice models. Our contribution is not a new identification strategy. To the best of our knowledge, however, the specific separable-index structure we consider is not stated in exactly this form in matzkin-92,matzkin-93,matzkin-94. For completeness, we state the identification result explicitly and provide a self-contained proof by adapting her proof steps.

A word on notation. For a function $f$ whose domain is a subset of $\mathbb{R}^d$, let \[ D^{(\lambda)} f(x) = \frac{\partial^{\lambda_1}}{\partial x_1^{\lambda_1}} \frac{\partial^{\lambda_2}}{\partial x_2^{\lambda_2}} \cdots \frac{\partial^{\lambda_d}}{\partial x_d^{\lambda_d}}, \] where $\lambda = (\lambda_1,\lambda_2,\cdots,\lambda_d)'$ and its elements are nonnegative integers. For such multi-index $\lambda$, let $|\lambda| = \sum_{j=1}^{d} \lambda_j$. Let $D^{(0)} f = f$. For functions whose domain is a subset of $\mathbb{R}$, the $\lambda$-th derivative is denoted as $f^{(\lambda)}$ for any nonnegative integer $\lambda$, and let $f^{(0)} = f$. We write $P_W$ and $P_X$ for the distribution of $W$ and $X$.

Kernelized Non-Parametric Estimator

Let $\{Y_i,X_i\}_{i=1}^n$ denote $n$ independent observations on the dependent variable $Y$ and covariate vector $X = (V,W')'$, and let $Z_i =(Y_i,X_i')'$. In this section, we first define the proposed estimator in Section (ref), followed by the discussions for implementations in Section (ref), and practical implementation procedure in Section (ref).

Kernelized Non-Parametric Estimator

Motivated by Theorem (ref), we propose an estimator, which will be referred to as kernelized non-parametric (KNP) estimator, for $\theta_0 = (g_0,F_0)$ as follows.

The KNP estimator $\hat \theta = (\hat g, \hat F)$ is given by

where we choose the least squares loss\footnote{Here, we choose the least squares loss function in (ref), as its boundedness properties facilitate the proofs. The MLE objective function in (ref) could also be used with additional conditions, including assumptions controlling the tails of $\log p_{\theta}(x), \log (1- p_{\theta}(x))$ over $\theta \in \Theta, x\in \mathcal X$. In addition, for the plug-in estimator of the weighted average partial derivative, the MLE objective will lead to the same asymptotic distribution as in Theorem (ref), under an analogous set of conditions to Assumption (ref). }

and the sets $\mathcal G_n, \mathcal F_n$ for optimization are defined in (ref) and (ref) below.

\paragraph{Set $\mathcal G_n$} To describe the set $\mathcal G_n$, we first briefly introduce the RKHS notation. Let a kernel $k: \mathcal W \times \mathcal W \to \mathbb{R}$ be a symmetric function which is positive definite, in the sense that $\sum_{i=1}^N \sum_{j=1}^N a_i a_j k(s_i, s_j) \geq 0$ for any $a_i \in \mathbb{R}$, $s_i \in \mathcal W $, $i =1,\cdots,N$, and any positive integer $N$. Let $\mathbb G_k$ be the reproducing kernel Hilbert space (RKHS) with reproducing kernel $k$, defined as the completion of the linear span of $\{ k(\cdot,w) |w\in \mathcal W \}$ under the RKHS norm $\|\cdot\|_{\mathbb{G}_k}$, where $\|g\|_{\mathbb G_k} = \sqrt{\langle g,g\rangle_{\mathbb G_k}}$ is induced by the inner product $\left\langle \sum_{i=1}^N a_i k(\cdot,s_i), \sum_{j=1}^M b_j k(\cdot,t_j) \right \rangle_{\mathbb{G}_k} := \sum_{i=1}^N \sum_{j=1}^M a_i b_j k(s_i,t_j)$ for any $s_i,t_j\in \mathcal W, a_i,b_j \in \mathbb{R}$ and any integers $M, N$. In particular, $\|\sum_{i=1}^N a_i k(\cdot,s_i) \|_{\mathbb G_k}^2 = \sum_{i=1}^N \sum_{j=1}^N a_i a_j k(s_i, s_j) $. An example of a commonly used kernel $k$ is the class of Gaussian kernels: For a $\sigma>0$, \[ k(s,t) = \exp\!\left( - \frac{\|s-t\|^2}{2\sigma^2} \right), \quad s,t \in \mathcal W. \] For any prespecified $\sigma^2>0$, the resulting RKHS $\mathbb{G}_k$ is dense in the space of all continuous functions on $\mathcal{W}$ under the uniform norm.

With these notations, the set for the optimization of $g$ is chosen so that $g_0:\mathcal W \to \mathbb R$ with $g_0(w_\ast)=0$ is approximated based on functions within the balls of a RKHS $\mathbb G_k$. Specifically,

where $B_n$ is the radius of the RKHS ball. The form $g(w) = \tilde g(w) - \tilde g(w_\ast)$ is used to ensure that $g(w_\ast)=0$, which is location normalization for identification. Intuitively, one can view $\tilde g_0(\cdot)$ (so that $g_0(\cdot) = \tilde g_0(\cdot) - \tilde g_0(w_\ast)$) as being approximated by a linear combination of kernel basis functions from the class $\{k(w,\cdot)| w\in \mathcal W\}$, and the constraint $\|\tilde g\|_{\mathbb G_k} \leq B_n$ then controls the “roughness” of the approximating function and helps prevent over-fitting.

In practice, we may choose $k$ to be a Gaussian kernel. Moreover, $w_\ast$ is specified as a point in $\mathcal{W}$. In practice, it is convenient to set $w_\ast$ to zero, coupled with standardizing the observations $(W_i)_{i=1}^n$ to have zero mean by subtracting their sample average.

The radius $B_n$ governs the bias-variance tradeoff when estimating $g_0$ and is required to grow to infinity as $n\to \infty$ to ensure the consistency of the estimator. A larger $B_n$ reduces approximation error, since $g_0$ can be better approximated by functions in $\mathcal G_n$, but also makes $\hat g$ more variable because it is optimized over a larger class $\mathcal G_n$. A theoretically optimal rate for $B_n$, established later in Corollary (ref), depends on unknown constants, including the number of uniformly bounded derivatives of $g_0$ that exist. In practice, $B_n$ can be chosen via multi-fold cross-validation.

\paragraph{Set $\mathcal F_n$} Following fenton-gallant-96, fenton-gallant-96-joe, the class of functions used to approximate $F_0$ is

where $J_n$ is some positive integer growing to infinity as $n\to \infty$. The idea behind $\mathcal F_n$ is that the density of $F_0$ is approximated by $f(u; \tau) = \big( e^{-u^2/4} \sum_{j=0}^{J_n} \tau_j u^j \big)^2$, which is effectively the square of the product of the density of $N(0,2)$ and a polynomial of order $J_n$.

The polynomial order $J_n$ governs the bias-variance tradeoff when estimating $F_0$ and is required to grow to infinity as $n\to \infty$ to ensure the consistency of the estimator. A larger $J_n$ reduces approximation error, since $F_0$ can be better approximated by CDFs in the larger class $\mathcal F_n$, but the estimator $\hat F$ becomes more variable. A theoretically optimal rate for $J_n$, as established later in Corollary (ref), depends on unknown constants, including the number of derivatives of $F_0$'s density $f_0$ that exist. In practice, $J_n$ can be chosen via multi-fold cross-validation.

Explanations for the Implementation and the Estimator

Now we discuss how we optimize over $g\in \mathcal{G}_n$ and over $F \in \mathcal{F}_n$, in order to obtain an estimator $(\hat g, \hat F)$. As in Remark (ref), it suffices to find one solution that solves (ref), either exactly or approximately.

From optimization over RKHS to a finite-dimensional problem

Since $\mathcal{G}_n$ is infinite-dimensional, the optimization over $g\in \mathcal{G}_n$ is not directly solvable in practice. The following Proposition (ref) ensures that we can find a solution to $g$ in (ref) in a finite-dimensional subspace, by setting

and optimizing over $\delta = (\delta_0, \delta_1,\dots,\delta_n)'$, a finite-dimensional Euclidean space. Then the constraint $\|\tilde g\|_{\mathbb{G}_k} \leq B_n$ in the definition of $\mathcal G_n$ can be imposed via $\delta' K\delta \leq B_n^2$, where $\delta :=(\delta_0, \delta_1,\dots,\delta_n)'$ and $K$ is the $(n+1)\times (n+1)$ square matrix whose elements are given by $k(W_i,W_j)$ for $i,j=0,1,\dots,n$. This is because $\big\| \sum_{j=0}^n \delta_j k(W_j, \cdot) \big\|_{\mathbb G_k}^2 = \sum_{i,j=0}^n \delta_i \delta_j k(W_i,W_j)$.

Proposition (ref) shows that, for any $(g,F)\in \mathcal G_n\times \mathcal F_n$, the criterion value $\hat Q( (g,F))$ can be attained by some $(g_\ast,F)$ where $g_\ast(\cdot) = \sum_{j=0}^n \delta_j \big( k(W_j,\cdot) - k(W_j,w_\ast) \big) \in \mathcal G_n$. Therefore, when minimizing $\hat Q( (g,F))$ over $\mathcal G_n\times \mathcal F_n$, we can restrict, without loss of generality, the search for $g$ to a finite-dimensional class of functions with the form in (ref). This restriction leaves the minimum value of the criterion unchanged and yields a valid solution to (ref).

Therefore, a solution $\hat\theta$ in (ref) is given by

where evaluations of $F(\cdot;\tau)$ are computed using the closed-form expression given in Appendix (ref) for any $\tau =(\tau_1,\dots,\tau_{J_n})'\in\mathbb R^{J_n}$, and $(\hat \delta, \hat \tau)$ solves\footnote{Note that $g = \tilde g - \tilde g(w_\ast)$ for some $\tilde g(w) =\sum_{j=0}^n \delta_j k(W_j,w)$ implies that, for each observation $i=1,\dots,n$, $g(W_i)$ admits the form $g(W_i) = \sum_{j=0}^n \delta_j \big( k(W_j,W_i) - k(W_j,W_0) \big) = [K\delta]_{i+1} - [K\delta]_{1} $. }

where $K$ is the $(n+1)\times (n+1)$ Gram matrix whose elements are given by $k(W_i,W_j)$ for $i,j=0,1,\dots,n$, and $[K\delta]_{j}$ denotes the $j$-th element of the vector $K\delta$.

Dimension reduction by spectral cut-off

The optimization in (ref) is over $\delta \in \mathbb R^{n+1}, \tau\in \mathbb R^{J_n}$. Following the numerical evidence provided by fenton-gallant-96-joe, one may set $J_n=n^{1/5}$. We also find that using cross-validation in practice to select $J_n$ often results in a relatively small choice of $J_n$. When the sample size is large, the optimization over $\tau$ is generally not demanding since $J_n$ remains relatively small. However, the optimization over $\delta \in \mathbb R^{n+1}$ can be computationally intensive. To address this, we reduce the dimensionality by using the leading eigenvectors to approximate the $(n+1)\times (n+1)$ matrix $K$.

More specifically, let $(\hat \lambda_j)_{j=0}^n$ be the eigenvalues of $K$ in descending order and $(\hat u_j)_{j=0}^n$ be the associated orthonormal eigenvectors. Let $\hat U_m$ be the $(n+1)\times m$ matrix collecting the first $m$ columns of $\hat u_{j}$ for $j=0,1\dots,n$, and let $\hat \Lambda_m$ be the $m\times m$ diagonal matrix collecting the first $m$ eigenvalues of $(\hat u_j)_{j=0}^n$. Then, by the Eckart-Young-Mirsky theorem (see, e.g., Theorem 2.4.8 in golub-vanLoan-13), the best rank-$m$ approximation of the matrix $K$, under both the operator norm and the Frobenius norm, is given by \[ K \approx \hat U_m \hat \Lambda_m \hat U_m' . \] In particular, the approximation of $K$ under the operator norm implies that we can approximate $K\delta$ in the Euclidean norm by \[ K\delta \approx \hat U_m \hat \Lambda_m \hat U_m' \delta = \hat U_m \zeta_{\delta}, \quad \text{where} \quad \zeta_{\delta}= \hat \Lambda_m \hat U_m'\delta \in \mathbb R^{m}. \]

Motivated by the approximation described above, i.e., using a low-rank approximation of $K$ based on its first $m$ principal components (PC), a reduced problem of (ref) is given by

where $[\hat U_m \zeta]_j$ denotes the $j$-th element of $\hat U_m\zeta$. Let $(\hat \zeta_{pc}, \hat \tau_{pc})$ be a solution to (ref), and let $\hat \delta_{pc} = \hat U_m \hat \Lambda_m^{-1} \hat \zeta_{pc}$.\footnote{This is by the approximation $\delta = \sum_{j=0}^n \hat u_j \hat u_j' \delta \approx \hat U_m \hat U_m' \delta = \hat U_m \hat \Lambda_m^{-1} \hat \Lambda_m \hat U_m' \delta = \hat U_m \hat \Lambda_m^{-1} \zeta$. } The PC regularized KNP estimator $(\hat g_{pc}, \hat F_{pc})$ is given by $(\hat \delta_{pc}, \hat \tau_{pc})$ via

The optimization in (ref) reduces the dimension of the one in (ref) from $n+J_n$ to $m +J_n$. Here, $m$ is expected to increase with the sample size $n$, both in theory and in practice, but we suppress this dependence for notational simplicity.

\paragraph{Numerical optimization based on analytical gradients:} For implementation, we solve the constrained nonlinear minimization problems in (ref) or (ref) by standard numerical optimization routines in existing software.\footnote{ One can formally write a Lagrangian and consider an unconstrained penalized formulation with a Ridge-type penalty in place of the constraint in (ref) and (ref). However, because the objective is nonconvex, strong duality is not guaranteed, so the constrained and penalized formulations need not be equivalent. }$^{,}$ \footnote{The objectives in (ref) and (ref) are not convex in $(\delta,\tau)$ and $(\zeta,\tau)$, and may have multiple minimizers. As in Remark (ref), $\hat\theta$ is defined as any minimizer of $\hat Q(\theta)$, so any minimizer $(\hat\delta,\hat\tau)$ of (ref) yields an estimator via (ref). Furthermore, any minimizer $(\hat\zeta,\hat\tau)$ of the dimension reduced objective (ref) yields, through (ref), an approximate estimator that shares the same theoretical properties, provided $m$ is chosen appropriately. } In this paper, we use the fmincon routine in MATLAB with the default interior-point algorithm to handle the constraint, and we provide the routine with the analytical gradients to facilitate computation.\footnote{ We use zero as the initial value for the numerical optimization throughout, and the obtained estimates perform well in the Monte Carlo experiments. We also experimented with different starting values, and they gave similar objective values; see Remark (ref) for details. }

In practice, we choose $m$, along with $J_n, B_n$, based on multi-fold cross-validation. In the theoretical framework of this paper, the number $m$ of eigenvectors used to approximate $K$ needs to satisfy certain conditions to ensure that $\hat\theta_{pc}$ is a near minimum of the optimization problem (ref). See Theorems (ref), (ref) and Corollary (ref) for the specific conditions, which depend on the eigenvalue decay of the Gram matrix $K$ and the radius $B_n$.

KNP for CCP, APEs, and conditional APEs

The KNP estimation procedure provides a natural plug-in estimator for CCP, in the form $\hat p(x) = \hat F(v+\hat g(w))$. Moreover, the derivatives of $\hat p(x)$, i.e. $\frac{\partial \hat p(x)}{\partial x} = \hat f(v+\hat g(w)) \frac{\partial (v+\hat g(w) )}{\partial x}$, can be easily evaluated using the estimated density $\hat f$ indexed by $\hat \tau$ and the derivatives of the kernel function. In particular, from the form of $\hat g$ in (ref) or (ref), the derivative $\frac{\partial (v+\hat g(w) )}{\partial x}$ (and hence $\frac{\partial \hat p(x)}{\partial x} $) involves only derivatives of the kernel function. Accordingly, the computation remains essentially unchanged as the dimension $d_w$ of $W$ increases.\footnote{ The computation of partial effects here does not require constructing and differentiating a large set of basis functions whose number can grow rapidly with $d_w$, as in classical series/sieve methods. }

This facilitates the estimation of weighted average partial derivative estimators. In particular, the APE of $X$ and its KNP estimator are given by \[ APE_{x} = \mathbb E \frac{\partial }{ \partial x} p_0(X), \quad \widehat{APE}_x = \frac{1}{n} \sum_{i=1}^n \frac{\partial }{ \partial x} \hat p(X_i). \] The APE of $W$ conditioning on $X\in \mathcal S$ for some region $ \mathcal S\subset \mathcal X$ and its estimator are \[ cAPE_{x| \mathcal S} = \mathbb E \Big( \frac{\partial }{ \partial x} p_0(X) \Big| X\in S\Big), \quad \widehat{cAPE}_{x| \mathcal S} = \frac{ \frac{1}{n} \sum_{i=1}^n \left( 1\{X_i\in \mathcal S\} \frac{\partial }{ \partial x} \hat p(X_i) \right) }{ \frac{1}{n} \sum_{i=1}^n 1\{X_i\in \mathcal S\} }. \]

Practical Implementation

Let $\{(Y_i, V_i,W_i)\}_{i=1}^n$ be the sample. Specify a point $w_\ast \in \mathcal{W}$ for the location normalization $g(w_\ast)=0$. In practice, it is convenient to set $w_\ast =0$ and standardize $(W_i)_{i=1}^n$ to have sample mean zero by subtracting their sample average. Specify a kernel $k$, e.g., $k(s,t)=\exp(-\|s-t\|^2/2)$ used in this paper's Monte Carlo experiments and the empirical application.

Given $J_n, B_n$ and $m$, the PC-regularized KNP estimator is implemented as follows.

MATLAB code for the implementation of the proposed estimator, together with replication packages for the simulation studies and the empirical application, is available online on the author’s website.

Theoretical Properties

This section presents the main theoretical properties of the proposed KNP estimator, both with and without spectral cut-off regularization. In the following subsections, we establish the consistency of the estimator for $\theta_0 = (g_0, F_0)$, the convergence rate of the estimated CCP, and asymptotic normality of weighted average partial derivatives of the CCP.

For the theory of the PC regularized KNP estimator defined through (ref) and (ref), we provide a bound on the difference between the values of $\hat Q$ at $\hat \theta$ and at $\hat \theta_{pc}$, based on which we impose conditions so that $\hat \theta_{pc}$ can be viewed as a near optimum solution when optimizing $\hat Q(\cdot)$ over $\mathcal G_n \times \mathcal F_n$.

Recall that $\hat \lambda_m$ is the $m$-th eigenvalue of the Gram matrix $K$, and $B_n$ is the radius of the RKHS ball used when approximating $g_0$.

Consistency of the KNP estimator

To establish consistency, we need to impose some assumptions. We define the weighted Sobolev norm \[ \|f\|_{m_0+m_e,2,\eta_0} := \left( \sum_{0 \leq \lambda \leq m_0+m_e} \int \left| f^{(\lambda)} (u) \right|^2 (1+u^2)^{\eta_0} du \right)^{1/2}, \] where $m_0, m_e$ are positive integers, constant $\eta_0>1/2$, and we focus on $\eta \in (1/2,\eta_0)$.

Now we are ready to establish consistency of the proposed estimator for $\hat \theta = (\hat g,\hat F)$, as well as the estimator $\hat p $ of the true conditional probability function $p_0$ given by

We denote by $\hat p_{pc}$ when the PC regularized KNP estimator $\hat \theta_{pc}$ is used.

Rate of Convergence

Now we consider the convergence rates for the estimator $\hat p$ of the conditional probability function $p_0(x) = \mathbb{P}\{ Y=1| X=x \}$, under the $L_2(P_{X})$ norm.

We need to impose some technical conditions.

Define the space $L_2(X) := L_2(P_X) := \big\{ h:\mathcal X\to \mathbb R \big| \int h(x)^2 P_X(dx) <\infty \big\}$, where $P_X$ denotes the distribution of $X$, equipped with the norm \[ \|h\|_{L_2(X)} := \bigg( \int h(x)^2 P_X(dx) \bigg)^{1/2} . \] The following theorem provides the convergence rate of $\|\hat p - p_0 \|_{L_2(X)}$. In particular, $\|\hat p - p_0 \|_{L_2(X)} = \big( \int \big( \hat p(x) - p_0(x) \big)^2 P_X(dx) \big)^{1/2} = \big( \mathbb{E} \big( \hat p(X) - p_0(X) \big)^2 \big)^{1/2}$.

As in the sieve literature, the two terms in the rate $\delta_n$ can be interpreted as measures of variance and bias, respectively. Specifically, the first term, $\gamma_n^2$, increases with $B_n$ and $J_n$, reflecting the complexity of the sieve $\Theta_n = \mathcal G_n \times \mathcal F_n$, which arises from the covering numbers of $\mathcal G_n $ and $\mathcal F_n$. This term can be interpreted as a measure of variance. The second term $ \left(\log B_n\right)^{-m_w/2} +J_n^{-m_e}$ decreases with $B_n$ and $J_n$, which arises as the square of the deterministic approximation error when using functions in $\Theta_n$ to approximate $\theta_0 = (g_0, F_0)$. Choosing $B_n, J_n$ to balance these two terms in $\delta_n$ yields the following rate of convergence.

Asymptotic Normality of Weighted Average Derivatives

In this subsection, we establish the asymptotic normality of the weighted average partial derivatives of $\hat p$. The results show that the proposed KNP approach can be used to estimate other functionals of the CCP, including APEs and, when accounting for heterogeneity, conditional APEs, which are often of policy interest.

We consider the weighted average partial effect of the $j$-th element of $X$. For that, we define the functional $\gamma: \Theta \to \mathbb{R}$ given by

where $b(x)$ is a weighting function defined on $ \mathcal{X}:= \text{Supp}(X)$, i.e. $\int b(x) dx = 1$ and $b(x) \geq 0$. Assume here that $b(\cdot)$ is zero outside some compact set. Then integration by parts gives

and $f_X(\cdot)$ denotes the density of $X$. Assume that $f_X$ is bounded away from zero on the set where $b(x)$ is positive. Note that $p(x,\theta) = F(v+g(w))$ is smooth in $\theta$ due to the conditions imposed on the parameter space $\Theta = \mathcal G\times \mathcal{F}$; See Lemma (ref) in Appendix (ref) for its pathwise derivatives. Consequently, $\gamma(\theta)$ is smooth, although nonlinear.

Let $\gamma_0 = \gamma(\theta_0)$ denote the true weighted average derivative, and $\hat \gamma = \gamma(\hat \theta)$ be given by the KNP estimator $\hat \theta$, with or without PC regularization. The following theorem establishes the limit distribution of the estimator $\hat \gamma := \gamma(\hat\theta)$.

The limit distribution in Theorem (ref) is the same as in Theorem 3 in newey-97, which estimates CCP $p_0(x) = \mathbb{E} (Y | X=x)$ using series regression, and obtains the plug-in estimator for the weighted average derivatives. In our approach, we estimate the latent structure, including both the systematic function and the density of the error term, beyond the reduced form CCP. Theorem (ref) shows that, in terms of estimating the weighted average derivatives, the KNP procedure provides an estimator with the same asymptotic variance as in newey-97. However, since $\theta$ enters the objective function in a highly nonlinear manner compared to the cases when approximating a regression function, we need to impose more assumptions than newey-97 in Assumption (ref) to control the high-order terms in the expansions of the objective functions and the functional $\gamma(\theta)$.

Simulation Studies

Compared to parametric and semi-parametric estimation methods, the proposed estimator is expected to perform well in large samples, as it is robust to misspecification of both the systematic function of the covariates and the distribution of the error term. For the usefulness of the proposed estimator, below we examine its finite sample performance by a series of simulations, in order to see (i) if there exists serious issues of efficiency loss relative to a correct fully parametric specification, and (ii) if the proposed estimator is effective in situations where the correct parametric specification is unusual. In addition, we examine whether the bootstrap confidence intervals for APEs and cAPEs have reasonable finite-sample coverage.

We consider the model $Y = 1\{ V + g_0(W) -\varepsilon >0 \}$, where $\varepsilon$ is independent of $X = (V,W')'$. We let $V =_d \mathbb{N}(0,1)$.

Unidimensional $W$

We first focus on unidimensional $W$, and let $W =_d \text{Unif}[-2,2]$. We consider two specifications for $g_0$, where the first one corresponds to the commonly assumed linear index model, and the other is nonlinear.

Note that at $w_\ast = 0$, $g_0(w_\ast)=0$ under (I) or (II). Specification (II) is chosen so that $g_0$ does not lie in any Gaussian RKHS or any finite-order polynomial RKHS.\footnote{The reproducing kernel Hilbert spaces with Gaussian kernels do not contain any nonzero constant, nor any finite-order polynomials. See, e.g., Theorem 2 in minh-10. For the $q$-th order polynomial kernel $k(s,t) = (1+s't)^q$, its RKHS is effectively finitely dimensional, with basis functions consisting of all polynomials up to order $q$. } The error term $\varepsilon$ is given by one of the two following specifications.

where B indicates the mixture of two normal distributions $\mathbb N(-3,1)$ and $\mathbb N(2,1)$, i.e. with probability 1/4, $\varepsilon$ follows $\mathbb N(-3,1)$, and with probability 3/4, $\varepsilon$ follows $\mathbb N(2,1)$.

For simplicity, we refer to the specifications introduced above as (I) and (II) for the systematic function and (A) and (B) for the error distribution. Moreover, we will refer to as (IA), (IB), (IIA), and (IIB) the four cases that are given by the combinations of I and II with A and B.

We compare the KNP estimator with (a) Kernelized probit (KPB), which specifies the standard normal error term and approximates $g_0$ based on functions in RKHS as in KNP (b) Semi-Nonparametric (SNP), which approximates $F_0$ using gallant-nychka-87's method and specifies $g_0$ as a linear function, (c) Probit, which specifies standard normal $\varepsilon$ and linear function of $g_0$. In addition, we consider methods specifying standard normal $\varepsilon$ and approximating $g_0$ based on 2nd, 3rd, 4th polynomials, respectively. For RKHS-based methods, we use the Gaussian kernel $k(s,t) = \exp(-\|s-t\|^2/2 )$. For both KNP and KPB, the number $m$ of eigenvectors in the spectral cut-off regularization are selected based on 5-fold cross-validation.

Table (ref) compares the performance of these methods for estimating $g_0, p_0$ under each of the four specifications (IIB), (IIA), (IB), (IA), , based on Monte Carlo simulations with $Nsim = 1000$ replications and sample size $ntrain=2000$ observations. The table reports $\text{RMSE}(\hat g) = \sqrt{\mathbb E(\hat g(W)- g_0(W) )^2}, \text{MAD}(\hat g) = \mathbb E |\hat g(W)-g_0(W) |$, along with $\text{RMSE}(\hat p) = \sqrt{\mathbb E(\hat p(X)- p_0(X) )^2}, \text{MAD}(\hat p) = \mathbb E |\hat p(X)-p_0(X) |$, where the expectations are estimated using sample means of $ntest=10,000$ observations in test sample. More details of the simulation procedure are given in the footnote of Table (ref).

Table (ref) shows that under specification (IIB), KNP provides the best estimators for $g_0$ and $p_0$, which are much better than all of the other methods. This suggests that the proposed estimator is effective in situations where the correct parametric specification is unusual. Under specification (IIA), KPB performs best as expected, followed closely by KNP, whereas all of the other methods are much worse. Under specification (IB), SNP performs best as expected, and KNP again does the second best. In particular, KNP's estimates for $p_0$ are very close to that of SNP. In specification (IA), where the probit model is correctly specified, the KNP estimator performs comparably to the probit model. This suggests that the efficiency loss of using the proposed method relative to a correct fully parametric specification is rather mild.

While Table (ref) uses sample size $2000$, Tables (ref), (ref) report the comparisons using sample sizes $1000, 500$, respectively. The comparisons using the sample size $1000$ are the same as above. When using the sample size $500$ in Table (ref), one difference is that under (IIB), KNP's estimator for $g_0$ does slightly worse than that of KPB. However, in terms of estimating $p_0$, KNP's estimator is still the best and much better than all of the other methods.

As supplements to Table (ref), Fig. (ref) presents $Nsim=1000$ simulated estimates of $g_0$ given by each of the methods using sample size $ntrain=2000$ under specification (IA) and (IIB). We note that, compared to other methods, KNP fits the true function best when $g_0$ is nonlinear under (IIB) and also performs well when $g_0$ is linear under (IA).

\afterpage{

}

\afterpage{

}

10-Dimensional $W$

We now consider a case where $W$ is 10-dimensional to evaluate the performance of the proposed estimator in more complex settings. We let $W = (W_{1},\dots,W_{10})'$, where each $W_{j} =_d \text{Unif}[0,1]$ for $j=1,\cdots,10$ and is independent of each other. We consider two specifications for $g_0$, where the first one corresponds to the commonly assumed linear index model, and the other is nonlinear.

where $\beta = (0.63, 0.81, -0.75, 0.83, 0.26, -0.80, -0.44, 0.09, 0.92, 0.93)'$ \footnote{These numbers are generated as the first 10 numbers from $\text{Unif}[-1,1]$ using rng(`default') in Matlab. } The error term $\varepsilon$ is independent of $X$ and is given by specifications (A) and (B) as before. We will refer to as (IIIA), (IIIB), (IVA), and (IVB) the four cases that are given by the combinations of III and IV with A and B.

Similar as Table (ref), Table (ref) presents the simulation results for designs (IVB), (IVA), (IIIB), (IIIA) using $Nsim = 1000$ replications, sample size $ntrain \in\{2000,5000,10000\}$ for estimation and $ntest =1,000,000$ for out-of-sample evaluation. Table (ref) further demonstrates that the KNP estimator is robust to both misspecification of the systematic function of the covariates and misspecification of the density of error term. Moreover, for moderate sample sizes, the KNP estimator shows desirable properties: It effectively improves the finite sample performance in case of misspecification, and has a rather mild efficiency loss if the model is correctly specified.

One unexpected pattern in Table (ref) is that, in designs (IIIA) and (IIIB), $\text{RMSE}(\hat g)$ and $\text{MAD}(\hat g)$ are not monotone in $ntrain$, even though $\text{RMSE}(\hat p)$ and $\text{MAD}(\hat p)$ decrease. This pattern arises from the data-driven tuning, where $J_n,m, B_n$ are selected by 5-fold cross-validation to minimize the squared loss in the objective (ref) for $p(\cdot)$, which is a composition of $g$ and $F$ and thus targets the error of $\hat p$ rather than the error of $\hat g$. To check that this pattern is indeed driven by the data-driven model selection, we conduct an additional experiment in which $J_n$ is fixed at $J_n=1$ for design (IIIA) and $J_n = 3$ for design (IIIB), while $m$ is still chosen by cross-validation as before and $B_n=100$ is fixed. Table (ref) in Appendix (ref) reports the corresponding results. In this setting, $\text{RMSE}(\hat g)$ and $\text{MAD}(\hat g)$ decrease with the training sample size, and $\text{RMSE}(\hat p)$ and $\text{MAD}(\hat p)$ also decrease, which is consistent with the theoretical properties of the KNP estimator.

Coverage of Bootstrap Confidence Intervals for APE and cAPEs

Theorem (ref) establishes the asymptotic normality of weighted average partial derivatives. In the empirical application in Section (ref), we construct nonparametric bootstrap confidence intervals for inference, since the asymptotic variance in Theorem (ref) depends on the density of $X$ and its derivatives, which are challenging to estimate given that $X$ is 10-dimensional in that example.\footnote{While it may be possible to adopt methods such as kernel density estimation or the ones proposed in spady-stouli-20 to estimate the asymptotic variance, we opted for bootstrap intervals for practicality. Moreover, Theorem (ref) shows a standard root-$n$ asymptotically normal distribution centered at zero, so arguments as in chen-linton-vankeilegom-03 can, in principle, be adapted to show the bootstrap validity under suitable conditions. We leave a formal proof of bootstrap validity for future research. } The empirical application focuses on the APE of one component of $W$ and on two conditional APEs (cAPEs) conditioning on this component being below or above its average.

\afterpage{

}

Motivated by this empirical application, we examine the finite-sample coverage of bootstrap confidence intervals (CIs) in the 10-dimensional $W$ design (IVB). We focus on the APE of $W_1$ and two cAPEs conditioning on $W_1<\mathbb E W_1 =0.5$ and $W_1>0.5$, denoted by $APE_{W_1}, cAPE_{W_1|W_1<0.5}, cAPE_{W_1|W_1>0.5}$, respectively. We generate $Nsim=1000$ Monte Carlo replications, with sample size $ntrain\in\{2000, 5000, 10000\}$. For each replication, we generate a sample of size $ntrain$, fit the KNP estimator, compute the estimates for the APE and cAPEs, and draw $Nboot = 1000$ bootstrap samples and re-estimate APEs to construct 90%, 95%, and 99% bootstrap confidence intervals. Bootstrap confidence intervals are constructed using the basic (reverse-percentile) bootstrap method.

Table (ref) reports the empirical coverage probabilities. The results indicate that the bootstrap CIs have coverage generally close to the nominal levels for both the APE and the two cAPEs: although the CIs are slightly conservative at $ntrain=2000$, coverage tends to be closer to the nominal as the sample size increases. Table (ref) in Appendix (ref) also reports the average CIs length, which decreases with sample size.

Application: Temperature and Judge’s Decision

heyes-saberian-19, heyes-saberian-22 analyze the effect of outdoor temperature on the probability of an asylum application being granted. Based on a linear probability model (LPM) including other weather and pollution characteristics, heyes-saberian-19 finds that, in their preferred specification, a $10^{\circ}$F increase in case-day temperature reduces the grant probability by 1.075 percent. Their results suggest that high temperatures may damage decision consistency, even for experienced professional decision-makers who work indoors and are “protected" by climate control. The evidence that such socially and economically important high-stakes decisions can be affected by extraneous variables suggests a potential welfare loss.

How do temperature and other weather and pollution characteristics affect a judge’s grant decision? heyes-saberian-19 highlight that their findings are consistent with established links from temperature to mental function, decision-making, risk attitudes, and mood. Following their argument, we may naturally expect these environmental variables to affect judges’ perceived utility from the outdoor environment, thereby influencing their grant decisions.

Two further questions arise. First, is the effect of temperature on judges' perceived utility approximately linear, or does it vary nonlinearly with temperature? Second, do environmental variables (e.g., temperature, air pressure, dew point, PM2.5) affect utility in an additive and separable way, or do they interact with each other when affecting decision-makers' mental states? A fully nonparametric utility function of environmental variables is useful for addressing both questions.

In this section, we apply KNP to investigate the effect of temperature on judges' grant decisions, allowing for these important features. This exercise illustrates the empirical relevance of our proposed method. We find a bell-shaped relationship between temperature and judges' perceived utility, which is consistent with existing evidence that outdoor environments may affect decisions through decision-makers' mental states or moods. Moreover, the effects of temperature on judges' decisions are heterogeneous across temperature ranges: while the conditional APE at high temperatures is significantly negative, the conditional APE at low temperatures is not statistically significant at the 10% significance level. This heterogeneity can potentially inform more targeted policy interventions across temperature ranges.

Empirical Specification and Implementation

Let $Y_i= 1$ if the application case $i$ is granted, and 0 otherwise. We apply the proposed KNP estimation procedure to the following model:

where $Y_i^\ast$ is a latent index that can be viewed as an unobserved score of case $i$. For each case $i$, let $j$ denote its assigned judge, $t$ the decision time (year-month-day), and $c$ the applicant's country of nationality. $W_i$ is the vector of 9 outdoor environmental variables that case $i$'s assigned judge $j$ was exposed to at the decision time $t$, including mean daily temperature, air pressure, dew point, precipitation, wind speed, sky cover, ozone, CO, PM2.5; see Section II of heyes-saberian-19 for the definitions and description of $W_i$ and $Y_i$, and Table 1 therein for summary statistics. Note that $W_i$ depends only on the judge and the decision time; hence, we also write $W_i = W_{jt}$. Here $g(W_i)$, or $g(W_{jt})$ more specifically, may be interpreted as the utility given by the outdoor environment with variables $W_{jt}$.

Variable $V_i$ is chosen to be the log-odds of the mean approval rate for different types of applications from country $c$ over each month.\footnote{ As explained in heyes-saberian-19, “There are two types of cases in immigration courts: affirmative cases in which the applicant presents in the courts on her/his own and defensive cases in which the applicant is instructed to attend on the initiative of the immigration authorities."} Since the observed case characteristics in the data are limited, this choice allows us to account for some of the heterogeneity in the applicant's nationality, types of application, and time at the level of year-month. Fig. (ref) in Appendix (ref) plots the histogram and a kernel density estimate of $V_i$ based on 99{,}773 observations, using Silverman’s rule-of-thumb bandwidth. The empirical range of $V_i$ is $[-3.761, 2.772]$; the $0.1\%, 0.5\%, 1\%, 50\%$, $ 99\%, 99.5\%$ and $99.9\%$ quantiles are $-3.71, -3.32, -3.08, -0.45$, $1.61, 1.87$, and $2.20$, respectively. Fig. (ref) and the quantiles indicate that $V_i$ takes many distinct values and has a large variation in the data, supporting its use as a regressor with a large support for identification.

The error term $\varepsilon_{i}$ is the idiosyncratic case-level component of the latent score $Y_i^\ast$, after accounting for the baseline approval propensity $V_{i}$ and the environmental influence $g_0(W_{jt})$. We assume that $\varepsilon_{i}$ is independent of $V_i,W_i$, where $\varepsilon_i$ has an unknown CDF $F_0$. This condition is crucial for identification in our model and is treated as a maintained assumption in this illustration. By construction, $V_{i}$ is a proxy for baseline approval propensity by applicant nationality, calendar month (year-month), and case type. It is intended to absorb systematic heterogeneity along these dimensions, including at the $(c,t)$ level, so that the remaining variation in $\varepsilon_i$ reflects the residual within-cell idiosyncratic component. In addition, environmental exposure $W_{jt}$ is plausibly exogenous in this setting because, as stated in heyes-saberian-19, “the setting of dates for cases and the rostering of judges is done well in advance”, and weather on the decision date is not manipulable by judges or applicants.

For implementation, we estimate $(g_0, F_0)$ by the PC-regularized KNP method; that is, we solve (ref), after which the estimator is given in (ref), and then the APE and cAPEs are computed following Section (ref). Specifically, we standardize each component of $W_{jt}$ to have zero mean and unit variance and set $w_\ast = 0$, so the location normalization $g_0(w_\ast) = 0$ normalizes the utility to zero at the mean environmental variables (in the original units). We use the Gaussian kernel $k(s,t)=\exp(-\|s-t\|^2/2)$. The tuning parameters $B_n, J_n, m$ are selected by 5-fold cross-validation, which yields the choice $B_n=\sqrt{10}, J_n=4, m = 50$.\footnote{ The candidate sets are $B_n^2 \in \{5,10,20,\dots, 50, 75, 100, 200, 500\}$, $J_n\in\{0,1,\dots, 7\}$, and $m\in \{10,25, 50, 75,100, 150,200\}$. For each triplet $(B_n,J_n,m)$, we solve (ref) on the training folds and evaluate the squared-loss objective on the held-out fold. We choose the triplet that minimizes the average validation loss across folds. } We solve the constrained optimization in (ref) using the MATLAB fmincon routine with the default interior-point method to handle the constraint. We supply the routine with analytical gradients with respect to $\zeta$ and $\tau$; the closed-form expressions for the objective, constraint, and gradients are given in Appendix (ref). All numerical optimizations are initialized at zeros. When constructing bootstrap confidence intervals, we use a nonparametric bootstrap with $1000$ bootstrap replications, resampling cases $i$ with replacement. In each bootstrap replication, we re-estimate $(g_0,F_0)$ using the same procedure as above while fixing $B_n=\sqrt{10}, J_n=4, m = 50$ at the values selected by cross-validation in the original sample. Bootstrap confidence intervals are constructed using the basic (reverse-percentile) bootstrap method.\footnote{ Specifically, for any scalar target parameter $\gamma_0$, let $\hat\gamma$ be the estimate obtained from the original sample, and let $q^\ast_u$ denote the $u$-percentile of the centered bootstrap distribution $\{\hat\gamma^\ast_b - \hat\gamma\}_{b=1}^{1000}$, where $\hat\gamma^\ast_b$ denotes the estimate of $\gamma_0$ obtained from the $b$-th bootstrap replication. The $(1-\alpha)$ confidence interval for $\gamma_0$ is given by $[\hat \gamma - q^\ast_{1-\alpha/2}, \hat \gamma - q^\ast_{\alpha/2}]$. In our application, $\gamma_0$ may represent a pointwise value $g_0(w)$ at any prespecified evaluation point $w$, an APE or a cAPE. }

Results

Our empirical goals are twofold: (a) to characterize how temperature enters the utility function $g_0(\cdot)$ of outdoor environmental variables, and (b) to quantify how temperature affects judges' grant decisions.

For (a), Fig. (ref) presents the estimated utility as a function of temperature in the original units, with all other environmental variables fixed at the mean level. The shaded area indicates 90% pointwise bootstrap confidence band. The figure shows that, when the temperature is at a high level, increasing temperature decreases utility. In contrast, when the temperature is at a low level, increasing temperature will increase the utility.

For (b), we estimate the average partial effect (APE) and conditional average partial effects (cAPEs) of temperature $T$. In particular, other than the APE \[ APE_T = \mathbb{E} \frac{\partial}{\partial T} p(X), \] we consider the cAPEs conditioning on temperature higher (or lower) than $70^{\circ}$F, \[ cAPE_{T|T>70^{\circ}F} = \mathbb{E} \left( \frac{\partial}{\partial T} p(X) \bigg| T>70^{\circ}F \right), \quad cAPE_{T|T<70^{\circ}F} = \mathbb{E} \left( \frac{\partial}{\partial T} p(X) \bigg| T<70^{\circ}F \right). \] We choose $70^\circ$F as a benchmark close to typical indoor comfort temperature (“room temperature”) and use it to summarize heterogeneity in the temperature effect across low- and high-temperature ranges.

Table (ref) reports 90% bootstrap confidence intervals from KNP for the APE and cAPEs, based on 1000 replications. We compare the KNP results with those from an LPM benchmark that regresses $Y_i$ on the same covariates $(V_i,W_i)$ as in (ref). Under the LPM, the estimated APE and both cAPEs are significantly negative, including $cAPE_{T|T<70^{\circ}F}$. In contrast, KNP yields an insignificant $cAPE_{T|T<70^{\circ}F}$ and a significantly negative $cAPE_{T|T>70^{\circ}F}$, with the overall APE also significantly negative.

These results complement heyes-saberian-19,heyes-saberian-22 by providing two empirically relevant outputs: an estimate of $g_0$ as a fully nonparametric function of the environmental variables $W$, and cAPEs of temperature across temperature ranges. Notably, the estimated heterogeneous effects of temperature can potentially inform more targeted policy interventions. While an LPM can allow for nonlinearity by adding higher-order and interaction terms in environmental variables, this quickly yields many polynomial terms (e.g., a 4-th order polynomial expansion of the 9-dimensional $W$ yielding 714 terms) and hard-to-interpret specifications. In contrast, the estimation of cAPEs in our method is straightforward, with a computational procedure that remains essentially unchanged regardless of the dimension of $W$, since it does not require constructing or differentiating a large interaction basis.

In terms of the average effect, our estimated APE is significantly negative and qualitatively consistent with heyes-saberian-19, although the magnitude is smaller.\footnote{ heyes-saberian-19 reports a 1.075% decrease in grant probability per $10^{\circ}$F increase in temperature, while our 90% CI implies a change in $[-0.378\%, -0.065\%]$ per $16.9^{\circ}$F (1 standard deviation) increase. } However, because we use a nonparametric nonlinear model, whereas heyes-saberian-19 uses an LPM with a rich set of fixed effects, differences in magnitude should not be taken literally, as they may reflect the different specifications. While our approach has the advantage of capturing flexible nonlinear and interactive effects of the full environmental vector within a latent utility framework,\footnote{Under the maintained assumptions, the latent-utility formulation provides a structural mapping from $(V,W)$ to the choice probabilities, which can be used to analyze counterfactuals holding the distribution of the error term fixed. } the tradeoff is that it does not accommodate a very fine set of fixed effects in the same way as linear models; however, coarse controls (e.g., year indicators) can be incorporated by including dummies in the systematic component.

Conclusion

In this paper, we propose a new estimation procedure for a class of identified nonparametric binary choice models. Compared to other possible methods, our estimation procedure is amenable to easier computation, especially when the number of covariates is non-small which may lead to a large number of basis functions if using the commonly used sieve method.

We show the proposed estimator has desirable asymptotic properties, and simulation studies suggest that the KNP estimator works well in finite samples. We demonstrate the practical relevance of the proposed method by revisiting the effect of temperature on immigration judges’ latent utility, and thus, their decisions on asylum applications.

In future work, a natural extension is to allow for further structures, motivated by either economic application-specific assumptions or existing theory, on the nonparametric and distribution-free BCM, without losing the tractability of the KNP procedure in terms of computation.\footnote{See, e.g., chetverikov-santos-shaikh-18 and the reference therein for a review of the roles of shape restrictions in identification, estimation and inference.} For example, one may expect some covariates to affect the latent utility partially linearly. Such a partial linear structure of latent utility can be easily incorporated into the KNP estimation procedure. However, for more complex structures, such as monotonicity or concavity/convexity of the latent utility function in some covariates, the computation will be much more difficult under big data environments, since it involves more constraints on the function evaluations at the data points to restrict the shape of the function.

Appendices