Estimating High Dimensional Monotone Index Models by Iterative Convex Optimization1

Estimating High Dimensional Monotone Index Models by Iterative Convex Optimization

\vskip -.5in { \singlespace

}

{

}

\

{\bf Key Words} Monotone Index models, Convex Optimization, Kernel and Sieve Estimation.

\thispagestyle{empty}

\setcounter{page}{1}

{5pt} {5pt} {5pt} {5pt}

\setcounter{equation}{0}

Introduction

Monotone index models have received a great deal of attention in both the theoretical and applied econometrics literature, as many economic variables of interest are of a limited or qualitative nature. A leading special case in this class is the binary choice model which is usually represented by some variation of the following equation:

where $I[\cdot]$ is the usual indicator function, $y_i$ is the observed response variable, taking the values 0 or 1 and $x_i$ is an observed $p$ dimensional vector of covariates which effect the behavior of $y_i$. Both the scalar disturbance term $u_i$ with distribution function denoted by $G(\cdot)$, and the $p-$ dimensional vector $\beta_e^\ast$ are unobserved, the latter often being the parameter estimated from a random sample $(y_i,x_i') \ \ i=1,2,...n$.

The disturbance term $u_i$ is restricted in ways that ensure identification of $\beta_e^\ast$. Parametric restrictions specify the distribution of $u_i$ up to a finite dimensional parameter and assume that $u_i$ distributed independently of the covariates $x_i$. Under such a restriction, $\beta_e^\ast$ can be estimated (up to scale) using maximum likelihood or nonlinear least squares. Estimators that are robust to these parametric distributional assumptions have been proposed and analyzed resulting in a variety of estimation procedures for $\beta_e^\ast$.

An important class of semiparametric restrictions used in the literature were based on independence/index restrictions. Estimation procedures under this restriction include those proposed by han1987non, ichimura1993semiparametric, klein1993efficient. These cover but are not limited to the above binary response model. This class of index models have a robustness advantage over parametric approaches, but estimators within this class are difficult to compute\footnote{Other estimation of index models includes \citep*{stoker1986consistent,powell1989semiparametric}. While these are relatively easy to compute, such derivative based estimators cannot be applied unless all components of $x_i$ are continuously distributed.} due to nonconvexity and in some cases also nonsmoothness of their respective objective functions. Furthermore the difficulty increases with the dimension of $x_i$. Recent work which is motivated by computational concerns is \citet*{ahn2018simple}. However, their two step procedure involves a fully nonparametric estimator in the first stage, so is also not suitable for models with a large number of regressors.

A related drawback of all these procedures is that they are designed to estimate parameters in models of a small and {\em fixed} dimension. A relatively recent and thriving literature in econometrics and machine learning is recognizing the many advantages of allowing for large dimensional models or models with a large set of controls. This class is a special case of models that consider the situation when the dimension of $x_i$ is large, and this is now often modeled with its dimension increasing with the sample size. Due primarily to its empirical relevance, there has been a burgeoning literature on estimation and inference in certain econometric and statistics models with a large number of regressors or a large number of moment conditions. For a surevey of examples in economics and finance, see fanlvqi2011. Recent papers include neweywind2009, chernozhukov2017central,bellonietal2018, cattaneoetal2018a, cattaneoetal2018b,

Related to our work is the recent literature on estimating large dimensional binary choice or monotone index models in sur2019modern and fan2020rank. sur2019modern considers inference in a large dimensional logit model, relying on the logistic distribution of the disturbance term where it is shown that $\chi^2$ asymptotic approximations of the LR statistic are suspect when the dimension of $x$ is large. \citet*{fan2020rank} on the other hand estimate parameters by optimizing the objective function introduced in han1987non, but with the number parameters increasing with the sample size. Optimizing these rank based objective functions is unfortunately hard even with recent developments in algorithms and search methods for optimizing non smooth and/or non convex objective functions. See for example important recent work based on mixed integer programming (MIP) as in, e.g. fan2020rank and ShinTodorov2020.

Therefore, in light of the drawbacks in the existing literature, this paper proposes a new estimation procedure that is amenable to easier computattion. Specifically we aim to construct a computationally feasible estimator for a semiparametric binary choice and monotone index models with {\em increasing} dimension based on a convex objective function and then establish its asymptotic properties. As we will discuss in detail in the next section, our algorithm uses an iterative estimator based on a batched gradient descent (BGD) method, and we show how to use nonparametric methods to approximate the distribution in each stage of the iteration. One is the method of sieves\footnote{ See chen2007large who pioneered the use of sieve methods in econometrics.}, and the other is kernel regression.

\ \ \

\baselineskip = .75\baselineskip {\bf Notation:}

Throughout the rest of this paper, to facilitate the description and properties of estimation procedures we will be using the following notation. For any real sequences $\left\{ a_{n}\right\} _{n=1}^{\infty}$ and $\left\{ b_{n}\right\} _{n=1}^{\infty}$, we write $a_{n}=o\left(b_{n}\right)$ if $\lim\sup_{n\rightarrow\infty}\left|a_{n}/b_{n}\right|=0$, $a_{n}=O\left(b_{n}\right)$ if $\lim\sup_{n\rightarrow\infty}\left|a_{n}/b_{n}\right|<\infty$, and $a_{n}\sim b_{n}$ if both $a_{n}=O\left(b_{n}\right)$ and $b_{n}=O\left(a_{n}\right)$. For any random sequences $\left\{ a_{n}\right\} _{n=1}^{\infty}$ and $\left\{ b_{n}\right\} _{n=1}^{\infty}$, we write $a_{n}=O_{p}\left(b_{n}\right)$ if for any $0<\tau<1$ there are $N$ and $C>0$ such that $P\left\{ \left|a_{n}/b_{n}\right|>C\right\} <\tau$ holds for all $n\geq N$, we write $a_{n}=o_{p}\left(b_{n}\right)$ if for any $C>0$, $\lim_{n\rightarrow\infty}P\left\{ \left|a_{n}/b_{n}\right|>C\right\} \rightarrow0$. For any Borel sets $A\subseteq\mathbb{R}^{k}$, denote its Lebesgue measure as $m\left(A\right)$. For any symmetric matrix $A$, we write $A\succ0$ if $A$ is positive definite, and $A\succeq0$ if $A$ is positive semi-definite. For any symmetric matrices $A$ and $B$, we write $A\succ B$ if $A-B\succ0$ and $A\succeq B$ if $A-B\succeq0$. For any matrix $A$, we denote $\sigma\left(A\right)$ as its singular value, and denote $\overline{\sigma}\left(A\right)$ and $\underline{\sigma}\left(A\right)$ as its largest and smallest singular value. For any symmetric matrix $A$, we denote $\lambda\left(A\right)$ as its eigenvalue, and denote $\overline{\lambda}\left(A\right)$ and $\underline{\lambda}\left(A\right)$ as its largest and smallest eigenvalue. For any vector $\boldsymbol{x}=\left(x_{1},\cdots,x_{p}\right)^{\mathrm{T}}$, we denote its Euclidean norm as $\left\Vert \boldsymbol{x}\right\Vert =\sqrt{\sum_{i=1}^{p}x_{i}^{2}}$. For any matrices $A=\left(a_{ij}\right)_{n\times m}$, we denote $\left\Vert A\right\Vert =\sqrt{\sum_{i=1}^{n}\sum_{j=1}^{m}a_{ij}^{2}}$. Note that when $A$ is positive semi-definite, there holds $\left\Vert A\boldsymbol{x}\right\Vert \leq\overline{\lambda}\left(A\right)\cdot\left\Vert \boldsymbol{x}\right\Vert $; for general square matrix $A$, there holds $\left\Vert A\boldsymbol{x}\right\Vert \leq\overline{\sigma}\left(A\right)\cdot\left\Vert \boldsymbol{x}\right\Vert $. Finally, for any function $f\left(\boldsymbol{x}\right)$ with domain $D$, define $\left\Vert f\right\Vert _{\infty}=\sup_{\boldsymbol{x}\in D}f\left(\boldsymbol{x}\right)$.

\baselineskip = 1.333\baselineskip

The BGD Estimator

To provide some intuition for our semiparametric estimators that will be introduced in the following sections, in this section we consider a simplified version of the model where the cumulative distribution function $G\left(\cdot\right)$ is completely known. Under such setup, we explore the batch gradient descent estimator (BGD estimator) of $\boldsymbol{\beta}_{e}^{\star}$ when its dimensionality $p$ may increase, which is also important on its own right. Throughout the following analysis we assume that the data set satisfies the following assumption.

Given any loss function $\ell_{G}\left(\boldsymbol{\beta}_{e},\mathbf{X}_{e},y\right)$ that depends on $G$ and is a.s. differentiable with respect to $\boldsymbol{\beta}_{e}\in\mathcal{B}_{e}$, the BGD estimator of $\boldsymbol{\beta}_{e}^{\star}$ is based on the following iteration,

where $\delta_{k}>0$ is the learning rate. Note that $n^{-1}\sum_{i=1}^{n}\partial\ell_{G}\left(\boldsymbol{\beta}_{e},\mathbf{X}_{e,i},y_{i}\right)/\partial\boldsymbol{\beta}_{e}$ constitutes a sample analogue of the derivative $\partial\mathbb{E}\left[\ell_{G}\left(\boldsymbol{\beta}_{e},\mathbf{X}_{e},y\right)\right]/\partial\boldsymbol{\beta}_{e}$. Unlike the stochastic gradient descent (SGD) algorithm, in the BGD algorithm, in each round of update we evaluate the derivative of the loss function over all data points. This increases the computational burden but provides a more accurate estimator for the derivative of the expected loss function. Given the initial guess of the parameter, $\boldsymbol{\beta}_{e,1}$, we iterate based on ((ref)) until some terminating conditions are reached.

In this paper, we consider the following loss function

for some sufficiently large positive constant $A$. The loss function ((ref)) was also considered in agarwal2014least and has many nice properties. For instance, under some mild conditions, we can show that

and \[ \frac{\partial^{2}\mathbb{E}\left(\ell_{G}\left(\boldsymbol{\beta}_{e},\mathbf{X}_{e},y\right)\right)}{\partial\boldsymbol{\beta}_{e}\partial\boldsymbol{\beta}_{e}^{\mathrm{T}}}=\mathbb{E}\left\{ G^{\prime}\left(\boldsymbol{\mathbf{X}}_{e}^{\mathrm{T}}\boldsymbol{\beta}_{e}\right)\mathbf{X}_{e}\mathbf{X}_{e}^{\mathrm{T}}\right\} \succ0,\forall\boldsymbol{\beta}_{e}\in\mathcal{B}_{e}. \] So $\boldsymbol{\beta}_{e}^{\star}$ uniquely minimizes $\mathbb{E}\ell_{G}\left(\boldsymbol{\beta}_{e},\mathbf{X}_{e},y\right)$ over $\mathcal{B}_{e}$. Another desirable property of the loss function ((ref)) is that the derivative of ((ref)) with respect to $\boldsymbol{\beta}_{e}$, which is $\left(G\left(\mathbf{X}_{e}^{\mathrm{T}}\boldsymbol{\beta}_{e}\right)-y\right)\mathbf{X}_{e}$, depends only on $G\left(\cdot\right)$ instead of on its derivatives. So when we conduct a semiparametric iteration in the following sections, we only need to nonparametrically approximate $G\left(\cdot\right)$, which is generally more robust compared with approximating its derivatives. Based on loss function ((ref)), the BGD estimator is obtained based on the following iteration

\ \

We now describe the asymptotic properties of $\boldsymbol{\beta}_{e,k}$. We first make the following assumption.

For any $\boldsymbol{\beta}_{e}\in\mathcal{B}_{e}$, define $\Delta\boldsymbol{\beta}_{e}=\boldsymbol{\beta}_{e}-\boldsymbol{\beta}_{e}^{\star}$. Also define $\varepsilon_{i}=y_{i}-G\left(\mathbf{X}_{e,i}^{\mathrm{T}}\boldsymbol{\beta}_{e}^{\star}\right)$, where $\mathbb{E}\left[\left.\varepsilon_{i}\right|\mathbf{X}_{e,i}\right]=0$. When (ref) and (ref) hold, we have the following result.

When $p$ is fixed, (ref)(i) implies that $ \sup_{k\geq k_{1,n}^{BGD}+1}\left\Vert \Delta\boldsymbol{\beta}_{e,k}\right\Vert =O_{p}\left(1/\sqrt{n}\right), $ and (ref)(ii) implies that for $k$ sufficiently large, the BGD estimator is an asymptotically linear estimator, so there holds $ \sqrt{n}\Delta\boldsymbol{\beta}_{e,k+k_{1,n}^{BGD}}\rightarrow_{d}N\left(0,\Sigma_{1}^{\star}\right) $ by the central limit theorem. The asymptotic variance can be estimated based on (ref)(iii). The number of iterations required to obtain root-$n$ consistency, $k_{1,n}^{BGD}$, is determined by many factors including the sample size $n$, the distance between the true parameter and the initial guess $||\Delta\boldsymbol{\beta}_{e,1}||$, as well as the lower bound of the eigenvalues of $M_{n}\left(\boldsymbol{\beta}_{e}\right)$. In general, $k_{1,n}^{BGD}$ is of order $O\left(\log n\right)$, but in practice when we apply the above algorithm, the specific number of iteration is difficult to determine. For detailed discussion of the number of iterations, see (ref) at the end of Section (ref). The inference on $\boldsymbol{\beta}_{e}^{\star}$ based on the BGD estimator is given by (ref)(iv). Note that for any given vector $\rho$, we require that $\frac{1}{\sqrt{n}}\rho^{\mathrm{T}}M^{-1}\left(\boldsymbol{\beta}_{e}^{\star}\right)\sum_{i=1}^{n}\varepsilon_{i}\mathbf{X}_{e,i}$ is asymptotically normally distributed. An alternative approach is to apply the high-dimensional central limit theorem to $\frac{1}{n}\sum_{i=1}^{n}M^{-1}\left(\boldsymbol{\beta}_{e}^{\star}\right) \mathbf{X}_{e,i}\varepsilon_{i}$ chernozhukov2017central.

Before we conclude this section and move to semiparametric estimation, we further comment on (ref). Different from the stochastic gradient descent algorithm toulis2017asymptotic, we show in (ref) that the learning rate $\delta_{k}$ can be selected as a sufficiently small constant. Indeed, in the following results, we show that $\delta_{k}$ can decay to zero at any rate as long as $\sum_{k=1}^{\infty}\delta_{k}=\infty$ holds, and the choice of $\delta_{k}$ will not change the asymptotic results displayed in (ref). In particular, we have the following proposition.

(ref) shows that the choice of the learning rate basically does not affect the convergence rate as well as the asymptotic distribution of the BGD estimators. The main advantage of using a sequence of decaying learning rates is that we do not need to choose the constant $\delta$ as required in (ref), since for $k$ sufficiently large, $\delta_{k}\leq2/\left(3\overline{\lambda}_{e}\right)$ will automatically hold. However, the disadvantage of using decaying learning rates is that such procedure takes much longer time to converge because the magnitude of the update in the $k$-th round decreases as $k$ increases. For instance, suppose that we choose $\delta_{k}\sim k^{-\upsilon}$ for some $0\leq\upsilon<1$, we have that $\sum_{j=1}^{k}\delta_{j}\sim k^{1-\upsilon}$. Then to ensure that $\sum_{j=1}^{\widetilde{k}_{1,n}^{BGD}}\delta_{j}\geq\underline{\lambda}_{e}^{-1}\left(\log n+2\log\left\Vert \Delta\boldsymbol{\beta}_{e,1}\right\Vert \right)$, we need $\widetilde{k}_{1,n}^{BGD}\sim\left(\log n\right)^{\frac{1}{1-\upsilon}}$. Obviously, setting $\upsilon=0$ leads to $k\sim\log n$, which corresponds to the requirement in (ref)(i); when $\upsilon>0$, we can see that more rounds of iteration is needed compared with required in (ref)(i).

Semiparametric BGD Estimation

In the previous section, we focused on iterative estimators based on the BGD algorithm for the parametric binary choice models. We show that when the CDF of the error term is known, the iterative estimators based on the BGD algorithm are consistent and attain asymptotic normality under mild conditions. However, having prior knowledge of the form of $G$ is generally too strong an assumption. In most applications, the source of the individual shock $u$ in (ref) is difficult to justify, which makes it quite difficult, if not completely impossible, to know the exact expression of $G$. In this scenario, the algorithm proposed in the previous section is infeasible. To overcome such problem, this section generalizes the BGD estimator proposed in Section (ref) to the semiparametric setting where $G$ is unknown.

In this setup, to ensure identification we set $\beta_{0}^{\star}$ to be 1, so our estimation target is $\boldsymbol{\beta}^{\star}$. To simplify our notation, we denote the space of $\mathbf{X}$ as $\mathcal{X}$, and the corresponding parameter space of $\boldsymbol{\beta}$ as $\mathcal{B}$. Suppose that an initial guess for $\boldsymbol{\beta}^{\star}$ is given by $\boldsymbol{\beta}_{1}$. In the $k$-th round of iteration, to update $\boldsymbol{\beta}$ based on the BGD algorithm, we require the knowledge of $G$ as in Section (ref), which is infeasible when $G$ is unknown. A natural idea is that we can construct an estimator for $G$ based on the index constructed from the updated parameter in the previous round. More intuitively, suppose for a moment that in the $k$-th round of iteration, $\boldsymbol{\beta}_{k}$ happens to be identical to the unknown true parameter $\boldsymbol{\beta}^{\star}$, then we have that $ G\left(z\right) =\mathbb{E}\left[\left.y\right|X_{0}+\boldsymbol{\mathbf{X}}^{\mathrm{T}}\boldsymbol{\beta}^{\star}=z\right]=\mathbb{E}\left[\left.y\right|X_{0}+\boldsymbol{\mathbf{X}}^{\mathrm{T}}\boldsymbol{\beta}_{k}=z\right] $ for any $z\in R$.

This motivates semiparametric estimation by using nonparametric methods to estimate $G\left(\cdot\right)$. We consider kernel estimation and the method of sieves in each of the following subsections.

The KBGD Estimator

In this section we consider tkernel estimation to estimate $G\left(\cdot\right)$. The Nadaraya-Watson kernel estimator of $G\left(\cdot\right)$ is of the form

where $K_{h}\left(\cdot\right)=h^{-1}K\left(\cdot/h\right)$, $K\left(\cdot\right)$ is some kernel function, and $h_{n}$ is some bandwidth parameter depending on $n$. Given the estimated CDF $\widehat{G}\left(\left.\cdot\right|\boldsymbol{\beta}_{k}\right)$, we can update the parameter as if it were the true CDF $G\left(\cdot\right)$. In particular, $\boldsymbol{\beta}_{k}$ is updated as

Keep updating $\boldsymbol{\beta}_{k}$ based on ((ref)) and ((ref)), until some terminating conditions are reached. The resulting estimator is labeled as the kernel-based batch gradient descent estimator (KBGD estimator).

For any fixed $z$ and $\boldsymbol{\beta}$, under mild conditions there holds $ \widehat{G}\left(\left.z\right|\boldsymbol{\beta}\right)\rightarrow_{p}\mathbb{E}\left[\left.y\right|X_{0}+\mathbf{X}^{\mathrm{T}}\boldsymbol{\beta}=z\right]. $ Denote such limit as $L\left(z,\boldsymbol{\beta}\right)$. Obviously, $L\left(z,\boldsymbol{\beta}^{\star}\right)=G\left(z\right)$ holds for any $z\in\mathbb{R}$. Before we move to a formal description of the statistical properties of the KBGD estimator based on ((ref)), we first provide some further discussion on $L\left(z,\boldsymbol{\beta}\right)$. For simplicity, in the following we only focus on the case where all the covariates are continuous which permit continuous joint density function. We leave further discussion of the case where some covariates are discrete to (ref). We point that when there are discrete covariates, our algorithm can be directly applied without any modification, although some further assumptions will be required.

When all the covariates are continuous, denote the joint density of $\mathbf{X}_{e}$ and $\mathbf{X}$ as $f_{e}\left(\mathbf{X}_{e}\right)=f_{e}\left(X_{0},\mathbf{X}\right)$ and $f\left(\mathbf{X}\right)=\int f_{e}\left(X_{0},\mathbf{X}\right)dX_{0}$, respectively. Denote $z\left(\mathbf{X}_{e},\boldsymbol{\beta}\right)=X_{0}+\mathbf{X}^{\mathrm{T}}\boldsymbol{\beta}$. Also denote $f_{\mathbf{X},z}\left(\left.\mathbf{X},z\right|\boldsymbol{\beta}\right)$ as the joint density of $\mathbf{X}$ and $z\left(\mathbf{X}_{e},\boldsymbol{\beta}\right)$ given $\boldsymbol{\beta}$. Note that for any $\boldsymbol{x}$ and $z$,

This implies that the joint density of $\mathbf{X}$ and $z\left(\mathbf{X}_{e},\boldsymbol{\beta}\right)$ given $\boldsymbol{\beta}$ is given by

and the marginal density of $z\left(\mathbf{X}_{e},\boldsymbol{\beta}\right)$ is given by

Define $f_{\mathbf{X}|z}\left(\left.\mathbf{X}\right|z,\boldsymbol{\beta}\right)=f_{\mathbf{X},z}\left(\left.\mathbf{X},z\right|\boldsymbol{\beta}\right)/f_{z}\left(\left.z\right|\boldsymbol{\beta}\right)$ as the conditional density of $\boldsymbol{\mathbf{X}}$ given $z$ and $\boldsymbol{\beta}$, we have that

where $\Delta\boldsymbol{\beta}=\boldsymbol{\beta}-\boldsymbol{\beta}^{\star}$.

Based on the above notations, now we formally study the asymptotic properties of the KBGD estimator under increasing dimensions. We first introduce some further assumptions.

The following lemma will be useful in the proof of our theorem.

(ref) implies that $\frac{1}{n}\sum_{i=1}^{n}\widehat{G}\left(\left.Z\left(\mathbf{X}_{e,i},\boldsymbol{\beta}\right)\right|\boldsymbol{\beta}\right)\boldsymbol{\mathbf{X}}_{i}$ will be closer to $\mathbb{E}\left[L\left(z\left(\boldsymbol{\mathbf{X}}_{e,i},\boldsymbol{\beta}\right),\boldsymbol{\beta}\right)\boldsymbol{\mathbf{X}}_{i}\right]$ uniformly with respect to $\boldsymbol{\beta}$ as $n$ increases. Note that such uniform convergence results are free of trimming; we do not need to trim $\boldsymbol{\mathbf{X}}_{e,i}$ even when the density of $z\left(\boldsymbol{\mathbf{X}}_{e,i},\boldsymbol{\beta}\right)$ is small. So even when $\widehat{G}\left(\left.z\left(\boldsymbol{\mathbf{X}}_{e,i},\boldsymbol{\beta}\right)\right|\boldsymbol{\beta}\right)$ is a poor estimator for $L\left(z\left(\boldsymbol{\mathbf{X}}_{e,i},\boldsymbol{\beta}\right),\boldsymbol{\beta}\right)$ for some $\boldsymbol{\mathbf{X}}_{e,i}$ and $\boldsymbol{\beta}$, our results are still valid. While on the same time, the cost of not conducting any trimming is that our guaranteed convergence rate depends heavily on the dimensionality. As is required in (ref), the dimension $p$ must satisfy $p^{\frac{5p+1}{2\left(p+1\right)}}\psi^{\frac{1}{p+1}}\left(n,p, h_{n}\right)\rightarrow0$. Suppose that $p/n\rightarrow0$ and we choose $h_{n}=\left(\left(\log n\right)/n\right)^{1/6}$, we have that $\psi\left(n,p,h_{n}\right)\sim\left(\left(\log n\right)/n\right)^{1/3}$. This implies that when $p$ is fixed, the convergence rate in (ref) is $\left(\left(\log n\right)/n\right)^{1/3\left(p+1\right)}$. When $p$ increases with $n$, the dimension $p$ should satisfy $p\log p=O\left(\log n\right)$, implying that $p$ is allowed to increase only mildly with $n$. The restriction on $p$ basically comes from the fact that as $\mathbf{X}_{e,i}$ moves towards the boundary of $\mathcal{X}_{e}$, the density of random variable $z\left(\mathbf{X}_{e,i},\boldsymbol{\beta}\right)$ decreases faster towards zero given a larger $p$, which makes the convergence rate sensitive to the increase of $p$.

For notational simplicity, in the following we denote $z\left(\mathbf{X}_{e,i},\boldsymbol{\beta}_{k}\right)$ and $z\left(\mathbf{X}_{e,i},\boldsymbol{\beta}^{\star}\right)$ as $z_{i,k}$ and $z_{i}^{\star}$. Based on the results in (ref), we have that under all conditions as imposed in (ref), there holds

Note that $z_{i,k}=z_{i}^{\star}+\mathbf{X}_{i}^{\mathrm{T}}\Delta\boldsymbol{\beta}_{k}$ and $L\left(z_{i,k},\boldsymbol{\beta}_{k}\right)=\int_{\mathcal{X}}G\left(z_{i,k}-\mathbf{X}^{\mathrm{T}}\Delta\boldsymbol{\beta}_{k}\right)f_{\mathbf{X}|z}\left(\left.\mathbf{X}\right|z_{i,k},\boldsymbol{\beta}_{k}\right)d\mathbf{X}$, so $\left(L\left(z_{i,k},\boldsymbol{\beta}_{k}\right)-G\left(z_{i}^{\star}\right)\right)\cdot\mathbf{X}_{i}$ equals to

where the integration is understood to be element-wise. To further simplify our notation, define \[ W\left(\boldsymbol{\mathbf{X}}_{e},\widetilde{\boldsymbol{\mathbf{X}}}_{e},\boldsymbol{\beta}\right)=G^{\prime}\left(z\left(\boldsymbol{\mathbf{X}}_{e},\boldsymbol{\beta}^{\star}\right)+\left(\boldsymbol{\mathbf{X}}-\widetilde{\boldsymbol{\mathbf{X}}}\right)^{\mathrm{T}}\Delta\boldsymbol{\beta}\right)f_{\boldsymbol{X}|z}\left(\left.\widetilde{\boldsymbol{\mathbf{X}}},\right|z\left(\boldsymbol{\mathbf{X}}_{e},\boldsymbol{\beta}\right),\boldsymbol{\beta}\right), \] \[ V\left(\boldsymbol{\mathbf{X}}_{e},\widetilde{\boldsymbol{\mathbf{X}}}_{e},\boldsymbol{\beta}\right)=\left(\boldsymbol{\boldsymbol{\mathbf{X}}}\boldsymbol{\boldsymbol{\mathbf{X}}}^{\mathrm{T}}-\boldsymbol{\boldsymbol{\mathbf{X}}}\widetilde{\boldsymbol{\boldsymbol{\mathbf{X}}}}^{\mathrm{T}}\right)W\left(\boldsymbol{\mathbf{X}}_{e},\widetilde{\boldsymbol{\mathbf{X}}}_{e},\boldsymbol{\beta}\right), \] and \[ \varLambda\left(\boldsymbol{\beta}\right)=\mathbb{E}\left[\int_{\mathcal{X}}V\left(\boldsymbol{\mathbf{X}}_{e,i},\mathbf{X}_{e},\boldsymbol{\beta}\right)d\mathbf{X}\right], \] we have that \[ \mathbb{E}\left[\left(L\left(z_{i,k},\boldsymbol{\beta}_{k}\right)-G\left(z_{i}^{\star}\right)\right)\cdot\mathbf{X}_{i}\right]=\int_{0}^{1}\varLambda\left(\boldsymbol{\beta}^{\star}+t\Delta\boldsymbol{\beta}_{k}\right)\Delta\boldsymbol{\beta}_{k}dt, \] which indicates that \[ \Delta\boldsymbol{\beta}_{k+1}=\left\{ \int_{0}^{1}\left(I_{p}-\delta_{k}\varLambda\left(\boldsymbol{\beta}^{\star}+t\Delta\boldsymbol{\beta}_{k}\right)\right)dt\right\} \Delta\boldsymbol{\beta}_{k}+\delta_{k}\cdot\left(\text{small order terms}\right). \] To ensure that with probability going to 1 the above iteration shrinks $\left\Vert \Delta\boldsymbol{\beta}_{k}\right\Vert $, we make the following assumption.

Based on the above assumptions, we have the following result.

(ref) implies that the iterative estimator based on ((ref)) and ((ref)) is consistent under increasing dimensions, no matter whether the starting point is close to the unknown true parameter or not. However, the convergence speed heavily depends on the dimensionality of the problem, $p$, even when $p$ is fixed. This is not ideal under our single-index setup but is not surprising since our algorithm does not involve any trimming procedure as we have discussed in (ref).

We proceed to establish the asymptotic normality of the KBGD estimator. Due to technical difficulties, throughout the following analysis in this section we only consider the case where $p$ is fixed. As we can see in (ref), even in the case of fixed dimensionality, the guaranteed convergence rate of the KBGD estimator based on ((ref)) and ((ref)) is at best $\left(\left(\log n\right)/n\right)^{\frac{1}{3p+3}}$, which still depends on $p$. To obtain asymptotic normality, we need to slightly modify our algorithm to get rid of the dependence on dimensionality. In particular, we introduce trimming to our algorithm. When updating the parameter, we only use observations that fall into a pre-selected region as did in ichimura1993semiparametric. In particular, the algorithm is modified as,

where $\widehat{G}\left(\left.z_{i,k}\right|\boldsymbol{\beta}_{k}\right)=\widehat{G}\left(\left.z\left(\boldsymbol{\mathbf{X}}_{e,i},\boldsymbol{\beta}_{k}\right)\right|\boldsymbol{\beta}_{k}\right)$ is defined in ((ref)), $I_{i}^{\phi}=I\left(\boldsymbol{\mathbf{X}}_{e,i}\in\mathcal{X}_{e}^{\phi}\right)$, and $\mathcal{X}_{e}^{\phi}$ is a subset of $\mathcal{X}_{e}$ given by

for some $\phi>0$ whose value will be determined later. Different from ((ref)), the update of $\boldsymbol{\beta}_{k}$ based on ((ref)) uses only a subset of the whole sample for which the covariate vector $\boldsymbol{\mathbf{X}}_{e,i}$ falls into $\mathcal{X}_{e}^{\phi}$. The reason why we choose the trimming set as in ((ref)) is that, as we show in the (ref), for any $0<\phi<1$, there holds $ \inf_{\left(\mathbf{X}_{e},\boldsymbol{\beta}\right)\in\mathcal{X}_{e}^{\phi}\times\mathcal{B}}f_{z}\left(\left.z\left(\mathbf{X}_{e},\boldsymbol{\beta}\right)\right|\boldsymbol{\beta}\right)\geq C\phi^{p}p^{-p} $ for some constant $C>0$ that depends on $\phi$. When $p$ and $\phi$ are both fixed, $f_{z}\left(\left.z\left(\boldsymbol{\mathbf{X}}_{e},\boldsymbol{\beta}\right)\right|\boldsymbol{\beta}\right)$ is uniformly lower bounded from zero for any combination $\left(\mathbf{X}_{e},\boldsymbol{\beta}\right)\in\mathcal{X}_{e}^{\phi}\times\mathcal{B}$, so the uniform estimation accuracy of $L\left(z\left(\boldsymbol{\mathbf{X}}_{e,i},\boldsymbol{\beta}\right),\boldsymbol{\beta}\right)$ over $\boldsymbol{\mathbf{X}}_{e,i}$ and $\boldsymbol{\beta}$ will be improved. Note that trimming will cause some efficiency loss by dropping some observations, but such loss can be controlled to be small if we choose $\phi$ to be close to zero. We also point that trimming is only applied to the update of the parameter; when nonparametrically estimating $G$, we still use all the data points.

To simplify our following notation, given the trimming parameter $\phi$, we denote $I^{\phi} \cdot \mathbf{X}$ as $\boldsymbol{\mathbf{X}}^{\phi}$. We also define \[ \varLambda_{\phi}\left(\boldsymbol{\beta}\right)=\mathbb{E}\left[I_{i}^{\phi}\cdot\int_{\mathcal{X}}V\left(\boldsymbol{\mathbf{X}}_{e,i},\boldsymbol{\mathbf{X}}_{e},\boldsymbol{\beta}\right)d\boldsymbol{\mathbf{X}}\right]. \] The following theorem provides a counterpart to the results in (ref).

Note that when $p$ is fixed, $\psi\left(n,p,h_{n}\right)$ no longer depends on $p$ asymptotically. The improvement over the convergence rate basically comes from the improvement of the uniform convergence rate of the kernel estimator due to trimming. Also note that under trimming, the minimum number of iteration in (ref)(i), $\widetilde{k}_{1,n}^{KBGD}$, is of order $\log n$ as long as $nh_{n}\rightarrow\infty$. This implies that under trimming, a faster convergence rate is guaranteed with the minimum number of iterations being of the same magnitude as that of the estimator without trimming.

We now proceed to establish the asymptotic normality of $\boldsymbol{\beta}_{k}$. Define

\[ \boldsymbol{\xi}_{n}^{\phi}=\frac{1}{n}\sum_{i=1}^{n}\left(\widehat{G}\left(\left.z_{i}^{\star}\right|\boldsymbol{\beta}^{\star}\right)-y_{i}\right)\mathbf{X}_{i}^{\phi}. \] We note that

where the integration is understood to be element-wise. To understand the properties of the above algorithm, we need the following lemmas.

Now we are in a position to illustrate the results of the asymptotic normality of our KBGD estimator.

We introduce the estimator for the variance matrix, based on which the confidence interval of $\boldsymbol{\beta}^{\star}$ can be then constructed.

We finally provide some remarks for the KBGD estimators.

The SBGD Estimator

In the previous section, we introduced the KBGD algorithm, where the update of the parameter is based on a BGD-type procedure while the unknown CDF is replaced with its Nadaraya-Watson kernel estimator constructed by the initial parameter. In this section, we consider an alternative nonparametric approximation for the unknown CDF based on the method of sieves. Given a set of basis functions $\{r_{j}\left(z\right)\}_{j=0}^{\infty}$ that is complete in $C\left(\mathbb{R}\right)$ space, any smooth CDF $G$ can be represented by $G\left(z\right)=\sum_{j=0}^{\infty}\pi_{j}^{\star}r_{j}\left(z\right)$ for any $z\in R$, where $\{\pi_{j}^{\star}\}_{j=0}^{\infty}$ is the unknown coefficients of the basis functions. In practice, to make our algorithm tractable, we truncate the sequence of the basis functions and only use the first $q+1$ basis functions for approximation, where $q$ increases with sample size $n$ at some rate. To approximate $G$, it then remains to provide an estimator for the unknown coefficients of the basis functions $\{\pi_{j}^{\star}\}_{j=0}^{q}$. Our estimation procedure for $\{\pi_{j}^{\star}\}_{j=0}^{q}$ shares similar intuition as the one that motivates the Nadaraya-Watson kernel estimator in the previous section. In particular, suppose for a moment that in the $k$-th round of update, we start with $\boldsymbol{\beta}_{k}$, which happens to be identical to the unknown true parameter $\boldsymbol{\beta}^{\star}$. In this case, define $\boldsymbol{r}_{q}(z)=\left(r_{0}\left(z\right),\cdots,r_{q}\left(z\right)\right)^{\mathrm{T}}$ and $\boldsymbol{\pi}_{q}^{\star}=\left(\pi_{1}^{\star},\cdots,\pi_{q}^{\star}\right)^{\mathrm{T}}$, we have that \[ y_{i}=G\left(z_{i,k}\right)+\varepsilon_{i}\approx\boldsymbol{r}_{q}^{\mathrm{T}}\left(z_{i,k}\right)\boldsymbol{\pi}_{q}^{\star}+\varepsilon_{i}, \] where recall that $z_{i,k}=X_{0,i}+\mathbf{X}_{i}^{\mathrm{T}}\boldsymbol{\beta}_{k}$. The above relationship motivates the following OLS estimator for the sieve coefficients

Given the estimator of the sieve coefficients $\widehat{\boldsymbol{\pi}}_{q,n,k}$, the unknown CDF $G$ in the $k$-th round of update is approximated by

Based on the estimated CDF $\widehat{G}\left(\left.z\right|\boldsymbol{\beta}_{k}\right)$, the update of the parameter can be carried out based on ((ref)). We iterate sequentially based on ((ref)), ((ref)) and ((ref)) until some terminating conditions are satisfied. The resulting estimator is then labeled as the sieve-based batch gradient descent estimator (SBGD estimator).

Define $R_{q}\left(z\right)=G\left(z\right)-\boldsymbol{r}^{\mathrm{T}}\left(z\right)\boldsymbol{\pi}_{q}^{\star}$, $\Gamma_{q,n}\left(\boldsymbol{\beta}\right)=\frac{1}{n}\sum_{i=1}^{n}\boldsymbol{r}_{q}\left(X_{0,i}+\mathbf{X}_{i}^{\mathrm{T}}\boldsymbol{\beta}\right)\boldsymbol{r}_{q}^{\mathrm{T}}\left(X_{0,i}+\mathbf{X}_{i}^{\mathrm{T}}\boldsymbol{\beta}\right)$, $\Gamma_{q,n,k}=\Gamma_{q,n}\left(\boldsymbol{\beta}_{k}\right)$, and $\mathfrak{X}_{q,n}\left(z,\boldsymbol{\beta}\right)=\frac{1}{n}\sum_{i=1}^{n}\left(\boldsymbol{r}_{q}^{\mathrm{T}}\left(X_{0,i}+\mathbf{X}_{i}^{\mathrm{T}}\boldsymbol{\beta}\right)\Gamma_{q,n}^{-1}\left(\boldsymbol{\beta}\right)\boldsymbol{r}_{q}\left(z\right)\mathbf{X}_{i}\right).$ Through tedious algebra, we can show that the SBGD procedure has the following representation,

where recall that $z_{i}^{\star}=X_{0,i}+\mathbf{X}_{i}^{\mathrm{T}}\boldsymbol{\beta}^{\star}$. To study the properties of the above procedure, we introduce some additional assumptions.

For any $-\infty<z<\infty$, define the population counterpart of $\mathfrak{X}_{q,n}\left(z,\boldsymbol{\beta}\right)$ as

\[ \mathfrak{X}_{q}\left(z,\boldsymbol{\beta}\right)=\mathbb{E}\left(\boldsymbol{r}_{q}^{\mathrm{T}}\left(z\left(\mathbf{X}_{e},\boldsymbol{\beta}\right)\right)\Gamma_{q}^{-1}\left(\boldsymbol{\beta}\right)\boldsymbol{r}_{q}\left(z\right)\mathbf{X}\right). \] Then we have the following lemma.

Obviously, (ref) provides a parallel result to ((ref)). In particular, define \[ \Psi_{q}\left(t,\boldsymbol{\beta}\right)=\mathbb{E}\left[G^{\prime}\left(z\left(\mathbf{X}_{e},\boldsymbol{\beta}^{\star}\right)+t\mathbf{X}^{\mathrm{T}}\Delta\boldsymbol{\beta}\right)\left(\mathbf{X}\mathbf{X}^{\mathrm{T}}-\mathfrak{X}_{q}\left(z\left(\mathbf{X}_{e},\boldsymbol{\beta}\right),\boldsymbol{\beta}\right)\mathbf{X}^{\mathrm{T}}\right)\right], \] under all the conditions imposed in (ref), we have that

Obviously, ((ref)) is also a parallel result to ((ref)). As a result, to ensure that ((ref)) actually constitutes a contraction for $\left\Vert \Delta\boldsymbol{\beta}_{k}\right\Vert $, we impose the following assumption that is similar to (ref).

Based on the above assumptions, we have the following result.

According to (ref), when $\chi_{2,n} \rightarrow 0$ as $n\rightarrow\infty$, the SBGD estimator is consistent as long as the number of updates exceeds $k_{1,n}^{SBGD}$. Based on such consistent estimator, we are ready to establish the asymptotic normality of our SBGD estimator. Apply the mean value theorem to ((ref)), we have that

Define $\Psi_{q}^{\star}=\mathbb{E}\left[G^{\prime}\left(z\left(\mathbf{X}_{e},\boldsymbol{\beta}^{\star}\right)\right)\left(\mathbf{X}\mathbf{X}^{\mathrm{T}}-\mathfrak{X}_{q}\left(z\left(\mathbf{X}_{e},\boldsymbol{\beta}^{\star}\right),\boldsymbol{\beta}^{\star}\right)\mathbf{X}^{\mathrm{T}}\right)\right]$ and $\mathfrak{V}_{q}=\mathbb{E}\left(\mathbf{X}_{i}\boldsymbol{r}_{q}^{\mathrm{T}}\left(z_{i}^{\star}\right)\Gamma_{q}^{-1}\left(\boldsymbol{\beta}^{\star}\right)\right)$. Similar to (ref) and (ref), we provide two additional lemmas that are useful to understand the above algorithm.

Based on the above two lemmas, we are now ready to study the asymptotic distribution of the SBGD estimator.

We now provide the estimator for the variance.

We finally provide some remarks on the empirical applications of the SBGD estimator.

Monte Carlo Experiments

This section conducts Monte Carlo simulations to study the performance of our KBGD and SBGD estimators. We focus on two aspects of our estimators. First we study the finite-sample properties of the KBGD estimator, including the bias and the root mean squared error (RMSE). Let the $j$-th argument of the true parameter be $\beta_{j}^{\star}$, and the simulation is repeated $R$ times, where its estimator in the $r$-th round of simulation is $\widehat{\beta}_{j}^{r}$, then the bias and RMSE are respectively given by $\text{Bias}=|\frac{1}{R}\sum_{r=1}^{R}(\widehat{\beta}_{j}^{r}-\beta_{j}^{\star})|$ and $\text{RMSE}=\sqrt{\sum_{r=1}^{R}(\widehat{\beta}_{j}^{r}-\beta_{j}^{\star})^{2}/R}.$ We also investigate whether the confidence interval based on the asymptotic distribution has good coverage rate. We consider nominal coverage rate $\alpha=0.95$, so the confidence interval for $\beta_{j}^{\star}$ in the $r$-th round of repetition is given by $CI_{j}^{r}=[\widehat{\beta}_{j}^{r}-1.96\cdot\widehat{\text{std}}_{j}^{r},\widehat{\beta}_{j}^{r}+1.96\cdot\widehat{\text{std}}_{j}^{r}]$, where $\widehat{\text{std}}_{j}^{r}$ is the estimated standard deviation of $\widehat{\beta}_{j}^{r}$. The actual coverage rate is then given by $CR=\frac{1}{R}\sum_{r=1}^{R}I(\beta_{j}^{\star}\in CI_{j}^{r}).$

We are also interested in how sensitive our estimators are to the initial guess of the true parameter. In each repetition of our simulation, we consider three different initial guesses: the true parameter vector, the parameter vector estimated based on the Logit regression, and the parameter with all elements being zeros. If the estimation results starting from different initial guesses are close or even identical to each other, the estimation methods are insensitive to the initial guesses and thus are robust in terms of computation. Denote $\widehat{\boldsymbol{\beta}}_{T}^{r}$, $\widehat{\boldsymbol{\beta}}_{L}^{r}$, and $\widehat{\boldsymbol{\beta}}_{Z}^{r}$ as the estimators with starting points being true parameter, Logit estimator, and vector of zeros. We use $S_{L}=\sqrt{\frac{1}{R}\sum_{i=1}^{n}||\widehat{\boldsymbol{\beta}}_{L}^{r}-\widehat{\boldsymbol{\beta}}_{T}^{r}||^{2}}$ and $S_{Z}=\sqrt{\frac{1}{R}\sum_{i=1}^{n}||\widehat{\boldsymbol{\beta}}_{Z}^{r}-\widehat{\boldsymbol{\beta}}_{T}^{r}||^{2}}$ as the measurement of the sensitivity.

We consider data generating process $y_{i}=I(X_{0,i}+\beta_{1}^{\star}X_{1,i}\cdots+\beta_{10}^{\star}X_{10,i}-u_{i}>0),i=1,2,\cdots,n,$ where data are i.i.d over $i$, and $X_{0,i},X_{1,i},\cdots,X_{10,i},u_{i}$ are also independent. We set $\boldsymbol{\beta}^{\star}=(1,0.5,-0.5,1,-1,2,-2,4,-4,1.5,-1.5)^{\mathrm{T}}$, $X_{j,i}\sim N\left(0,1\right)$ for $0\leq j\leq8$, $X_{9,i}\sim\text{Bernoulli}\left(1/2\right)$, $X_{10,i}\sim\text{Poisson}\left(2\right)$, and $u_{i}\sim Cauchy$. We consider two sample sizes $n=2500$ and $5000$. Finally, for finite-sample performance, we repeat the simulation 500 times; for sensitivity analysis, we repeat 100 times.

(ref) reports the finite-sample properties of our estimators. It can be seen that our estimators works well in finite sample cases. Both estimators have small bias, whose RMSE decrease with sample size. Moreover, the confidence interval constructed based on the asymptotic variance and normal approximation has actual coverage rate that is quite close to the nominal rate $0.95$.

(ref) reports the sensitivity of our estimators to the starting points. We can see that for both estimators, $S_{L}$ and $S_{Z}$ are close to zero, indicating that the resulting estimators starting from Logit estimator or zeros are almost identical to the ones starting from the unknown true parameter. Such a result demonstrates that our algorithms are robust to different initial guesses. We also find that compared with SBGD, our KBGD estimator takes much longer time to converge.

The robustness of our algorithm might be sensitive to the setups of coefficients. To check whether this is the case, instead of using the fixed parameters specified before, in each round of simulation we randomly draw true parameter $\boldsymbol{\beta}^{\star}$ as follows $\beta_{1}^{\star},\beta_{2}^{\star},\beta_{9}^{\star},\beta_{10}^{\star}\sim N\left(0,1\right)$, $\beta_{3}^{\star},\beta_{4}^{\star},\beta_{5}^{\star},\beta_{6}^{\star}\sim2N\left(0,1\right)$, and $\beta_{7}^{\star},\beta_{8}^{\star}\sim4N\left(0,1\right)$. The simulation results are reported in (ref). We can see that the results are similar to those under fixed parameters, indicating that our algorithm is robust to initial point under different parameter setups.

Empirical Application

As an empirical illustration of our new methods, this section applies our KBGD and SBGD estimation procedures to study how education affects the risk aversion. In the existing researches, it's extensively documented that, on the individual level, risk aversion is significantly correlated with the level of education, although the directions of correlation are mixed, see outreville2015relationship for a comprehensive review. In this study, we investigate how educational background of the family affects the risk aversion of the household as well as household-level investing behaviors. We use the national survey data from 2019 China Household Financial Survey Project (CHFS) gan2014data, which provides household-level information over demographics, asset and debt, income and consumption, social security and insurance, and various household's subjective preferences. The dependent variable we are interested in is the degree of risk aversion of the household. In particular, $y_{i}$ is constructed to take value of $0$ if the $i$-th household is completely against any form of risks and thus is described as being extremely risk averse; it takes value of $1$ if the family is willing to bear some form of risks when making investments. We study how the probability of $y_{i}=1$ is affected by a set of factors based on the binary choice model. The key factor that we are particularly interested in is the educational backgrounds, which is defined as the year of education of the head of the household. We also consider a set of other control variables including gender, ethnicity, health conditions, marital status, region of residence, economic knowledge, total income and total asset, whose impacts on the risk aversion are of interest on their own right. See Yao (2023) for detailed discussion on the construction of the data sets.

Before estimation, we normalize all the continuous variables so that the resulting variables all have zero mean and unity variance. To provide a comparison to the semiparametric estimation results, we first conduct parametric Logit regression and report the normalized coefficients in regression (I) in (ref). We then conduct KBGD and SBGD estimation and report the estimated coefficients of education in (II) and (III). As we can see from (ref), no matter which estimation methods we use, the coefficient of educational background is estimated to be positive with significance at $1\%$ level. This implies that, holding other conditions fixed, on average an increase in the year of education of the head in the households leads to the increase of willingness to bear risks. Comparing the semiparametric estimation results with that of Logit regression, we can see that the KBGD and SBGD estimators are close to each other, which are both smaller than that of Logit regression, indicating that parametric estimation might suffer from model misspecification and lead to an overestimation of the impacts of education on risk aversion. We finally compare the computation time of each method. We can see that both KBGD and SBGD estimators take much longer to converge compared with the parametric estimation. Comparatively, the SBGD algorithm is significantly faster than the KBGD algorithm, which takes over two hours to converge. This result supports the use of SBGD algorithm when there are data of large scale.

Conclusions

\setcounter{equation}{0} In this paper, we proposed new estimation procedures for binary choice and monotonic index models with increasing dimensions. Existing semiparametric estimation procedures for this model cannot be implemented in practice when the number of regressors is large. In contrast, our algorithmic based procedures can be used for many regressor models as it involves convex optimization at each iteration of the procedure. We show this iterative procedure also has desirable asymptotic properties when the number of regressors increases with the sample size in ways that are standard in big data literature.