The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
103,584 characters
Estimating High Dimensional Monotone Index Models by Iterative Convex Optimization
\maketitle
\vskip -.5in
{\footnotesize \singlespace
\begin{abstract}
In this paper we propose new approaches to estimating large dimensional monotone index models. This class of models has been popular in the applied and theoretical econometrics literatures as it includes discrete choice, nonparametric transformation, and duration models. A main advantage of our approach is computational. For instance, rank estimation procedures such as those proposed in \cite{han1987non} and \cite{cavanagh1998rank} that optimize a nonsmooth, non convex objective function are difficult to use with more than a few regressors and so limits their use in with economic data sets. For such monotone index models with increasing dimension, we propose to use a new class of estimators based on {\em batched gradient descent (BGD) } involving nonparametric methods such as kernel estimation or sieve estimation, and study their asymptotic properties. The BGD algorithm uses an iterative procedure where the key step exploits a strictly convex objective function, resulting in computational advantages. A contribution of our approach is that our model is large dimensional and semiparametric and so does not require the use of parametric distributional assumptions.
\end{abstract}
}
\noindent
{\small
}
\
\noindent
{\bf Key Words} Monotone Index models, Convex Optimization, Kernel and Sieve Estimation.
\thispagestyle{empty}
\newpage
\setcounter{page}{1}
\setlength{\abovedisplayshortskip}{5pt}
\setlength{\belowdisplayshortskip}{5pt}
\setlength{\abovedisplayskip}{5pt}
\setlength{\belowdisplayskip}{5pt}
\setcounter{equation}{0}
\section{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:
\begin{equation} y_i=I[x_i'\beta^\ast_e -u_i \geq 0 ] \end{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 \cite{han1987non}, \cite{ichimura1993semiparametric}, \cite{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 \cite{fanlvqi2011}. Recent papers include
\cite{neweywind2009},
\cite{chernozhukov2017central},\cite{bellonietal2018}, \cite{cattaneoetal2018a}, \cite{cattaneoetal2018b},
Related to our work is the recent literature on estimating large dimensional binary choice or monotone index models in \cite{sur2019modern} and \cite{fan2020rank}.
\cite{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 \cite{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. \cite{fan2020rank} and \cite{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 \cite{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
\section{\label{section2}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 \textit{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.
\begin{assumption}
\label{assu1}An i.i.d. data set $\mathscr{D}_{n}=\left\{ \left(\mathbf{X}_{e,i},y_{i}\right)\right\} _{i=1}^{n}$
of sample size $n$ is observed, where $y_{i}$ is generated \footnote{Here we are decomposing the vector $\mathbf{X}_{e,i}$ into a scalar component
$X_{0,i}$ and the vector $\mathbf{X}_i$, and decomposing the vector of parameters $\boldsymbol{\beta}^{\star}_e$ into the scalar term $\beta_0^\star$ and the vector $\boldsymbol{\beta}^{\star}$. As we will see this is done for notational convenience when imposing scale normalizations. } by $
y_i=I\left(X_{0,i}\beta_{0}^{\star}+\mathbf{X}_i^{\mathrm{T}}\boldsymbol{\beta}^{\star}-u_i>0\right)$ with unobserved shock $u_{i}$ that is independent
of $\mathbf{X}_{e,i}$ and has CDF $G\left(\cdot\right)$.
\end{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,
\begin{equation}
\boldsymbol{\beta}_{e,k+1}=\boldsymbol{\beta}_{e,k}-\frac{\delta_{k}}{n}\sum_{i=1}^{n}\partial\ell_{G}\left(\boldsymbol{\beta}_{e,k},\mathbf{X}_{e,i},y_{i}\right)/\partial\boldsymbol{\beta}_{e},\label{known_loss}
\end{equation}
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{known_loss}) until some terminating conditions are reached.
In this paper, we consider the following loss function
\begin{equation}
\ell_{G}\left(\boldsymbol{\beta}_{e},\mathbf{X}_{e},y\right)=\int_{-A}^{\mathbf{X}_{e}^{\mathrm{T}}\boldsymbol{\beta}_{e}}G\left(z\right)dz-y\mathbf{X}_{e}^{\mathrm{T}}\boldsymbol{\beta}_{e},\label{loss_function}
\end{equation}
for some sufficiently large positive constant $A$. The loss function (\ref{loss_function}) was also considered in \citet{agarwal2014least}
and has many nice properties. For instance, under some mild conditions, we can show that
\begin{align*}
\frac{\partial\mathbb{E}\left(\ell_{G}\left(\boldsymbol{\beta}_{e}^{\star},\mathbf{X}_{e},y\right)\right)}{\partial\boldsymbol{\beta}_{e}} & =\mathbb{E}\left\{ \left(G\left(\mathbf{X}_{e}^{\mathrm{T}}\boldsymbol{\beta}_{e}^{\star}\right)-y\right)\mathbf{X}_{e}\right\} \\
& =\mathbb{E}\left\{ \left(G\left(\mathbf{X}_{e}^{\mathrm{T}}\boldsymbol{\beta}_{e}^{\star}\right)-\mathbb{E}\left(\left.y\right|\mathbf{X}_{e}\right)\right)\mathbf{X}_{e}\right\} =0,
\end{align*}
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{loss_function})
is that the derivative of (\ref{loss_function}) 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{loss_function}), the BGD estimator is
obtained based on the following iteration
\begin{equation}
\boldsymbol{\beta}_{e,k+1}=\boldsymbol{\beta}_{e,k}-\frac{\delta_{k}}{n}\sum_{i=1}^{n}\left(G\left(\mathbf{X}_{e,i}^{\mathrm{T}}\boldsymbol{\beta}_{e,k}\right)-y_{i}\right)\mathbf{X}_{e,i}.\label{known_G}
\end{equation}
\ \
\begin{remark}
Key to the above approach is the construction of a convex objective function that facilitates computation even with high dimensions. This transformed convex objective works for any monotone model. In particular, for any model of the form $y_i =G(x_i'\beta) + \epsilon$ with $E[\epsilon_i|x_i]=0$ and monotone $G(.)$, a similar {\it convex} criterion as in (\ref{loss_function}) can be used for inference on $\beta.$
\end{remark}
We now describe the asymptotic properties of $\boldsymbol{\beta}_{e,k}$.
We first make the following assumption.
\begin{assumption}
\label{assump:2}(i) $\mathcal{X}_{e}=\left[-1,1\right]^{p+1}$; (ii)
$\mathcal{B}_{e}$ is convex, and there exists some constant $B_{0}>0$
such that for any $\boldsymbol{\beta}_{e}\in\mathcal{B}_{e}$, $\left|\beta_{j}\right|\leq B_{0}$
for any $0\leq j\leq p$; (iii) there exists integer $\upsilon_{G}$
such that $G$ has up to $\upsilon_{G}$-th bounded derivatives; (iv)
Define $M_{n}\left(\boldsymbol{\beta}_{e}\right)=\frac{1}{n}\sum_{i=1}^{n}G^{\prime}\left(\mathbf{X}_{e,i}^{\mathrm{T}}\boldsymbol{\beta}_{e}\right)\mathbf{X}_{e,i}\mathbf{X}_{e,i}^{\mathrm{T}}$
and $M\left(\boldsymbol{\beta}_{e}\right)=\mathbb{E}[M_{n}\left(\boldsymbol{\beta}_{e}\right)]$.
For any $\boldsymbol{\beta}_{e}\in\mathcal{B}_{e}$, there holds $0<\underline{\lambda}_{e}\leq\underline{\lambda}\left(M\left(\boldsymbol{\beta}_{e}\right)\right)\leq\overline{\lambda}\left(M\left(\boldsymbol{\beta}_{e}\right)\right)\leq\overline{\lambda}_{e}<\infty$.
\end{assumption}
\begin{remark}
\autoref{assump:2}(i) and \autoref{assump:2}(ii) are convenient
normalizations that facilitate the assessment of our model.
Note that to ensure that $\boldsymbol{\beta}_{e,k}$ falls into a
compact set for each $k$, some form of truncation on $\boldsymbol{\beta}_{e,k+1}$
in (\ref{known_G}) is needed. While according to our results below, as long as $\mathcal{B}_{e}$
is sufficiently large, it can be
shown that $\boldsymbol{\beta}_{e,k}$ will fall into $\mathcal{B}_{e}$
for all $k$ with probability going to 1. We then assume that
$\boldsymbol{\beta}_{e,k}\in\mathcal{B}_{e}$ for all $k$.
\autoref{assump:2}(iii) imposes some smoothness conditions on $G$, where the requirement on $\upsilon_G$ will be stated in the following propositions and theorems.
\autoref{assump:2}(iv) requires that the eigenvalue of $M_{n}\left(\boldsymbol{\beta}_{e}\right)$
is bounded from both below and above uniformly over $\mathcal{B}_{e}$.
\end{remark}
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 \autoref{assu1} and \autoref{assump:2} hold, we have the
following result.
\begin{theorem}
\label{prop:Known_G}Suppose that \autoref{assu1} and \autoref{assump:2}
hold with $\upsilon_{G}=3$,
that $p^{5}\left(\log p\right)^{2}n^{-1}\rightarrow0$, that the learning rate is chosen such that $\delta_{k}=\delta\leq2/\left(3\overline{\lambda}_{e}\right)$, and that $\boldsymbol{\beta}_e$ is updated under (\ref{known_G}). We have
that
(i) Define
\[
k_{1,n}^{BGD}=\frac{\log\left\Vert \Delta\boldsymbol{\beta}_{e,1}\right\Vert +\frac{1}{2}\log\left(n/\left(p\log p\right)\right)}{-\log\left(1-\underline{\lambda}_{e}\delta/2\right)},
\]
we then have
\[
\sup_{k\geq k_{1,n}^{BGD}+1}\left\Vert \Delta\boldsymbol{\beta}_{e,k}\right\Vert =O_{p}\left(\sqrt{p\left(\log p\right)/n}\right);
\]
(ii) Define $k_{2,n}^{BGD}$ such that $\left(1-\underline{\lambda}_e\delta\right)^{k_{2,n}^{BGD}}\sqrt{p\log p}\rightarrow 0$, we have
\[
\sup_{k\geq k_{2,n}^{BGD}+1}\left\Vert \Delta\boldsymbol{\beta}_{e,k+k_{1,n}^{BGD}}-M^{-1}\left(\boldsymbol{\beta}_{e}^{\star}\right)\frac{1}{n}\sum_{i=1}^{n}\varepsilon_{i}\mathbf{X}_{e,i}\right\Vert =o_{p}\left(1/\sqrt{n}\right);
\]
(iii) For any $k\geq k_{1,n}^{BGD} + k_{2,n}^{BGD} + 1$, define
$\widehat{\boldsymbol{\beta}}_{e}=\widehat{\boldsymbol{\beta}}_{k}$.
Also define
\[
\Sigma_{1}^{\star}=M^{-1}\left(\boldsymbol{\beta}_{e}^{\star}\right)\mathbb{E}\left[G_{i}^{\star}\left(1-G_{i}^{\star}\right)\boldsymbol{\mathbf{X}}_{e,i}\boldsymbol{\mathbf{X}}_{e,i}^{\mathrm{T}}\right]M^{-1}\left(\boldsymbol{\beta}_{e}^{\star}\right),
\]
and
\[
\widehat{\Sigma}_{1,n}=M_n^{-1}\left(\widehat{\boldsymbol{\beta}}_{e}\right)\left\{ \frac{1}{n}\sum_{i=1}^{n}\widehat{G}_{i}\left(1-\widehat{G}_{i}\right)\boldsymbol{\mathbf{X}}_{e,i}\boldsymbol{\mathbf{X}}_{e,i}^{\mathrm{T}}\right\} M_n^{-1}\left(\widehat{\boldsymbol{\beta}}_{e}\right),
\]
where $G_{i}^{\star}=G\left(\boldsymbol{\mathbf{X}}_{e,i}^{\mathrm{T}}\boldsymbol{\beta}_{e}^{\star}\right)$
and $\widehat{G}_{i}=G\left(\boldsymbol{\mathbf{X}}_{e,i}^{\mathrm{T}}\widehat{\boldsymbol{\beta}}_{e}\right)$. Suppose further that $\mathbb{E}\left(\mathbf{X}_{e,i}\mathbf{X}_{e,i}^{\mathrm{T}}\right)$ has uniformly (with respect to $p$) upper bounded eigenvalues,
there holds
\[
\left\Vert \widehat{\Sigma}_{1,n}-\Sigma_{1}^{\star}\right\Vert \rightarrow_{p} 0.
\]
(iv) For any $p+1$ vector $\rho$ such that $\lim_{n\rightarrow\infty}\left\Vert \rho\right\Vert <\infty$,
$\lim_{n\rightarrow\infty}\rho^{\mathrm{T}}\Sigma_{1}^{\star}\rho=\sigma^{2}\left(\rho\right)$,
and that $\rho^{\mathrm{T}}M^{-1}\left(\boldsymbol{\beta}_{e}^{\star}\right)\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\varepsilon_{i}\mathbf{X}_{e,i}\rightarrow_{d}N\left(0,\sigma^{2}\left(\rho\right)\right)$,
we have that
\[
\rho^{\mathrm{T}}\Delta\widehat{\boldsymbol{\beta}}_{e}/\sqrt{\widehat{\sigma}^2\left(\rho\right)/n}\rightarrow_{d}N\left(0,1\right),
\]
where $\widehat{\sigma}^{2}\left(\rho\right)=\rho^{\mathrm{T}}\widehat{\Sigma}_{1,n}\rho$.
\end{theorem}
\begin{proof}[Proof of \autoref{prop:Known_G}]
See \autoref{appendixB}.
\end{proof}
When $p$ is fixed, \autoref{prop:Known_G}(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 \autoref{prop:Known_G}(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 \autoref{prop:Known_G}(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 \autoref{rem4} at the end of Section \ref{section4}. The
inference on $\boldsymbol{\beta}_{e}^{\star}$ based on the BGD estimator
is given by \autoref{prop:Known_G}(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}$ \citep[e.g.,][]{chernozhukov2017central}.
Before we conclude this section and move to semiparametric estimation, we further comment on \autoref{prop:Known_G}.
Different from the stochastic gradient descent algorithm \citep[e.g.,][]{toulis2017asymptotic}, we show in \autoref{prop:Known_G}
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 \autoref{prop:Known_G}. In particular,
we have the following proposition.
\begin{theorem}
\label{prop2}Suppose that all the conditions in \autoref{prop:Known_G}
hold and that $\boldsymbol{\beta}_e$ is updated under (\ref{known_G}). For any sequence of tuning parameters $\left\{ \delta_{k}\right\} _{k=1}^{\infty}$
satisfying $\delta_{k}\geq0$, $\delta_{k}\rightarrow0$, $\limsup_{k\rightarrow\infty}\delta_{k-1}/\delta_{k}<\infty$,
and $\sum_{k=1}^{\infty}\delta_{k}=\infty$, we have that
(i) Define $\widetilde{k}_{1,n}^{BGD}$ such that
$
\sum_{k=1}^{\widetilde{k}_{1,n}^{BGD}}\delta_{k}\geq\underline{\lambda}_{e}^{-1}\left\{ \log\left(n/p\left(\log p\right)\right)+2\log\left\Vert \Delta\boldsymbol{\beta}_{e,1}\right\Vert \right\} ,
$
and that
$
\sup_{k\geq\widetilde{k}_{1,n}^{BGD}+1}\delta_{k}\leq2/\underline{\lambda}_{e},
$
then there holds
\[
\sup_{k\geq \widetilde{k}_{1,n}^{BGD}+1}\left\Vert \Delta\boldsymbol{\beta}_{e,k}\right\Vert =O_{p}\left(\sqrt{p\left(\log p\right)/n}\right);
\]
(ii) Define $\widetilde{k}_{2,n}^{BGD}$ such that $\sum_{k=\widetilde{k}_{1,n}^{BGD}+1}^{k=\widetilde{k}_{2,n}^{BGD}}\delta_k /\log p \rightarrow \infty$, then we have that
\[
\sup_{k\geq \widetilde{k}_{2,n}^{BGD}+1}\left\Vert \Delta\boldsymbol{\beta}_{e,k+\widetilde{k}_{1,n}^{BGD}}-M^{-1}\left(\boldsymbol{\beta}_{e}^{\star}\right)\frac{1}{n}\sum_{i=1}^{n}\varepsilon_{i}\mathbf{X}_{e,i}\right\Vert =o_{p}\left(1/\sqrt{n}\right);
\]
(iii) For any $k\geq \widetilde{k}_{1,n}^{BGD} + \widetilde{k}_{2,n}^{BGD} + 1$, define
$\widehat{\boldsymbol{\beta}}_{e}=\widehat{\boldsymbol{\beta}}_{k}$.
We have that \autoref{prop:Known_G}(iii) and (iv) hold.
\end{theorem}
\begin{proof}[Proof of \autoref{prop2}]
See \autoref{appendixB}.
\end{proof}
\autoref{prop2} 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 \autoref{prop:Known_G}, 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 \autoref{prop:Known_G}(i); when $\upsilon>0$,
we can see that more rounds of iteration is needed compared with required
in \autoref{prop:Known_G}(i).
\section{\label{section3}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 \autoref{assu1} 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{section2} 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 \textbf{$\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{section2}, 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.
\subsection{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
\begin{equation}
\widehat{G}\left(\left.z\right|\boldsymbol{\beta}_{k}\right)=\frac{\sum_{j=1}^{n}K_{h_{n}}\left(z-X_{0,j}-\boldsymbol{\mathbf{X}}_{j}^{\mathrm{T}}\boldsymbol{\beta}_{k}\right)y_{j}}{\sum_{j=1}^{n}K_{h_{n}}\left(z-X_{0,j}-\boldsymbol{\mathbf{X}}_{j}^{\mathrm{T}}\boldsymbol{\beta}_{k}\right)},z\in R,\label{kernel_estimator}
\end{equation}
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
\begin{equation}
\boldsymbol{\beta}_{k+1}=\boldsymbol{\beta}_{k}-\frac{\delta_{k}}{n}\sum_{i=1}^{n}\left(\widehat{G}\left(\left.X_{0,i}+\mathbf{X}_{i}^{\mathrm{\mathrm{T}}}\boldsymbol{\beta}_{k}\right|\boldsymbol{\beta}_{k}\right)-y_{i}\right)\mathbf{X}_{i}.\label{update}
\end{equation}
Keep updating $\boldsymbol{\beta}_{k}$
based on (\ref{kernel_estimator}) and (\ref{update}), until some
terminating conditions are reached. The resulting estimator is labeled as
the \textit{kernel-based batch gradient descent estimator} (KBGD estimator).
\begin{remark}
In essence, the KBGD estimator can not be classified as a BGD estimator
based on a semiparametric loss function. In the semiparametric setup,
given any loss function $\ell_G\left(\boldsymbol{\beta},\mathbf{X}_{e},y\right)$
(quadratic distance in \citet{ichimura1993semiparametric}, log-likelihood
in \citet{klein1993efficient}, or loss function given in (\ref{loss_function}))
with unknown function $G$, it's a common practice to replace $G$
with its nonparametric estimator $\widehat{G}$ and then minimize
(or maximize) the resulting loss function to obtain the
estimator of $\boldsymbol{\beta}$. Note that under the single-index framework, $\widehat{G}$ usually
involves the unknown parameter $\boldsymbol{\beta}$, which is
of the form $\widehat{G}\left(\cdot\right)=\widehat{G}\left(\left.\cdot\right|\boldsymbol{\beta}\right)$.
In this scenario, the BGD estimator is constructed by the following
iteration
\[
\boldsymbol{\beta}_{k+1}^{BGD}=\boldsymbol{\beta}_{k}^{BGD}-\frac{\delta_{k}}{n}\sum_{i=1}^{n}\frac{\partial\ell_{\widehat{G}\left(\left.\cdot\right|\boldsymbol{\beta}_{k}^{BGD}\right)}\left(\boldsymbol{\beta}_{k}^{BGD},\mathbf{X}_{e,i},y_{i}\right)}{\partial\boldsymbol{\beta}},
\]
where $\partial\ell_{\widehat{G}\left(\left.\cdot\right|\boldsymbol{\beta}_{k}^{BGD}\right)}\left(\boldsymbol{\beta}_{k}^{BGD},\mathbf{X}_{e,i},y_{i}\right)/\partial\boldsymbol{\beta}$
involves $\partial\widehat{G}\left(\left.\cdot\right|\boldsymbol{\beta}_k\right)/\partial\boldsymbol{\beta}$,
a complicated functions of $\boldsymbol{\beta}_k$. In particular, the BGD estimator under loss function (\ref{loss_function}) is given by
\[
\boldsymbol{\beta}_{k+1}=\boldsymbol{\beta}_{k}-\frac{\delta_{k}}{n}\sum_{i=1}^{n}\left(\widehat{G}\left(\left.X_{0,i}+\mathbf{X}_{i}^{\mathrm{\mathrm{T}}}\boldsymbol{\beta}_{k}\right|\boldsymbol{\beta}_{k}\right) + \int_{-\infty}^{X_{0,i}+\mathbf{X}_i^{\mathrm{T}}\boldsymbol{\beta}_k}\frac{\partial\widehat{G}\left(\left.z\right|\boldsymbol{\beta}_k\right)}{\partial \boldsymbol{\beta}}dz-y_{i}\right)\mathbf{X}_{i}.
\]
Obviously, an additional term is introduced compared with (\ref{update}). On the contrary,
during the construction (\ref{update}), we take $G$ as given
when taking the first order derivative of the loss function and then replace the unknown $G$
with its non-parametric estimator in the derivative. More specifically,
the KBGD estimator is updated as follows
\[
\boldsymbol{\beta}_{k+1}=\boldsymbol{\beta}_{k}-\frac{\delta_{k}}{n}\sum_{i=1}^{n}\left.\frac{\partial\ell_G\left(\boldsymbol{\beta}_{k},\mathbf{X}_{e,i},y_{i}\right)}{\partial\boldsymbol{\beta}}\right|_{G\left(\cdot\right)=\widehat{G}\left(\left.\cdot\right|\boldsymbol{\beta}_{k}\right)},
\]
so additional terms involving $\partial\widehat{G}\left(\left.\cdot\right|\boldsymbol{\beta}_k\right)/\partial\boldsymbol{\beta}$ are avoided. Finally, as we discussed in Section \ref{section2}, the derivative of loss function
(\ref{loss_function}) with respect to $\boldsymbol{\beta}$
depends only on $G$, so we also avoid approximating the derivative of $G$, which has poorer finite-sample performance compared with approximating $G$. Such update also ensures contraction map under some conditions, see \autoref{assu:5}.
\end{remark}
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{update}), 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
\autoref{rem5}. 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$,
\begin{align*}
P\left[\mathbf{X}\leq\boldsymbol{x},z\left(\mathbf{X}_{e},\boldsymbol{\beta}\right)\leq z\right] & =\int_{\widetilde{\mathbf{X}}\leq\boldsymbol{x},\widetilde{X}_{0}+\widetilde{\mathbf{X}}^{\mathrm{T}}\boldsymbol{\beta}\leq z}f_{e}\left(\widetilde{X}_{0},\widetilde{\mathbf{X}}\right)d\widetilde{X}_{0}d\widetilde{\mathbf{X}}\\
& =\int_{\widetilde{\mathbf{X}}\leq\boldsymbol{x}}\left[\int_{\widetilde{X}_{0}\leq z-\widetilde{\mathbf{X}}^{\mathrm{T}}\boldsymbol{\beta}}f_{e}\left(\widetilde{X}_{0},\widetilde{\mathbf{X}}\right)d\widetilde{X}_{0}\right]d\widetilde{\mathbf{X}}.
\end{align*}
This implies that the joint density of $\mathbf{X}$ and $z\left(\mathbf{X}_{e},\boldsymbol{\beta}\right)$
given $\boldsymbol{\beta}$ is given by
\begin{equation}
f_{\mathbf{X},z}\left(\left.\mathbf{X},z\right|\boldsymbol{\beta}\right)=f_{e}\left(z-\mathbf{X}^{\mathrm{T}}\boldsymbol{\beta},\mathbf{X}\right),\label{eq:joint_pdf_X_Z}
\end{equation}
and the marginal density of $z\left(\mathbf{X}_{e},\boldsymbol{\beta}\right)$
is given by
\begin{equation}
f_{z}\left(\left.z\right|\boldsymbol{\beta}\right)=\int_{\mathcal{X}}f_{\mathbf{X},z}\left(\left.\mathbf{X},z\right|\boldsymbol{\beta}\right)d\mathbf{X} = \int_{\mathcal{X}}f_{e}\left(z-\mathbf{X}^{\mathrm{T}}\boldsymbol{\beta},\mathbf{X}\right)d\mathbf{X}.\label{eq:marginal_pdf_Z}
\end{equation}
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
\begin{align}
L\left(z,\boldsymbol{\beta}\right) & =\mathbb{E}\left(\left.G\left(z-\boldsymbol{\mathbf{X}}^{\mathrm{T}}\Delta\boldsymbol{\beta}\right)\right|z\left(\boldsymbol{\mathbf{X}}_{e},\boldsymbol{\beta}\right)=z\right)\nonumber \\
& =\int_{\mathcal{X}}G\left(z-\boldsymbol{\mathbf{X}}^{\mathrm{T}}\Delta\boldsymbol{\beta}\right)f_{\mathbf{X}|z}\left(\left.\mathbf{X}\right|z,\boldsymbol{\beta}\right)d\mathbf{X},\label{eq:L}
\end{align}
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.
\begin{assumption}
\label{assu:3}The kernel function $K\left(\cdot\right)$ satisfies:
(i) $K$ is bounded and twice continuously differentiable with bounded
first and second derivatives, and the second derivative satisfies
Lipschitz condition on the whole real line; (ii) $\int K\left(s\right)ds=1$;
(iii) there exists positive integer $\upsilon_{K}$ such that $\int s^{\upsilon}K\left(s\right)du=0$
for $1\leq\upsilon\leq\upsilon_{K}-1$ and $\int u^{\upsilon_{K}}K\left(u\right)du\neq0$;
(iv) $K\left(s\right)=0$ for $\left|s\right|>1$.
\end{assumption}
\begin{assumption}
\label{assu:4}(i) There exists some constant $\zeta>1$ such that
$\zeta^{-1}\leq f_{e}\left(\mathbf{X}_{e}\right)\leq\zeta$ holds
for all $\mathbf{X}_{e}\in\mathcal{X}_{e}$; (ii) there exists positive
integer $\upsilon_{f}$ such that $f_{e}\left(\mathbf{X}_{e}\right)$
has bounded up to $\upsilon_{f}$-th derivatives.
\end{assumption}
\begin{remark}
\autoref{assu:4}(i) together with \autoref{assump:2}(i)
is a commonly-used assumption in the machine learning literature
\citep[e.g.,][]{wager2018estimation}. It basically requires that the joint density
of $\mathbf{X}_{e}$ is uniformly bounded from both above and below
over $\mathcal{X}_{e}$, so the density does not degenerate over $\mathcal{X}_{e}$.
\autoref{assu:4}(i) basically allows us to construct a subset of
$\mathcal{X}_{e}$ such that $f_z\left(z\left(\mathbf{X}_{e},\boldsymbol{\beta}\right)|\boldsymbol{\beta}\right)$
is uniformly lowered bounded from zero over such subset.
\end{remark}
The following lemma will be useful in the proof of our theorem.
\begin{lemma}
\label{lem:3.1}Suppose that \autoref{assu1}, \autoref{assump:2}(i)-(iii),
\autoref{assu:3}, and \autoref{assu:4} hold with $\upsilon_{G}=3$,
$\upsilon_{K}=2$, and $\upsilon_{f}=3$. Define $\psi\left(n,p,h\right)=h^{-1}\sqrt{\log\left(pnh^{-1}\right)/n}+h^{2}.$
If $h_{n}\rightarrow0$ and $p^{\frac{5p+1}{2\left(p+1\right)}}\psi^{\frac{1}{p+1}}\left(n,p,h_{n}\right)\rightarrow0$
further hold, we have that
\begin{align*}
\sup_{\boldsymbol{\beta}\in\mathcal{B}}\left\Vert \frac{1}{n}\sum_{i=1}^{n}\widehat{G}\left(\left.z\left(\boldsymbol{\mathbf{X}}_{e,i},\boldsymbol{\beta}\right)\right|\boldsymbol{\beta}\right)\boldsymbol{\mathbf{X}}_{i}-\mathbb{E}\left[L\left(z\left(\boldsymbol{\mathbf{X}}_{e,i},\boldsymbol{\beta}\right),\boldsymbol{\beta}\right)\boldsymbol{\mathbf{X}}_{i}\right]\right\Vert = & O_{p}\left(p^{\frac{5p+1}{2\left(p+1\right)}}\psi^{\frac{1}{p+1}}\left(n,p,h_{n}\right)\right).
\end{align*}
\end{lemma}
\begin{proof}[Proof of \autoref{lem:3.1}]
See \autoref{appendixA}.
\end{proof}
\label{rem3.2} \autoref{lem:3.1} 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 \autoref{lem:3.1},
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
\autoref{lem:3.1} 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
\autoref{lem:3.1}, we have that under all conditions as imposed in \autoref{lem:3.1}, there holds
\begin{align}\label{kernel_representation}
\boldsymbol{\beta}_{k+1} & =\boldsymbol{\beta}_{k}-\delta_{k}\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]+\delta_{k}\cdot\left(\text{small order terms}\right).
\end{align}
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
\begin{align}\label{kernel_representation_mean_value}
& \left\{ \int_{\mathcal{X}}\left[G\left(z_{i}^{\star}+\mathbf{X}_{i}^{\mathrm{T}}\Delta\boldsymbol{\beta}_{k}-\mathbf{X}^{\mathrm{T}}\Delta\boldsymbol{\beta}_{k}\right)-G\left(z_{i}^{\star}\right)\right]f_{\mathbf{X}|z}\left(\left.\mathbf{X}\right|z_{i,k},\boldsymbol{\beta}_{k}\right)d\mathbf{X}\right\} \cdot\mathbf{X}_{i} \nonumber \\
& =\int_{0}^{1}\int_{\mathcal{X}}\left[G^{\prime}\left(z_{i}^{\star}+t\left(\mathbf{X}_{i}-\mathbf{X}\right)^{\mathrm{T}}\Delta\boldsymbol{\beta}_{k}\right)f_{\mathbf{X}|z}\left(\left.\mathbf{X}\right|z_{i,k},\boldsymbol{\beta}_{k}\right)\left(\mathbf{X}_{i}\boldsymbol{\mathbf{X}}_{i}^{\mathrm{T}}-\mathbf{X}_{i}\mathbf{X}^{\mathrm{T}}\right)\right]\Delta\boldsymbol{\beta}_{k}d\mathbf{X}dt,
\end{align}
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.
\begin{assumption}
\label{assu:5}There hold
\[
\sup_{\boldsymbol{\beta}\in\mathcal{B}}\overline{\lambda}\left(\varLambda\left(\boldsymbol{\beta}\right)+\varLambda^{\mathrm{T}}\left(\boldsymbol{\beta}\right)\right)\leq\overline{\lambda}_{\varLambda}<\infty,
\]
and
\[
\inf_{\boldsymbol{\beta}\in\mathcal{B}}\underline{\lambda}\left(\varLambda\left(\boldsymbol{\beta}\right)+\varLambda^{\mathrm{T}}\left(\boldsymbol{\beta}\right)\right)\geq\underline{\lambda}_{\varLambda}>0.
\]
\end{assumption}
Based on the above assumptions, we have the following result.
\begin{theorem}
\label{thm:3.1}Suppose that \autoref{assu1}, \autoref{assump:2}(i)--(iii), \autoref{assu:3}--\autoref{assu:5} hold with $\upsilon_{G}=3$, $\upsilon_{K}=2$,
and $\upsilon_{f}=3$, $\delta_{k}=\delta$ such that $\delta<\min\left\{ 1/\left(2\underline{\lambda}_{\varLambda}\right),1/\left(4p^{2}\left\Vert G^{\prime}\right\Vert _{\infty}\right)\right\} $, and that $\boldsymbol{\beta}$ is updated under (\ref{kernel_estimator})
and (\ref{update}).
Define
\[
k_{1,n}^{KBGD}=\frac{\log\left(\left\Vert \Delta\boldsymbol{\beta}_{1}\right\Vert \right)-\log\left(p^{\frac{5p+1}{2\left(p+1\right)}}\psi^{\frac{1}{p+1}}\left(n,p,h_{n}\right)\right)}{-\log\left(1-\delta\underline{\lambda}_{\varLambda}/4\right)}.
\]
Then if $h_{n}\rightarrow0$ and $p^{\frac{5p+1}{2\left(p+1\right)}}\psi^{\frac{1}{p+1}}\left(n,p,h_{n}\right)\rightarrow0$
hold, we have that
\[
\sup_{k\geq k_{1,n}^{KBGD}+1}\left\Vert \Delta\boldsymbol{\beta}_{k}\right\Vert =O_{p}\left(p^{\frac{5p+1}{2\left(p+1\right)}}\psi^{\frac{1}{p+1}}\left(n,p,h_{n}\right)\right).
\]
In particular, if $h_{n}$ is chosen such that $h_{n}=\left(\left(\log n\right)/n\right)^{1/6}$,
then
\[
\sup_{k\geq k_{1,n}^{KBGD}+1}\left\Vert \Delta\boldsymbol{\beta}_{k}\right\Vert =O_{p}\left(p^{\frac{5p+1}{2\left(p+1\right)}}\left(\frac{\log n}{n}\right)^{\frac{1}{3p+3}}\right).
\]
\end{theorem}
\begin{proof}[Proof of \autoref{thm:3.1}] See \autoref{appendixB}.
\end{proof}
\autoref{thm:3.1} implies that the iterative estimator based
on (\ref{kernel_estimator}) and (\ref{update}) 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
\autoref{rem3.2}.
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 \autoref{thm:3.1}, even in the
case of fixed dimensionality, the guaranteed convergence rate of the KBGD estimator
based on (\ref{kernel_estimator}) and (\ref{update}) 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 \citet{ichimura1993semiparametric}. In particular, the algorithm is modified as,
\begin{equation}
\boldsymbol{\beta}_{k+1}=\boldsymbol{\beta}_{k}-\frac{\delta_{k}}{n}\sum_{i=1}^{n}I_{i}^{\phi}\cdot\left(\widehat{G}\left(\left.z_{i,k}\right|\boldsymbol{\beta}_{k}\right)-y_{i}\right)\boldsymbol{\mathbf{X}}_{i},\label{eq:truncation_update}
\end{equation}
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{kernel_estimator}), $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
\begin{equation}
\mathcal{X}_{e}^{\phi}=\left\{ \boldsymbol{\mathbf{X}}_{e}\in\mathcal{X}_{e}:\left|X_{j}\right|\leq1-\phi,0\leq j\leq p\right\} \label{eq:truncation_seet}
\end{equation}
for some $\phi>0$ whose value will be determined later. Different
from (\ref{update}), the update of $\boldsymbol{\beta}_{k}$ based
on (\ref{eq:truncation_update}) 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{eq:truncation_seet}) is that, as we
show in the \autoref{appendixA}, 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
\autoref{thm:3.1}.
\begin{theorem}
\label{thm:3.2}Suppose that all the assumptions and conditions
on $\upsilon_{G}$, $\upsilon_{K}$, and $\upsilon_{f}$ in
\autoref{thm:3.1} hold. Suppose moreover that $h_{n}\rightarrow0$, $\delta_{k}=\delta<\min\left\{ 1/\left(2\underline{\lambda}_{\varLambda}\right),1/\left(4p^{2}\left\Vert G^{\prime}\right\Vert _{\infty}\right)\right\} $,
$\phi<\delta\underline{\lambda}_{\varLambda}/\left(16p^{2}\left\Vert G^{\prime}\right\Vert _{\infty}\zeta\right)$, and that $\boldsymbol{\beta}$ is updated under (\ref{kernel_estimator})
and (\ref{eq:truncation_update}). Define
\[
\widetilde{k}_{1,n}^{KBGD}=\frac{\log\left(\left\Vert \Delta\boldsymbol{\beta}_{1}\right\Vert \right)-\log\left(\psi\left(n,p,h_{n}\right)\right)}{-\log\left(1-\delta\underline{\lambda}_{\varLambda}/8\right)},
\]
then there holds
\[
\sup_{k\geq\widetilde{k}_{1,n}^{KBGD}+1}\left\Vert \Delta\boldsymbol{\beta}_{k}\right\Vert =O_{p}\left(\psi\left(n,p,h_{n}\right)\right).
\]
\end{theorem}
\begin{proof}[Proof of \autoref{thm:3.2}]
See \autoref{appendixB}.
\end{proof}
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 \autoref{thm:3.1}(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
\begin{align}
\Delta\boldsymbol{\beta}_{k+1} & =\Delta\boldsymbol{\beta}_{k}-\frac{\delta_{k}}{n}\sum_{i=1}^{n}\left(\widehat{G}\left(\left.z_{i,k}\right|\boldsymbol{\beta}_{k}\right)-y_{i}\right)\mathbf{X}_{i}^{\phi},\nonumber \\
& =\Delta\boldsymbol{\beta}_{k}-\frac{\delta_{k}}{n}\sum_{i=1}^{n}\left(\widehat{G}\left(\left.z_{i,k}\right|\boldsymbol{\beta}_{k}\right)-\widehat{G}\left(\left.z_{i}^{\star}\right|\boldsymbol{\beta}^{\star}\right)\right)\mathbf{X}_{i}^{\phi}-\delta_{k}\boldsymbol{\xi}_{n}^{\phi}\nonumber \\
& =\int_{0}^{1}\left\{ I_{p}-\frac{\delta_{k}}{n}\sum_{i=1}^{n}\left[\mathbf{X}_{i}^{\phi}\left.\frac{\partial\widehat{G}\left(\left.z\left(\boldsymbol{X}_{e,i},\boldsymbol{\beta}\right)\right|\boldsymbol{\beta}\right)}{\partial\boldsymbol{\beta}^{\mathrm{T}}}\right|_{\boldsymbol{\beta}=\boldsymbol{\beta}^{\star}+t\Delta\boldsymbol{\beta}_{k}}\right]\right\} dt\Delta\boldsymbol{\beta}_{k}-\delta_{k}\boldsymbol{\xi}_{n}^{\phi},\label{eq:gradient_update}
\end{align}
where the integration is understood to be element-wise. To understand
the properties of the above algorithm, we need the following lemmas.
\begin{lemma}
\label{lem3.2}Suppose that all the assumptions in \autoref{thm:3.1}
hold with $\upsilon_{G}=4$, $\upsilon_{K}=3$, and $\upsilon_{f}=4$.
For any sequence of subset $\left\{ \mathcal{B}_{n}\right\} _{n=1}^{\infty}$
with $\mathcal{B}_{n}\subseteq\mathcal{B}$, we have that
\[
\sup_{\boldsymbol{\beta}\in\mathcal{B}_{n}}\left\Vert \frac{1}{n}\sum_{i=1}^{n}\mathbf{X}_{i}^{\phi}\frac{\partial\widehat{G}\left(\left.z\left(\mathbf{X}_{e,i},\boldsymbol{\beta}\right)\right|\boldsymbol{\beta}\right)}{\partial\boldsymbol{\beta}^{\mathrm{T}}}-\varLambda_{\phi}\left(\boldsymbol{\beta}\right)\right\Vert =O_{p}\left(h_{n}^{-2}\sqrt{\left(\log\left(nh_{n}^{-1}\right)\right)/n}+h_{n}^{3}+\sup_{\boldsymbol{\beta}\in\mathcal{B}_{n}}\left\Vert \Delta\boldsymbol{\beta}\right\Vert \right).
\]
\end{lemma}
\begin{proof}[Proof of \autoref{lem3.2}]
See \autoref{appendixA}.
\end{proof}
\begin{lemma}
\label{lem3.3}Suppose that all the assumptions in \autoref{thm:3.1}
hold with $\upsilon_{G}=4$, $\upsilon_{K}=3$, and $\upsilon_{f}=4$.
If $h_{n}$ is chosen such that $h_{n}^{6}n\rightarrow0$, we have
that $\sqrt{n}\boldsymbol{\xi}_{n}^{\phi}\rightarrow_{d}N\left(0,\Sigma_{\boldsymbol{\xi}}^{\phi}\right)$,
where
\[
\Sigma_{\boldsymbol{\xi}}^{\phi}=\mathbb{E}\left[\left(1-G\left(z_{i}^{\star}\right)\right)G\left(z_{i}^{\star}\right)\left(\mathbf{X}_{i}^{\phi}-\mathbb{E}\left(\left.\mathbf{X}_{i}^{\phi}\right|z_{i}^{\star}\right)\right)\left(\mathbf{X}_{i}^{\phi}-\mathbb{E}\left(\left.\mathbf{X}_{i}^{\phi}\right|z_{i}^{\star}\right)\right)^{\mathrm{T}}\right].
\]
\end{lemma}
\begin{proof}[Proof of \autoref{lem3.3}]
See \autoref{appendixA}.
\end{proof}
Now we are in a position to illustrate the results of the asymptotic
normality of our KBGD estimator.
\begin{theorem}
\label{thm:3.3}Suppose that all the assumptions in \autoref{thm:3.1}
hold with $\upsilon_{G}=4$, $\upsilon_{K}=3$, and $\upsilon_{f}=4$.
Suppose moreover that
$\delta_{k}=\delta<\min\left\{ 1/\left(2\underline{\lambda}_{\varLambda}\right),1/\left(4p^{2}\left\Vert G^{\prime}\right\Vert _{\infty}\right)\right\} $, $\phi<\delta \underline{\lambda}_{\varLambda}/\left(16p^{2}\left\Vert G^{\prime}\right\Vert _{\infty}\zeta\right)$,
$h_{n}$ is chosen such that $nh_{n}^{6}\rightarrow0$ and $h_{n}^{4}n/\left(\log n\right)^{2}\rightarrow\infty$,
and that $\boldsymbol{\beta}$ is updated under (\ref{kernel_estimator})
and (\ref{eq:truncation_update}). Then
(i) There holds
\[
\sup_{k\geq\widetilde{k}_{1,n}^{KBGD}+k_{2,n}^{KBGD}+1}\left\Vert \Delta\boldsymbol{\beta}_{k}\right\Vert =O_{p}\left(n^{-1/2}\right),
\]
where $k_{2,n}^{KBGD}$ is given by
\[
k_{2,n}^{KBGD}=\frac{\log\left(n^{1/2}\right)+\log\left(\psi\left(n,p,h_{n}\right)\right)}{-\log\left(1-\delta\underline{\lambda}_{\varLambda}/16\right)};
\]
(ii) Define $\widehat{\boldsymbol{\beta}}=\widehat{\boldsymbol{\beta}}_{k}$
for any $k -\widetilde{k}_{1,n}^{KBGD}-k_{2,n}^{KGBD} \rightarrow\infty$, we have that
\[
\sqrt{n}\left(\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}^{\star}\right)\rightarrow N\left(0,\Sigma_{\boldsymbol{\beta}}^{\phi}\right),
\]
where $\Sigma_{\boldsymbol{\beta}}^{\phi}=\varLambda_{\phi}^{-1}\left(\boldsymbol{\beta}^{\star}\right)\Sigma_{\boldsymbol{\xi}}^{\phi}\left(\varLambda_{\phi}^{-1}\left(\boldsymbol{\beta}^{\star}\right)\right)^{\mathrm{T}}$.
\end{theorem}
\begin{proof}[Proof of \autoref{thm:3.3}]
See \autoref{appendixB}.
\end{proof}
We introduce the estimator for the variance matrix, based on which the confidence interval of $\boldsymbol{\beta}^{\star}$ can be then constructed.
\begin{theorem}\label{thm:3.4}
Suppose that all the assumptions and conditions in \autoref{thm:3.3}
hold. Suppose also that $\widehat{\boldsymbol{\beta}}$ is defined
as in \autoref{thm:3.3}. Define
\[
\widehat{\Sigma}_{\boldsymbol{\xi}}^{\phi}=\frac{1}{n}\sum_{i=1}^{n}\left(\widehat{G}_{i}\left(1-\widehat{G}_{i}\right)\left(\mathbf{X}_{i}^{\phi}-\widehat{\mathbb{E}}\left(\left.\mathbf{X}_{i}^{\phi}\right|\widehat{z}_{i}\right)\right)\left(\mathbf{X}_{i}^{\phi}-\widehat{\mathbb{E}}\left(\left.\mathbf{X}_{i}^{\phi}\right|\widehat{z}_{i}\right)\right)^{\mathrm{T}}\right),
\]
and
\[
\widehat{\varLambda}_{\phi}\left(\widehat{\boldsymbol{\beta}}\right)=\frac{1}{n}\sum_{i=1}^{n}\mathbf{X}_{i}^{\phi}\frac{\partial\widehat{G}\left(\left.z\left(\mathbf{X}_{e,i},\widehat{\boldsymbol{\beta}}\right)\right|\widehat{\boldsymbol{\beta}}\right)}{\partial\boldsymbol{\beta}^{\mathrm{T}}},
\]
where
\[
\widehat{G}_{i}=\frac{\sum_{j=1}^{n}K_{h_{n}}\left(\widehat{z}_{i}-\widehat{z}_{j}\right)y_{j}}{\sum_{j=1}^{n}K_{h_{n}}\left(\widehat{z}_{i}-\widehat{z}_{j}\right)},\ \widehat{\mathbb{E}}\left(\left.\mathbf{X}_{i}^{\phi}\right|\widehat{z}_{i}\right)=\frac{\sum_{j=1}^{n}K_{h_{n}}\left(\widehat{z}_{i}-\widehat{z}_{j}\right)\mathbf{X}_{j}^{\phi}}{\sum_{j=1}^{n}K_{h_{n}}\left(\widehat{z}_{i}-\widehat{z}_{j}\right)},
\]
and $\widehat{z}_{i}=X_{0,i}+\mathbf{X}_{i}^{\mathrm{T}}\widehat{\boldsymbol{\beta}}$.
Then we have that
\[
\left\Vert \widehat{\varLambda}_{\phi}^{-1}\left(\widehat{\boldsymbol{\beta}}\right)\widehat{\Sigma}_{\boldsymbol{\xi}}^{\phi}\left(\widehat{\varLambda}_{\phi}^{-1}\left(\widehat{\boldsymbol{\beta}}\right)\right)^{\mathrm{T}}-\Sigma_{\boldsymbol{\beta}}^{\phi}\right\Vert \rightarrow_{p}0.
\]
\end{theorem}
\begin{proof}[Proof of \autoref{thm:3.4}]
See \autoref{appendixB}.
\end{proof}
We finally provide some remarks for the KBGD estimators.
\begin{remark}
\label{rem4}We first provide some remarks on the implementation
of our KBGD estimator. The KBGD estimator might be sensitive to the
data magnitude. So when implementing such an estimator, we recommend
first standardizing the data so that each covariate has zero mean
and unit variance. Note that when constructing the KBGD estimator,
we normalize the coefficient of $X_{0,i}$ to 1, indicating that the
coefficients of $\mathbf{X}_{e,i}$ can not all be zeros. So we need
to test whether at least one covariate affects the conditional probability
of $y_{i}=1$. One option is to run a Logit or Probit regression and
test whether all the coefficients are equal to zero.
When applying our algorithm, it is also crucial to determine the learning rate $\delta$, bandwidth of kernel estimator $h_{n}$, and terminating
conditions of the algorithm. In \autoref{thm:3.3}, the tuning
parameter $\delta$ is required to be smaller than $1/\left(2\underline{\lambda}_{\varLambda}\right)$
and $1/\left(4p^{2}\left\Vert G^{\prime}\right\Vert _{\infty}\right)$,
neither of which is known. So we recommend setting $\delta$ to be
1 in the first place, and gradually shrink it if the iteration does
not converge. For the choice of the bandwidth $h_{n}$, \autoref{thm:3.3}
requires that $h_{n}$ is chosen such that $nh_{n}^{6}\rightarrow0$
and $nh_{n}^{4}/\left(\log n\right)^{2}\rightarrow\infty$. As a rule
of thumb, we recommend choosing $h_{n}=C\cdot n^{-1/5}$. For the
choice of the constant $C$, we can choose $C=C_{k}=\text{std}\left(z_{i,k}\right)$
for the $k$-th round of iteration and $C=\text{std}\left(\widehat{z}_{i}\right)$
when estimating the variance $\Sigma_{\boldsymbol{\beta}}^{\phi}$.
We finally discuss the terminating conditions. As we show in
\autoref{thm:3.3}, to obtain root-$n$ consistency and asymptotic normality,
the iteration number is required to be only of order $\log\left(n\right)$.
However, such rule can not be directly applied to determine the number
of iterations since the initial distance $\left\Vert \Delta\boldsymbol{\beta}_{1}\right\Vert $
as well as the lower bounded on the eigenvalues $\underline{\lambda}_{\varLambda}$
are both unknown. We recommend the terminating condition
$\max_{1\leq j\leq p}|\widehat{\beta}_{j,k+1}-\widehat{\beta}_{j,k}|<\varrho$
for some predetermined tolerance $\varrho$. During the simulation,
we choose $\varrho=10^{-5}$. Note that in many cases, $\max_{1\leq j\leq p}|\widehat{\beta}_{j,k+1}-\widehat{\beta}_{j,k}|$
may not be monotonically decreasing with $k$; in some extreme cases,
$\max_{1\leq j\leq p}|\widehat{\beta}_{j,k+1}-\widehat{\beta}_{j,k}|$
may even be oscillating and does not shrink to zero. On these condition,
we recommend decreasing $\delta$ or choosing $h_{n}=C\cdot n^{-1/5}$
with $C=1$ when iterating. If the maximum distance still oscillates,
we recommend stop iteration when the maximum distance achieves its
minimum value.
\end{remark}
\begin{remark}
\label{rem5}Our previous discussion has be confined to the case where
all the covariates are continuously distributed, while our algorithm
can be directly applied to the case where there are discrete covariates
without any modifications. The basic reason is that, in contrast to the
average derivative approach \citep{stoker1986consistent,powell1989semiparametric} that uses the differentiation with respect to covariates,
the KBGD estimator performs differentiation with respect to the parameters,
so it does not impose requirements on the continuity of the covariates.
It should be noted that we do require at least one continuous covariate
to guarantee identification of the parameters. For simplicity, we
recommend choosing a continuous covariate as the standardization covariate
$X_{0}$. Finally, we point out that stronger assumption should be
imposed to make our results valid when there are discrete covariates.
In particular, suppose that $\mathbf{X}_{e}=\left(\mathbf{X}_{c}^{\mathrm{T}},\mathbf{X}_{d}^{\mathrm{T}}\right)^{\mathrm{T}}$,
where $\mathbf{X}_{c}$ is the collection of all the continuous covariates,
whereas $\mathbf{X}_{d}$ is the collection of all the discrete covariates.
Also denote the density function of $\mathbf{X}_{c}$ conditional
on $\mathbf{X}_{d}$ as $f_{\mathbf{X}_{c}|\mathbf{X}_{d}}\left(\mathbf{X}_{c}|\mathbf{X}_{d}\right)$.
Then we require that all the conditions imposed on the $f_{e}\left(\mathbf{X}_{e}\right)$
hold for $f_{\mathbf{X}_{c}|\mathbf{X}_{d}}\left(\mathbf{X}_{c}|\mathbf{X}_{d}\right)$
for any realizations of $\mathbf{X}_{d}$.
\end{remark}
\subsection{\label{section4} 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
\begin{equation}
\widehat{\boldsymbol{\pi}}_{q,n,k}=\left(\sum_{i=1}^{n}\boldsymbol{r}_{q}\left(z_{i,k}\right)\boldsymbol{r}_{q}^{\mathrm{T}}\left(z_{i,k}\right)\right)^{-1}\left(\sum_{i=1}^{n}\boldsymbol{r}_{q}\left(z_{i,k}\right)y_{i}\right).\label{OLS_update}
\end{equation}
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
\begin{equation}
\widehat{G}\left(\left.z\right|\boldsymbol{\beta}_{k}\right)=\boldsymbol{r}_{q}^{\mathrm{T}}\left(z\right)\widehat{\boldsymbol{\pi}}_{n,q,k},\ -\infty<z<\infty.\label{sivev_estimator}
\end{equation}
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{update}).
We iterate sequentially based on (\ref{OLS_update}), (\ref{sivev_estimator})
and (\ref{update}) until some terminating conditions are satisfied.
The resulting estimator is then labeled as the \textit{sieve-based batch gradient
descent estimator} (SBGD estimator).
\begin{remark}\label{rem6}
In the above SBGD procedure, we update the sieve parameter based on the
OLS-type estimation. An alternative procedure can be based on the flexible
Logit regression proposed by \citet{hirano2003efficient}. The advantage
of using flexible Logit regression is that the estimated CDF $\widehat{G}\left(\left.z\right|\boldsymbol{\beta}_{k}\right)$
always falls between 0 and 1 for all $z$, which makes the update more stable. While
the disadvantage of such update is that the flexible Logit regression
is based on MLE, which does not allow for an analytical solution. Using
numerical optimization to solve for the sieve coefficients in each
round of update will add to additional computational burdens.
\end{remark}
\begin{remark}\label{rem7}
Compared with the KBGD algorithm, the SBGD procedure has at least
two advantages. On the one side, the sieve-based approximation for
the unknown CDF is global and guarantees uniform approximation error
rate. This allows us to update the parameter without performing any
form of trimming as we did for the KBGD estimator. Moreover, this allows us to develop the asymptotic distribution of the SBGD estimator for the case of increasing dimensionality. On the otherhand,
the KBGD procedure relies on the kernel estimation of CDF $G$ at
$n$ data points, whose computational complexity of each update is
of order $O\left(n^{2}\right)$. While the most time-consuming part
of the SBGD procedure is the OLS procedure (\ref{OLS_update}), whose
computational complexity is of order $O\left(nq^{2}+q^{3}\right)$.
When $q/\sqrt{n}\rightarrow0$, the computational burden of SBGD estimator
will be substantially lower than that of KBGD estimator.
\end{remark}
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,
\begin{align}
\boldsymbol{\beta}_{k+1} & =\boldsymbol{\beta}_{k}-\frac{\delta_{k}}{n}\sum_{i=1}^{n}\left(\mathbf{X}_{i}-\mathfrak{X}_{q,n}\left(z_{i,k},\boldsymbol{\beta}_{k}\right)\right)\left(G\left(z_{i,k}\right)-G\left(z_{i}^{\star}\right)\right)\nonumber \\
& -\frac{\delta_{k}}{n}\sum_{i=1}^{n}\mathbf{X}_{i}\boldsymbol{r}_{q}^{\mathrm{T}}\left(z_{i,k}\right)\Gamma_{q,n,k}^{-1}\left(\frac{1}{n}\sum_{j=1}^{n}\boldsymbol{r}_{q}\left(z_{j,k}\right)R_{q}\left(z_{j,k}\right)+\frac{1}{n}\sum_{i=1}^{n}\boldsymbol{r}_{q}\left(z_{j,k}\right)\varepsilon_{j}\right)\nonumber \\
& +\frac{\delta_{k}}{n}\sum_{i=1}^{n}\left(R_{q}\left(z_{i,k}\right)\mathbf{X}_{i}+\varepsilon_{i}\mathbf{X}_{i}\right),\label{SBGD_expression}
\end{align}
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.
\begin{assumption}
\label{assu:6}(i) There holds $\max_{0\leq j\leq q}\left\Vert r_{j}\right\Vert _{\infty}\leq D_{q,0}$,
$\max_{0\leq j\leq q}\left\Vert r_{j}^{\prime}\right\Vert _{\infty}\leq D_{q,1}$,
and $\max_{0\leq j\leq q}\left\Vert r_{j}^{\prime\prime}\right\Vert _{\infty}\leq D_{q,2}$;
(ii) Define $\Gamma_{q}\left(\boldsymbol{\beta}\right)=\mathbb{E}\left(\boldsymbol{r}_{q}\left(X_{0}+\mathbf{X}^{\mathrm{T}}\boldsymbol{\beta}\right)\boldsymbol{r}_{q}^{\mathrm{T}}\left(X_{0}+\mathbf{X}^{\mathrm{T}}\boldsymbol{\beta}\right)\right),$
there hold $\inf_{\boldsymbol{\beta}\in\mathcal{B}}\underline{\lambda}\left(\Gamma_{q}\left(\boldsymbol{\beta}\right)\right)\geq\underline{\lambda}_{\Gamma}>0$
and $\sup_{\boldsymbol{\beta}\in\mathcal{B}}\overline{\lambda}\left(\Gamma_{q}\left(\boldsymbol{\beta}\right)\right)\leq\overline{\lambda}_{\Gamma}<\infty$
for all $q$; (iii) There hold $\sup_{z\in R} \left|G\left(z\right) - \boldsymbol{r}^{\mathrm{T}}\left(z\right)\boldsymbol{\pi}_q^{\star}\right| \leq \mathcal{E}_{q,0}$ and $\sup_{z\in R} \left|G^{\prime}\left(z\right) - \left(\boldsymbol{r}^{\prime }\left(z\right)\right)^\mathrm{T}\boldsymbol{\pi}_q^{\star}\right| \leq \mathcal{E}_{q,1}$, where $\boldsymbol{r}^{\prime}(z) = \left(r^{\prime}_0(z), \cdots, r^{\prime}_q(z)\right)^{\mathrm{T}}$.
\end{assumption}
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.
\begin{lemma}
\label{lem4.1}Define
$
\chi_{1,n}=\sqrt{pq^{2}D_{q,0}^{4}\log\left(pqD_{q,0}D_{q,1}n\right)/n},
$
and
$
\chi_{2,n}= \sqrt{p}qD_{q,0}^{2}\left(\chi_{1,n}+ \mathcal{E}_{q,0}\right).
$
Suppose that \autoref{assu1}, \autoref{assump:2}(i)-(iii),
and \autoref{assu:6} hold, and moreover, $\upsilon_G\geq 1$ and the combination of $p$, $q$ and $\upsilon_{G}$ guarantees that $\chi_{1,n}\rightarrow0$
as $n\rightarrow\infty$. Then the following holds,
\[
\boldsymbol{\beta}_{k+1}=\boldsymbol{\beta}_{k}-\delta_{k}\mathbb{E}\left[\left(\mathbf{X}-\mathfrak{X}_{q}\left(z\left(\mathbf{X}_{e},\boldsymbol{\beta}_{k}\right),\boldsymbol{\beta}_{k}\right)\right)\left(G\left(z\left(\mathbf{X}_{e},\boldsymbol{\beta}_{k}\right)\right)-G\left(z\left(\mathbf{X}_{e},\boldsymbol{\beta}^{\star}\right)\right)\right)\right]+\delta_{k}\mathfrak{R}_{n,k},
\]
where $\sup_{k\geq1}\left\Vert \mathfrak{R}_{n,k}\right\Vert =O_{p}\left(\chi_{2,n}\right)$.
\end{lemma}
\begin{proof}[Proof of \autoref{lem4.1}]
See \autoref{appendixA}.
\end{proof}
Obviously, \autoref{lem4.1} provides a parallel result to (\ref{kernel_representation}). 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 \autoref{lem4.1}, we have that
\begin{align}
\Delta\boldsymbol{\beta}_{k+1} & =\left\{ \int_{0}^{1}\left(I_{p}-\delta_{k}\Psi_{q}\left(t,\boldsymbol{\beta}_{k}\right)\right)dt\right\} \Delta\boldsymbol{\beta}_{k}+\delta_{k}\mathfrak{R}_{n,k}.\label{contraction_sieve}
\end{align}
Obviously, (\ref{contraction_sieve}) is also a parallel result to (\ref{kernel_representation_mean_value}). As a result, to ensure that (\ref{contraction_sieve})
actually constitutes a contraction for $\left\Vert \Delta\boldsymbol{\beta}_{k}\right\Vert $,
we impose the following assumption that is similar to \autoref{assu:5}.
\begin{assumption}
\label{assu:7}For any $q\geq0$, there hold
\[
\inf_{0\leq t\leq1,\boldsymbol{\beta}\in\mathcal{B}}\underline{\lambda}\left(\Psi_{q}\left(t,\boldsymbol{\beta}\right)+\Psi_{q}^{\mathrm{T}}\left(t,\boldsymbol{\beta}\right)\right)\geq\underline{\lambda}_{\Psi}>0,
\]
\[
\sup_{0\leq t\leq1,\boldsymbol{\beta}\in\mathcal{B}}\underline{\lambda}\left(\Psi_{q}\left(t,\boldsymbol{\beta}\right)+\Psi_{q}^{\mathrm{T}}\left(t,\boldsymbol{\beta}\right)\right)\geq\overline{\lambda}_{\Psi}<\infty.
\]
\end{assumption}
Based on the above assumptions, we have the following result.
\begin{theorem}
\label{thm4.1}Suppose that \autoref{assu1}, \autoref{assump:2}(i)-(iii),
\autoref{assu:6} and \autoref{assu:7} hold, $\upsilon_G\geq 1$, and the combination of $p$, $q$ and $\upsilon_{G}$ guarantees that
$\chi_{1,n}\rightarrow0$ as $n\rightarrow \infty$. Suppose moreover that the learning rate is chosen such that $\delta_{k}=\delta$
with
$
0<\delta<\min\left\{ 1/\left(2\underline{\lambda}_{\Psi}\right),\underline{\lambda}_{\Psi}/\left(2\left\Vert G^{\prime}\right\Vert _{\infty}^{2}p^{2}\left\{ 1+\underline{\lambda}_{\Gamma}^{-1}qD_{q,0}^{2}\right\} ^{2}\right)\right\}
$, and that $\boldsymbol{\beta}$ is updated based on (\ref{OLS_update}), (\ref{sivev_estimator})
and (\ref{update}).
Define
\[
k_{1,n}^{SBGD}=\frac{\log\left(\left\Vert \Delta\boldsymbol{\beta}_{1}\right\Vert \right)-\log\left(\chi_{2,n}\right)}{-\log\left(1-\underline{\lambda}_{\Psi}\delta/4\right)},
\]
then we have that
\[
\sup_{k\geq k_{1,n}^{SBGD}+1}\left\Vert \Delta\boldsymbol{\beta}_{k}\right\Vert =O_{p}\left(\chi_{2,n}\right).
\]
\end{theorem}
\begin{proof}[Proof of \autoref{thm4.1}]
See \autoref{appendixB}.
\end{proof}
According to \autoref{thm4.1}, 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{SBGD_expression}), we have that
\begin{align*}
\Delta\boldsymbol{\beta}_{k+1} & = \left\{I_p-\delta_{k} \int^{1}_0\frac{1}{n}\sum_{i=1}^{n}G^{\prime}\left(z_i^{\star} + t\mathbf{X}_i^{\mathrm{T}}\Delta\boldsymbol{\beta}_k\right)\left(\mathbf{X}_{i}\mathbf{X}_{i}^{\mathrm{T}}-\mathfrak{X}_{q,n}\left(z_{i,k},\boldsymbol{\beta}_{k}\right)\mathbf{X}_{i}^{\mathrm{T}}\right)dt\right\}\Delta\boldsymbol{\beta}_{k} \\
& -\frac{\delta_{k}}{n}\sum_{i=1}^{n}\mathbf{X}_{i}\boldsymbol{r}_{q}^{\mathrm{T}}\left(z_{i,k}\right)\Gamma_{q,n,k}^{-1}\left(\frac{1}{n}\sum_{j=1}^{n}\boldsymbol{r}_{q}\left(z_{j,k}\right)R_{q}\left(z_{j,k}\right)+\frac{1}{n}\sum_{i=1}^{n}\boldsymbol{r}_{q}\left(z_{j,k}\right)\varepsilon_{j}\right) \\
& +\frac{\delta_{k}}{n}\sum_{i=1}^{n}\left(R_{q}\left(z_{i,k}\right)\mathbf{X}_{i}+\varepsilon_{i}\mathbf{X}_{i}\right).
\end{align*}
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 \autoref{lem3.2} and \autoref{lem3.3},
we provide two additional lemmas that are useful to understand the above
algorithm.
\begin{lemma}
\label{lem4.2}Suppose that \autoref{assu1}, \autoref{assump:2}(i)-(iii),
and \autoref{assu:6} hold, $\upsilon_G\geq 2$ and the combination of $p$, $q$ and $\upsilon_{G}$ guarantees that
$\chi_{1,n}\rightarrow0$ as $n\rightarrow \infty$. Then for any sequence $\left\{ \mathcal{B}_{n}\right\} _{n=1}^{\infty}$
with \textup{$\mathcal{B}_{n}\subseteq\mathcal{B}$} we have that
\begin{align*}
& \sup_{0\leq t\leq1,\boldsymbol{\beta}\in\mathcal{B}_{n}}\left\Vert \frac{1}{n}\sum_{i=1}^{n}G^{\prime}\left(z_{i}^{\star}+t\mathbf{X}_{i}^{\mathrm{T}}\Delta\boldsymbol{\beta}\right)\left(\mathbf{X}_{i}\mathbf{X}_{i}^{\mathrm{T}}-\mathfrak{X}_{q,n}\left(z\left(\mathbf{X}_{e,i},\boldsymbol{\beta}\right),\boldsymbol{\beta}\right)\mathbf{X}_{i}^{\mathrm{T}}\right)-\Psi_{q}^{\star}\right\Vert \\
& =O_{p}\left(pqD_{q,0}^{2}\chi_{1,n} + \sqrt{p^3}q^{2}D_{q,0}^{3}D_{q,1}\sup_{\boldsymbol{\beta}\in\mathcal{B}_{n}}\left\Vert \Delta\boldsymbol{\beta}\right\Vert \right).
\end{align*}
\end{lemma}
\begin{proof}[Proof of \autoref{lem4.2}]
See \autoref{appendixA}.
\end{proof}
\begin{lemma}
\label{lem4.3}Suppose that \autoref{assu1}, \autoref{assump:2}(i)-(iii), \autoref{assu:6}, and \autoref{assu:7} hold, and the combination of $p$, $q$ and $\upsilon_{G}$ guarantees that
$\chi_{1,n}\rightarrow0$ as $n\rightarrow \infty$. Define $\boldsymbol{r}_{q,i,k}=\boldsymbol{r}_{q}\left(z_{i,k}\right)$, and
$R_{q,i,k}=R_q\left(z_{i,k}\right)$.
Also define
\[
\chi_{3,n}=\sqrt{p^{2}qD_{q,1}^{2}\log\left(pqD_{q,2}n\right)/n},
\]
then we have that
\begin{align*}
\sup_{k\geq k_{1,n}^{SBGD}+1} & \left\Vert \frac{1}{n}\sum_{i=1}^{n}\mathbf{X}_{i}\boldsymbol{r}_{q,i,k}^{\mathrm{T}}\Gamma_{q,n,k}^{-1}\left(\frac{1}{n}\sum_{j=1}^{n}\boldsymbol{r}_{q,j,k}R_{q,j,k}+\frac{1}{n}\sum_{j=1}^{n}\boldsymbol{r}_{q,j,k}\varepsilon_{j}\right)+\right.\\
&
\left.\frac{1}{n}\sum_{i=1}^{n}R_{q}\left(z_{i,k}\right)\mathbf{X}_{i} -\frac{1}{n}\sum_{i=1}^{n}\mathfrak{X}_{q}\left(z_{i}^{\star},\boldsymbol{\beta}^{\star}\right)\varepsilon_{j}\right\Vert =O_{p}\left(\chi_{4,n}\right),
\end{align*}
where $\chi_{4,n}=\sqrt{p}qD_{q,0}^{2}\mathcal{E}_{q,0}+\sqrt{pq}D_{q,0}\chi_{2,n}\chi_{3,n}+\chi_{2,n}\sqrt{p^2q^4D_{q,0}^{6}D_{q,1}^2\left(\log q\right)/n}$.
\end{lemma}
\begin{proof}[Proof of \autoref{lem4.3}]
See \autoref{appendixA}.
\end{proof}
Based on the above two lemmas, we are now ready to study the asymptotic distribution of the SBGD estimator.
\begin{theorem}
\label{thm4.2}Suppose that \autoref{assu1}, \autoref{assump:2}(i)-(iii),
\autoref{assu:6} and \autoref{assu:7} hold, $\upsilon_G\geq 2$, the combination of $p$, $q$ and $\upsilon_{G}$ guarantees that
$\chi_{1,n}\rightarrow0$ as $n\rightarrow \infty$, and that $\boldsymbol{\beta}$ is updated based on (\ref{OLS_update}), (\ref{sivev_estimator})
and (\ref{update}). We have
that
(i) There holds
\[
\Delta\boldsymbol{\beta}_{k+1}=\left(I_{p}-\delta\Psi_{q}^{\star}\right)\Delta\boldsymbol{\beta}_{k}+\frac{\delta}{n}\sum_{i=1}^{n}\left(\mathbf{X}_{i} - \mathfrak{X}_{q}\left(z_{i}^{\star},\boldsymbol{\beta}^{\star}\right) \right)\varepsilon_{i}+\widetilde{\mathfrak{R}}_{n,k},
\]
where $\sup_{k\geq k_{1,n}^{SBGD}+1}\left\Vert \widetilde{\mathfrak{R}}_{n,k}\right\Vert =O_{p}\left(\chi_{5,n}\right)$
with
\[
\chi_{5,n}=\sqrt{p}qD_{q,0}^{2}\left(p+qD_{q,0}D_{q,1}\right)\chi_{2,n}^{2}+\chi_{4,n};
\]
(ii) Define $\widehat{\boldsymbol{\beta}}=\boldsymbol{\beta}_{k+k_{1,n}^{SBGD}+k_{2,n}^{SBGD}+1}$
with
\[
k_{2,n}^{SBGD}=\frac{-\log\chi_{2,n}+\log\sqrt{n}}{-\log\left(1-\underline{\lambda}_{\Psi}\delta/4\right)},
\]
and any $k\geq1$. If the combination of $p$, $q$ and $\upsilon_{G}$ further guarantees that $\sqrt{n}\chi_{5,n}\rightarrow0$ as $n\rightarrow\infty$,
we have that
\[
\sqrt{n}\left(\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}^{\star}\right)=\Psi_{q}^{\star-1}\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left( \mathbf{X}_{i} - \mathfrak{X}_{q}\left(z_{i}^{\star},\boldsymbol{\beta}^{\star}\right) \right)\varepsilon_{i}+o_{p}\left(n^{-\frac{1}{2}}\right).
\]
Then for any $p\times1$ vector $\rho$ such that $\left\Vert \rho\right\Vert <\infty$
and $\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\rho^{\mathrm{T}}\Psi_{q}^{\star-1}\left(\mathbf{X}_{i} - \mathfrak{X}_{q}\left(z_{i}^{\star},\boldsymbol{\beta}^{\star}\right) \right)\varepsilon_{i}\rightarrow_{d}N\left(0,\sigma_{S}^{2}\left(\rho\right)\right)$
with
\[
\sigma_{S}^{2}\left(\rho\right)=\lim_{n\rightarrow\infty}\rho^{\mathrm{T}}\Psi_{q}^{\star-1}\mathbb{E}\left\{ G\left(z_{i}^{\star}\right)\left(1-G\left(z_{i}^{\star}\right)\right)\left(\mathbf{X}_{i} - \mathfrak{X}_{q}\left(z_{i}^{\star},\boldsymbol{\beta}^{\star}\right) \right)\left( \mathbf{X}_{i} - \mathfrak{X}_{q}\left(z_{i}^{\star},\boldsymbol{\beta}^{\star}\right) \right)^{\mathrm{T}}\right\} \left(\Psi_{q}^{\star-1}\right)^{\mathrm{T}}\rho,
\]
there holds \[\sqrt{n}\rho^{\mathrm{T}}\left(\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}^{\star}\right)\rightarrow_{d}N\left(0,\sigma_{S}^{2}\left(\rho\right)\right)\].
\end{theorem}
\begin{proof}[Proof of \autoref{thm4.2}]
See \autoref{appendixB}.
\end{proof}
We now provide the estimator for the variance.
\begin{theorem}\label{thm4.3}
Suppose that all the conditions listed in \autoref{thm4.2} hold and $pq^{2}D_{q,0}^{4}\mathcal{E}_{q,1}\rightarrow 0$ as $n\rightarrow0$. Let
$\widehat{\boldsymbol{\beta}}$ be as defined as in \autoref{thm4.2}.
Define $\widehat{\boldsymbol{r}}_{q,i}=\boldsymbol{r}_{q}\left(z\left(\mathbf{X}_{e,i},\widehat{\boldsymbol{\beta}}\right)\right)$,
$\widehat{\boldsymbol{r}}_{q,i}^{\prime}=\boldsymbol{r}_{q}^{\prime}\left(z\left(\mathbf{X}_{e,i},\widehat{\boldsymbol{\beta}}\right)\right)$,
$\widehat{\boldsymbol{\pi}}_{q}=\left(\sum_{i=1}^{n}\widehat{\boldsymbol{r}}_{q,i}\widehat{\boldsymbol{r}}_{q,i}^{\mathrm{T}}\right)^{-1}\left(\sum_{i=1}^{n}\widehat{\boldsymbol{r}}_{q,i}y_{i}\right),$
$\widehat{G}_{i}=\widehat{\boldsymbol{r}}_{q,i}^{\mathrm{T}}\widehat{\boldsymbol{\pi}},\ \widehat{G}_{i}^{\prime}=\widehat{\boldsymbol{r}}_{q,i}^{\prime\mathrm{T}}\widehat{\boldsymbol{\pi}}_{q},$
$\widehat{\Psi}_{q,i}^{\star}=\frac{1}{n}\sum_{i=1}^{n}\widehat{G}_{i}^{\prime}\cdot\left(\mathbf{X}_{i}\mathbf{X}_{i}^{\mathrm{T}}-\mathfrak{X}_{q,n}\left(\widehat{z}_{i},\widehat{\boldsymbol{\beta}}\right)\mathbf{X}_{i}^{\mathrm{T}}\right)$, $\widehat{\mathfrak{X}}_{q,i}=\frac{1}{n}\sum_{j=1}^{n}\mathbf{X}_{j}\widehat{\boldsymbol{r}}_{q,j}^{\mathrm{T}}\Gamma_{q,n}^{-1}\left(\widehat{\boldsymbol{\beta}}\right)\widehat{\boldsymbol{r}}_{q,i},$
and
\[
\widehat{\sigma}_{S}^{2}\left(\rho\right)=\rho^{\mathrm{T}}\widehat{\Psi}_{q}^{\star-1}\frac{1}{n}\sum_{i=1}^{n}\left\{ \widehat{G}_{i}\left(1-\widehat{G}_{i}\right)\left(\mathbf{X}_{i} - \widehat{\mathfrak{X}}_{q,i} \right)\left( \mathbf{X}_{i} - \widehat{\mathfrak{X}}_{q,i} \right)^{\mathrm{T}}\right\} \left(\widehat{\Psi}_{q}^{\star-1}\right)^{\mathrm{T}}\rho,
\]
Then for any $p\times1$ vector $\rho$ such that $\left\Vert \rho\right\Vert <\infty$,
there holds
\[
\left|\widehat{\sigma}_{S}^{2}\left(\rho\right)-\sigma_{S}^{2}\left(\rho\right)\right|\rightarrow_{p}0.
\]
\end{theorem}
\begin{proof}[Proof of \autoref{thm4.3}]
See \autoref{appendixB}.
\end{proof}
We finally provide some remarks on the empirical applications of the SBGD estimator.
\begin{remark}\label{rem8}
For the choice of sieve functions, we can use polynomial series for the case where the error term $u_i$ has bounded support and Hermite polynomials for the case where $u_i$ has unbounded support. Note that when using polynomial series $\left\{1, z, z^2, \cdots, z^q\right\}$, the correlation between the sieve functions increases as the approximation order $q$ increases, which may lead to a violation of \autoref{assu:6}(ii). To improve the finite sample performance of our method, we recommend using Chebyshev or Legendre polynomials. Moreover, in the case where $u_i$ has unbounded support, following \citet{bierens2014consistency}, we recommend first conducting the following transformation $G\left(z\right)$ = $\widetilde{G}\left(T\left(z\right)\right)$, where $T:R \mapsto [-1,1]$ is a differentiable function, and then using standard Chebyshev or Legendre polynomials to approximate $\widetilde{G}$. For example, in our following simulations and empirical applications in Section \ref{section6}, we use $T\left(z\right) = 2\pi^{-1}\arctan\left(z\right)$. For the uniform error bound of truncated Legendre polynomials, see \citet{wang2012convergence}.
\end{remark}
\section{\label{section4}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.
\begin{table}
\begin{centering}
\caption{\label{tab1}Finite Sample Performance of KBGD and SBGD Estimators}
\begin{tabular}{>{\centering}p{0.6cm}>{\centering}p{0.9cm}>{\centering}p{0.9cm}>{\centering}p{0.9cm}>{\centering}p{0.9cm}>{\centering}p{0.9cm}>{\centering}p{0.9cm}c>{\centering}p{0.9cm}>{\centering}p{0.9cm}>{\centering}p{0.9cm}>{\centering}p{0.9cm}>{\centering}p{0.9cm}>{\centering}p{0.9cm}}
\hline
& \multicolumn{2}{c}{Bias} & \multicolumn{2}{c}{RMSE} & \multicolumn{2}{c}{CR} & & \multicolumn{2}{c}{Bias} & \multicolumn{2}{c}{RMSE} & \multicolumn{2}{c}{CR}\tabularnewline
\hline
& K{\small{}BGD} & S{\small{}BGD} & K{\small{}BGD} & S{\small{}BGD} & K{\small{}BGD} & S{\small{}BGD} & & K{\small{}BGD} & S{\small{}BGD} & K{\small{}BGD} & S{\small{}BGD} & K{\small{}BGD} & S{\small{}BGD}\tabularnewline
\hline
\multicolumn{7}{c}{$n=2500$} & & \multicolumn{6}{c}{$n=5000$}\tabularnewline
\hline
$\beta_{1}$ & 0.0024 & 0.0031 & 0.1193 & 0.1240 & 0.9600 & 0.9680 & & 0.0047 & 0.0005 & 0.0844 & 0.0867 & 0.9500 & 0.9600\tabularnewline
$\beta_{2}$ & 0.0002 & 0.0055 & 0.1255 & 0.1336 & 0.9480 & 0.9500 & & 0.0031 & 0.0074 & 0.0846 & 0.0878 & 0.9520 & 0.9540\tabularnewline
$\beta_{3}$ & 0.0136 & 0.0260 & 0.1544 & 0.1791 & 0.9480 & 0.9460 & & 0.0004 & 0.0074 & 0.1053 & 0.1112 & 0.9320 & 0.9320\tabularnewline
$\beta_{4}$ & 0.0093 & 0.0213 & 0.1551 & 0.1706 & 0.9500 & 0.9440 & & 0.0012 & 0.0095 & 0.1035 & 0.1117 & 0.9600 & 0.9500\tabularnewline
$\beta_{5}$ & 0.0257 & 0.0482 & 0.2511 & 0.2968 & 0.9540 & 0.9400 & & 0.0007 & 0.0168 & 0.1648 & 0.1889 & 0.9400 & 0.9480\tabularnewline
$\beta_{6}$ & 0.0236 & 0.0477 & 0.2502 & 0.2860 & 0.9480 & 0.9580 & & 0.0121 & 0.0269 & 0.1723 & 0.1931 & 0.9540 & 0.9360\tabularnewline
$\beta_{7}$ & 0.0500 & 0.0964 & 0.4513 & 0.5416 & 0.9640 & 0.9420 & & 0.0051 & 0.0352 & 0.3083 & 0.3525 & 0.9440 & 0.9420\tabularnewline
$\beta_{8}$ & 0.0447 & 0.0920 & 0.4662 & 0.5441 & 0.9360 & 0.9520 & & 0.0098 & 0.0394 & 0.3121 & 0.3477 & 0.9420 & 0.9440\tabularnewline
$\beta_{9}$ & 0.0242 & 0.0454 & 0.2921 & 0.3303 & 0.9480 & 0.9500 & & 0.0072 & 0.0048 & 0.1840 & 0.1909 & 0.9540 & 0.9560\tabularnewline
$\beta_{10}$ & 0.0168 & 0.0338 & 0.1881 & 0.2223 & 0.9520 & 0.9440 & & 0.0030 & 0.0147 & 0.1247 & 0.1402 & 0.9440 & 0.9380\tabularnewline
\hline
\end{tabular}
\par\end{centering}
{\footnotesize{}NOTE: For KBGD estimator, we use fourth-order Epanechinikov
kernel to construct the Nadaraya-Watson estimator. We choose $\delta=1$.
In each round of iteration, the bandwidth $h_{n}$ is chosen as $h_{n}=\sigma_{\widehat{z}}\cdot n^{-1/5}$,
where $n$ is sample size, $\sigma_{\widehat{z}}$ is the standard
deviation of $z_{i,k}$, and $z_{i,k}=X_{0,i}+\mathbf{X}_{i}^{\mathrm{T}}\boldsymbol{\beta}_{k}$.
For SBGD estimator, we choose $q=9$ and use Legendre polynomials
with transformation discussed in \autoref{rem8}. For both estimators, the
stopping rule is either $\max_{1\leq j\leq p}|\widehat{\beta}_{j,k+1}-\widehat{\beta}_{j,k}|<10^{-5}$
or $k\geq20000$. The above also applies to our empirical analysis
in Section \ref{section6}. Trimming is ignored during all the simulations.
Due to the outliers of the simulation, we trim out the lower and upper
2\% simulation results and calculate the bias and RMSE. }{\footnotesize\par}
\end{table}
\autoref{tab1} 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$.
\begin{table}
\begin{centering}
\caption{\label{tab2}Sensitivity of KBGD and SBGD Estimators: Fixed Coefficients}
\begin{tabular}{c>{\centering}p{1.3cm}>{\centering}p{1.2cm}>{\centering}p{1.2cm}c>{\centering}p{1.2cm}>{\centering}p{1.2cm}>{\centering}p{1.2cm}}
\hline
& & \multicolumn{2}{c}{Sensitivity} & & \multicolumn{3}{c}{Running Time}\tabularnewline
\hline
& Method & $S_{L}$ & $S_{Z}$ & & True & Logit & Zeros\tabularnewline
\hline
\multirow{2}{*}{$n=2500$} & KBGD & 0.0242 & 0.0198 & & 113.21 & 79.120 & 158.91\tabularnewline
& SBGD & 0.0175 & 0.0259 & & 0.9504 & 0.9482 & 1.1587\tabularnewline
\hline
\multirow{2}{*}{$n=5000$} & KBGD & 0.0241 & 0.0175 & & 157.48 & 87.954 & 230.07\tabularnewline
& SBGD & 0.0189 & 0.0282 & & 1.4644 & 1.4722 & 1.9074\tabularnewline
\hline
\end{tabular}
\par\end{centering}
{\footnotesize{}NOTE: The running time is in seconds. Due to the outliers
of the simulation, we trim out the lower and upper 2\% simulation
results and calculate the corresponding results. The above also applies
to \autoref{tab3}.}{\footnotesize\par}
\end{table}
\autoref{tab2} 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.
\begin{table}
\centering{}\caption{\label{tab3}Sensitivity of KBGD and SBGD Estimators: Random Coefficients}
\begin{tabular}{c>{\centering}p{1.3cm}>{\centering}p{1.2cm}>{\centering}p{1.2cm}c>{\centering}p{1.2cm}>{\centering}p{1.2cm}>{\centering}p{1.2cm}}
\hline
& & \multicolumn{2}{c}{Sensitivity} & & \multicolumn{3}{c}{Running Time}\tabularnewline
\hline
& Method & $S_{L}$ & $S_{Z}$ & & True & Logit & Zeros\tabularnewline
\hline
\multirow{2}{*}{$n=2500$} & KBGD & 0.0270 & 0.0214 & & 122.00 & 74.433 & 166.94\tabularnewline
& SBGD & 0.0123 & 0.0246 & & 1.0132 & 0.8252 & 1.2044\tabularnewline
\hline
\multirow{2}{*}{$n=5000$} & KBGD & 0.0234 & 0.0232 & & 163.74 & 91.449 & 247.49\tabularnewline
& SBGD & 0.0077 & 0.0234 & & 1.5529 & 1.4377 & 1.9217\tabularnewline
\hline
\end{tabular}
\end{table}
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 \autoref{tab3}. 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.
\section{\label{section6}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 \citet{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) \citep{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 \autoref{empirical_table}. We
then conduct KBGD and SBGD estimation and report the estimated coefficients
of education in (II) and (III). As we can see from \autoref{empirical_table},
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.
\begin{table}
\caption{\label{empirical_table}Estimation Results}
\begin{centering}
\begin{tabular}{>{\raggedright}p{4cm}>{\centering}m{2.5cm}>{\centering}m{2.5cm}>{\centering}m{2.5cm}}
\hline
& (I) & (II) & (III)\tabularnewline
\hline
Estd. Coefficients & $2.5543^{***}$
$(0.1070)$ & $2.4832^{***}$
$(0.3638)$ & $2.4647^{***}$
$(0.3239)$\tabularnewline
Num. of Obs. & $26906$ & $26906$ & $26906$\tabularnewline
Estimation Methods & Logit & KBGD & SBGD\tabularnewline
Running Time & 1.4276 & 8573.1 & 40.9941\tabularnewline
Num. of Iteration & -- & 14996 & 12986\tabularnewline
\hline
\end{tabular}
\par\end{centering}
Note: For Logit regression, we report the coefficient of education
divided by that of total asset. For semiparametric estimation, we
normalize the coefficient of total asset to be 1. The standard deviations
are reported in the brackets below the coefficients. $^{***}$ indicates
significance at $1\%$ level. For both KBGD and SBGD estimators, we
choose $\delta_{k}=1$. For KBGD estimator, we choose $h_{n}=C\cdot n^{-1/5}$
with $C=C_{k}=\text{std}(z_{i,k})$, and use the fourth-order Epanechinikov
kernel. For SBGD estimator, we choose $q=9$ and use Legendre polynomials
with transformation discussed in \autoref{rem8}. The starting point of iteration
for both KBGD and SBGD estimators is chosen as the origin point with
all arguments being 0. The stopping rule is set as $\max_{1\leq j\leq p}|\widehat{\beta}_{j,k+1}-\widehat{\beta}_{j,k}|<\varrho$
with $\varrho=10^{-5}$. Finally, the running time is in second.
\end{table}
\pagebreak
\section{Conclusions}\label{conclude}
\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.
\newpage{}