Stochastic Learning of Semiparametric Monotone Index Models with Large Sample Size

Stochastic Learning of Semiparametric Monotone Index Models with Large Sample Size

Introduction

With the rapid development of technology in data collection and data storage, it's becoming more and more common nowadays for data analysts to deal with data set with extraordinary amount of observations. This offers the researchers unprecedented opportunities to more precisely understand the potential mechanism lurking behind the data, while on the same time brings about a series of new challenges. Among others, the key challenge is the heavy computational burdens that make the existing statistical methods numerically prohibitive. For example, when estimating a model using gradient-based iterative optimization procedure, the gradient of some objective function is repeatedly evaluated at a sequence of candidate parameters so that the optimal point can be numerically found. When the sample size is extremely large, even a single evaluation of the gradient would cost a huge amount of computation time, let along evaluating repeatedly at many points, making model estimation practically infeasible. Consequently, it's more urgent than ever before to study estimation methods that is applicable in the big-data era.

This paper studies semiparametric estimation of monotone index models in a large $n$ scenario. To fix idea, throughout this paper I will focus on the following binary choice model

where $\mathds{1}\left(\cdot\right)$ is indicator function, $\mathbf{X}_{e}=\left(X_{0},\mathbf{X}^{\mathrm{T}}\right)^{\mathrm{T}}=\left(X_{0},X_{1},\cdots,X_{p}\right)^{\mathrm{T}}\in\mathcal{X}_{e}$ is $(p+1)\times 1$ covariate vector, $\boldsymbol{\beta}_{e}^{\star}=\left(\beta_{0}^{\star},\boldsymbol{\beta}^{\star\mathrm{T}}\right)^{\mathrm{T}}=\left(\beta_{0}^{\star},\beta_{1}^{\star},\cdots,\beta_{p}^{\star}\right)^{\mathrm{T}}\in\mathcal{B}_{e}$ is the unknown true parameter vector, and $u$ is the unobserved individual shock with CDF $G\left(\cdot\right)$. Binary choice model is a leading example of the class of monotone index models, which has a wide range of applications in many areas such as economics, business, and biostatistics. I also point out that all of the conclusions obtained under such setup can be trivially extended to more general class of monotone index models.

When the CDF $G(\cdot)$ in ((ref)) is known, parametric estimation method such as maximum likelihood estimation can be applied. However, as I have discussed before, even under such setup estimation can be computationally costly when the data size is massive. To deal with the “large $n$” issue, subsample-based estimation strategy are widely applied. For example, when applying the gradient descent algorithm to iteratively search for the maxima of the log-likelihood function, instead of using the full sample, it's generally proposed to use a random subsample whose sample size is much smaller than $n$ to perform the update, which is known as the mini-batch gradient descent algorithm (MBGD, see bottou2018optimization,ruder2016overview). The batch size can be chosen as small as 1, in which case the algorithm is known as the stochastic gradient descent toulis2017asymptotic. For another example, forneron2022estimation studies stochastic optimization based on Newton-Raphson and quasi Newton iterations for a general class of parametric objective functions, and proposes subsample-based estimation and inference procedure for the unknown parameters.

In this paper, I focus on the semiparametric estimation of $\boldsymbol{\beta}_{e}^{\star}$. In other words, I seek to estimate $\boldsymbol{\beta}_{e}^{\star}$ without specifying the functional form of $G\left(\cdot\right)$. The main advantages of semiparametric specification are model flexibility as well as tractability. In the existing literature, semiparametric estimation for monotone index models and binary choice model in particular has been extensively studied. The methods can be roughly classified into two categories: M-estimation approach and direct construction approach. For the first category, the estimator is obtained by optimizing some objective functions. The standing estimators include maximum score estimator manski1975maximum,manski1985semiparametric,horowitz1992smoothed, maximum rank correlation estimator han1987non,sherman1993limiting,cavanagh1998rank,fan2020rank, semiparametric least squares estimator hardle1993optimal,ichimura1993semiparametric and semiparametric maximum likelihood estimator cosslett1983distribution,klein1993efficient. Apart from M-estimation, the second class of estimation methods features direct construction of the estimators, which includes average derivative estimator stoker1986consistent,powell1989semiparametric,horowitz1996direct,hristache2001direct, special regressor approach lewbel2000semiparametric and eigenvalue approach ahn2018simple.

The key feature that distinguishes my paper from the existing literature is that I try to estimate the model in a scenario where the sample size $n$ is extremely large. Large sample size $n$ imposes computational challenges to model estimation even in the parametric setup, and such issue turns out to be far more serious in the semiparametric setup. In his famous paper, ichimura1993semiparametric pointed out that for semiparametric least square estimator, “the computation time is roughly $n$ times more than with smooth parametric nonlinear regression estimation”. So if I estimate the semiparametric model based on a data set of millions of observations, the estimation time would be roughly millions of times longer than parametric estimation, say, Logit or Probit regression. This makes semiparametric estimation almost computationally infeasible when $n$ is extremely large. Indeed, for many semiparametric M-estimators such as ichimura1993semiparametric's semiparametric least squares estimator and klein1993efficient's semiparametric maximum likelihood estimator, the unknown CDF (or monotonic link function for more general monotone index models) $G(\cdot)$ in the objective function is replaced with its Nadaraya-Watson kernel estimator. So evaluating the objective function (or its gradient) generally involves calculating kernel estimators (or their gradients) at $n$ points. Since each kernel estimator (or its gradient) requires computational complexity of order $O(n)$, a single assessment of the objective function (or its gradient) requires computational time of order $O(n^2)$, which increases fast with the sample size $n$. This makes the conventional semiparametric estimation method not applicable even for data set with only tens of thousands of observation points. Apart from intensive computational burdens, there are many other crucial limitations that prohibit the use of existing semiparametric estimation methods\footnote{For M-estimation approach, the objective functions involved are usually heavily discontinuous and/or non-convex with respect to the parameter. In this case, even looking for a local optimum is generally NP-Hard murty1987some, let alone the global optimum. This makes the optimization procedure computationally infeasible. On the other side, the direct construction approach generally imposes more structure on the covariates. For example, the average derivative approach requires that the covariates are all continuous, so can not be directly applied to discrete covariates such as dummy variables. Moreover, the application of such method usually involves nonparametric estimation of the density functions or their partial derivative of some random variables conditional on the covariates. Such estimation becomes an intractable problem even when the number of covariates is modest. Although there have been some attempts to reduce the dimensionality of conditional density estimation (e.g., hall2004cross), the methods are still computationally-intensive, which may not be applicable in a data-rich environment, see ouyang2023 and references therein.}.

In this paper, I propose a novel semiparametric estimation procedure for ((ref)) that can be easily implemented with very fast speed even on a regular laptop when the sample size $n$ is extremely large. My method is motivated by the MBGD algorithm. For any random variable $Z$, parameter $\theta$, and loss function $L(Z,\theta)$, given a sequence of realizations $Z_1, \cdots, Z_n$ of $Z$, to search for the optimal point $\theta^{\star}$ that minimizes the population loss function $\mathbb{E}_Z\left(L(Z, \theta)\right)$, MBGD conducts the following iteration,

where $\theta_1$ is some initial guess, $\delta_k >0$ is the learning rate, and $\mathfrak{I}_{k}$ is the subsample used in the $k$-th round of iteration. In other words, the MBGD algorithm updates the parameter based on the gradient of the loss function at observation points that fall into the subsample $\mathfrak{I}_k$. Compared with the full-sample-based batch gradient descent (BGD) that uses gradient at all the data points to perform the update, MBGD update is less accurate\footnote{When using the full sample to conduct update, the gradient of the empirical loss function $L_n(\theta) = \frac{1}{n}\sum_{i=1}^nL(Z_i, \theta)$ is accurately evaluated at each candidate parameter $\theta$ because the gradient of the loss function at each data point $Z_i$ is evaluated. While when using subsample-based update, the gradient of the empirical loss function is only approximated by the gradients at a subsample of observations.} but significantly alleviates the computational burden when $\left|\mathfrak{I}_k\right|\ll n$. Typically, the MBGD algorithm applies only to the parametric setup where the loss function $L\left(\cdot, \cdot\right)$ is fully known. While when estimating the binary choice model ((ref)), the loss function generally has form $L\left(\left. \cdot,\cdot\right|G\right)$, so depends on the link function $G\left(\cdot\right)$\footnote{For example, the quadratic loss function is given by $L(\mathbf{X},y,\boldsymbol{\beta}|G) = (y - G(\mathbf{X}^{\mathrm{T}}\boldsymbol{\beta}))^2$ and the log-likelihood loss function is given by $L(\mathbf{X},y,\boldsymbol{\beta}|G) = -(y\log(G(\mathbf{X}^{\mathrm{T}}\boldsymbol{\beta})) + (1-y)\log(1-G(\mathbf{X}^{\mathrm{T}}\boldsymbol{\beta}))) $. }. In the semiparametric setup where $G$ is unknown, $L\left(\left. \cdot,\cdot\right|G\right)$ is then not fully specified, which makes the above MBGD update no longer feasible.

To make ((ref)) feasible, I consider a two-step updating procedure. In the $k$-th round of update, I first nonparametrically estimate the unknown function $G\left(\cdot\right)$, whose estimator is denoted as $\widehat{G}_k$. Then in the second step, I plug the first-step estimator $\widehat{G}_k$ into the loss function $L\left(\left.Z,\theta\right|G\right)$ and perform the update ((ref)) based on the estimated loss function $L\left(\left.Z,\theta\right|\widehat{G}\right)$ as if it were the true loss function. The key difficulty of such two-step update in the large $n$ scenario lies in the heavy computational burden caused by nonparametric estimation of $G\left(\cdot\right)$. Indeed, conventional nonparametric estimator such as Nadaraya-Watson kernel estimator requires computational complexity of order $O(n)$ to evaluate $\widehat G_k$ at a single point. So if I use a subsample of size $B$ to perform the update, I need to evaluate $\widehat G_k $ at a total of $B$ points, and the computational burden of each single update is of order $O(Bn)$. This is too large to be practical if I choose $B \gg 1/\sqrt{n}$\footnote{Indeed, this is required if I pursue $1/\sqrt{n}$-consistency and asymptotic normality of the estimator, see (ref).} and update hundreds of thousands of times. The main novelty of this paper is that instead of using conventional nonparametric estimator based on the full sample, I propose to use subsample to construct the Nadaraya-Watson kernel estimator, so that the above two-step update is fully subsample-based. The idea behind such subsample-based nonparametric estimation is intuitive: if I believe that using subsample for iteration leads to relatively accurate update, then the subsample-based nonparametric estimator should also be reasonably close to the one based on the full sample. When the subsample size is $B$, evaluating $B$ subsample-based Nadaraya-Watson kernel estimators requires computational complexity $O(B^2)$. This will be much smaller than $O(n^2)$ if I choose $B\ll n$. Indeed, I will show that as long as I properly choose the kernel function, $B$ can be chosen sufficiently close to $1/\sqrt{n}$, so the computational burden of update can be made close to $O(n)$, which is almost linear in $n$. This makes semiparametric estimation of monotone index models practically feasible when the sample size $n$ is large.

khan2022 (KLTY hereafter) also consider a similar two-step updating procedure. While the main difference between my method and theirs lies in that in KLTY, both the first-step nonparametric estimation and the second-step update are based on the full sample. Full-sample-based update increases the update accuracy, but as I discussed before, it leads to heavy computational burdens so is only applicable when the sample size is modest. Comparatively, the main novelty of my method lies in that I propose a fully subsample-based update which substantially improves the computation speed and can be easily applied when the sample size is extremely large. Roughly speaking, the relationship between my method and KLTY's method is similar to that between mini-batch gradient descent and batch gradient descent. Finally, similar to KLTY's method, my proposed method also overcomes the optimization issue of the M-estimator, see KLTY for more discussion.

I also develop the statistical properties of the above fully subsample-based two-step updating algorithm. Under some regularity conditions, I show that the proposed alogithm yields an asymptotically consistent estimator. However, its guaranteed convergence rate is slower than the parametric rate $1/\sqrt{n}$ if I choose $B\ll n$ to improve computational speed. Indeed, the guaranteed convergence rate will be even slower than rate $1/\sqrt{B}$, which is the convergence rate of conventional MBGD estimators. Such slower convergence rate is mainly caused by subsample-based nonparametric estimation in the first step. The subsample-based nonparametric estimator is no longer an unbiased estimator for the one based on the full sample, and such bias dampens the $1/\sqrt{B}$-convergence. I then decompose the bias. I find that the first-order bias have $1/\sqrt{n}$-trivial conditional mean (conditioned on the subsamples in the previous updates and the data set), while the second-order bias are uniformly $1/\sqrt{n}$-trivial as long as I update sufficiently many times. This motivates me to follow polyak1992acceleration and use average to eliminate the first-order bias and accelerate the convergence rate. In particular, after some burn-in rounds of updates, all the estimators produced during the following updates are averaged. I show that as long as the numbers of burn-in and follow-up updates are both large enough, the averaged estimator will converge at $1/\sqrt{n}$ rate and is asymptotically normally distributed. Such a result demonstrates that our subsample-based method not only improves the computational speed, it also maintains the estimation accuracy on the same time.

Since the subsample-based estimator is asymptotically normally distributed after averaging, inference on the true parameter can be conducted if some consistent estimator of the asymptotic covariance matrix is available. Unfortunately, when sample size $n$ is extremely large, estimating the covariance matrix based on the full sample also requires large amount of time because it involves evaluating a large number of nonparametric estimators. To faciliate the inference, I also propose a subsample-based estimator of the covariance matrix, which subtantially improves the computation speed. I show that the subsample-based estimator is a consistent estimator of the unknown covariance matrix, so the inference using such subsample-based estimator will be asymptotically valid.

The main contribution of this paper to the econometric literature is that I propose a computationally friendly algorithm that can be used to semiparametrically estimate the monotone index models when the sample size $n$ is extremely large. My new algorithm essentially generalizes the mini-batch estimation method to the semiparametric setup. It can be easily applied when there are hundreds of covariates and hundreds of thousands of or even millions of data points. Essentially, it bridges the gap between semiparametric estimation theories and empirical applications in the data-rich environment.

As an empirical illustration of my new method, I revisit the empirical results in helpman2008estimating. In their paper, helpman2008estimating use a parametric Probit model to study how the conditional probability of one country exporting to another is affected by a set of coutry-pair factors, and such estimation results are further embedded into a second-step estimation of the gravity equation. The full data set they use contains a total of 248060 observation points and 337 covariates including large number of country and year fixed effects, which features both large $n$ and $p$. Given that Probit estimation assumes that the random shock in the binary choice model has tail that decays at a fast speed, the estimation results could be biased if the true random shock has heavier tails, and in that case, the subsequent inference of the true parameter will also be invalid. Above discussion motivates semiparametric estimation, but given the size of the data set, the conventional semiparametric estimation are practically infeasible. In this paper I apply the proposed KMBGD estimation procedure to revisit the estimation results. The estimation and inference based on my method take around 8 hours and 0.8 hours respectively, which is practically feasible. Interestingly, compared with Probit distribution, I find that semiparametric estimation results are more in favor of a Logit distributed random shock in the sense that the KMBGD estimator is close to Logit estimator while differs significantly from Probit estimator. Such a result also highlights the use of semiparametric estimation as opposed to parametric estimation in applications.

The remainder of the paper is arranged as follows. In (ref), I formally introduce the two-step fully subsample-based updating algorithm. In (ref), I develop the asymptotic properties of the proposed algorithm. Then in (ref), I propose a subsample-based inference procedure. In (ref), I study the finite-sample performance of the proposed algorithm by conducting some Monte Carlo simulations. In (ref), I apply my new algorithm to revisit helpman2008estimating's Probit estimation results. Finally, (ref) concludes. All the proofs of the lemmas and theorems are arranged to the Appendix.

Notations

For any real sequences $\left\{ a_{n}\right\} _{n=1}^{\infty}$ and $\left\{ b_{n}\right\} _{n=1}^{\infty}$, I 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}$, I write $a_{n}=O_{p}\left(b_{n}\right)$ if for any $0<\tau<1$ there exist $N$ and $C>0$ such that $P\left\{ \left|a_{n}/b_{n}\right|>C\right\} <\tau$ holds for all $n\geq N$, I write $a_{n}=o_{p}\left(b_{n}\right)$ if for any $C>0$, $P\left( \left|a_{n}/b_{n}\right|>C\right)\rightarrow0$. For any Borel set $A\subseteq\mathbb{R}^{k}$, denote its Lebesgue measure as $m\left(A\right)$. Denote $I_p$ as the $p$-dimensional identity matrix. 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$, I write $A\succ B$ if $A-B\succ0$ and $A\succeq B$ if $A-B\succeq0$. For any matrix $A$, I 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$, I 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}}$, I 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}$, I denote $\left\Vert A\right\Vert =\sqrt{\sum_{i=1}^{n}\sum_{j=1}^{m}a_{ij}^{2}}$.

The KMBGD Algorithm

This section formally introduces the subsample-based learning algorithm for binary choice models. To make my illustration more intuitive, I will start with a special case where the CDF function $G(\cdot)$ is known. Given any loss function $L\left(\mathbf{X}_{e},y,\boldsymbol{\beta}_{e}|G\right)$ that depends on $G(\cdot)$ and is differentiable with respect to $\boldsymbol{\beta}_{e}\in\mathcal{B}_{e}$, the conventional MBGD estimator of $\boldsymbol{\beta}_{e}^{\star}$ is constructed based on the following iteration bottou2018optimization,ruder2016overview,

where $\boldsymbol{\beta}_1$ is given, $B$ is a positive integer and is the sbusample size. For each $k$, $\delta_{k}>0$ is the learning rate, and

is an index set that is randomly drawn from $\{1,2,\cdots,n\}$ with replacement and is independent over $k$. In other words, under MBGD algorithm, in each iteration I randomly draw a subset of size $B$, and then update the estimator based on such subsample.

Given a choice of the subsample size $B$, to apply the MBGD algorithm ((ref)) to estimate $\boldsymbol{\beta}^{\star}$, it remains to choose the loss function. Following agarwal2014least and khan2022, I consider the loss function

for some sufficiently large positive constant $A$. khan2022 show that loss function ((ref)) has many properties such as global minimization at true parameter $\boldsymbol{\beta}^{\star}$ and positive definite Hessian matrix with respect to $\boldsymbol{\beta}_e$. Based on the MBGD updating rule ((ref)) and loss function ((ref)), the MBGD estimator of $\boldsymbol{\beta}_{e}^{\star}$ is constructed based on the following iteration procedure:

Now I turn to the case of semiparametric estimation, which is the main focus of this paper. To ensure identification, I set $\beta_{0}^{\star}$ to be 1, so the estimation target now is $\boldsymbol{\beta}^{\star}$. To simplify notation, denote the space of $\mathbf{X}$ as $\mathcal{X}$, and the corresponding parameter space of $\boldsymbol{\beta}$ as $\mathcal{B}$.

Note that the MBGD algorithm ((ref)) relies on the nonparametric component $G\left(\cdot\right)$ as a key input, which is unavailable in the current semiparametric setup. So the conventional MBGD algorithm is infeasible. To make the update feasible, a natural idea is to replace the unknown component with its nonparametric estimator. Intuitively, suppose that in the $k$-th round of iteration, the starting point $\boldsymbol{\beta}_{k}$ is close to the unknown true parameter $\boldsymbol{\beta}^{\star}$, then there holds \[G\left(z\right) =\mathbb{E}\left(\left.y\right|X_{0}+\boldsymbol{\mathbf{X}}^{\mathrm{T}}\boldsymbol{\beta}^{\star}=z\right)\approx\mathbb{E}\left(\left.y\right|X_{0}+\boldsymbol{\mathbf{X}}^{\mathrm{T}}\boldsymbol{\beta}_{k}=z\right),\] for any $z$. This immediately motivates the following Nadaraya-Watson kernel estimator for $G\left(\cdot\right)$,

where $K_{h}\left(\cdot\right)=h^{-1}K\left(\cdot/h\right)$, $K\left(\cdot\right)$ is kernel function, and $h_{n}$ is bandwidth parameter depending on $n$. Given the estimated CDF $\widehat{G}\left(\left.\cdot\right|\boldsymbol{\beta}_{k}\right)$, we can directly plug it back to ((ref)) and perform the update as if it were the true CDF $G\left(\cdot\right)$. Note that a potential issue for ((ref)) is that it's based on the full data set, so evaluating its value has computational complexity of order $O(n)$ for each input $z$. If I use $B$ data points to perform the update, then a total of $B$ kernel estimators need to be evaluated in each update, which leads to computational burden of order $O(nB)$. The computational speed can be improved if I choose $B\ll n$, but note that to obtain an estimator with $1/\sqrt{n}$-consistency, it is generally required that $B\sim \sqrt{n}$, see forneron2022estimation. Indeed, in the current semiparametric setup, the order of $B$ has to be chosen even slightly larger, see the following (ref). In this case, the computational burden will be of order at least $O(n\sqrt{n})$, which is far from being linear in $n$.

The key philosophy of my new algorithm is that, if I trust that using $B$ data points provides relatively accurate updates, then the kernel estimation based on such $B$ points should also be reasonably close to that based on the full sample for all input $z$. Such an idea motivates me to use only the randomly-drawn subset to construct the kernel estimator. In particular, consider the following Nadaraya-Watson kernel estimator of $G(z)$ constructed based on the data points in subsample $\mathfrak{I}_{B,k}$,

where $K_h$, $K$ and $h_n$ are all similarly defined as before, and $\underline{c}_{f}>0$ is some sufficiently small constant. Basically, the subsample-based estimator ((ref)) is constructed as if I only observe the random subsample $\left\{ \left(\mathbf{X}_{e,i},y_{i}\right)\right\} _{i\in\mathfrak{I}_{B,k}}$. The computational complexity for evaluating $\widehat{G}\left(\left.z\right|\boldsymbol{\beta},\mathfrak{I}_{B,k},\underline{c}_{f}\right)$ is obviously of order $O\left(B\right)$.

Given the subsample-based kernel estimator, I can formally illustrate my subsample-based learning algorithm. At the beginning of the $k$-th update, the initial point $\boldsymbol{\beta}_{k}$ is given. Then using the subsample-based kernel estimator of $G\left(z\right)$ given in ((ref)), I consider the following updating algorithm,

where $\mathbf{X}_{i}^{\phi} = \mathbf{X}_{i}\cdot \mathds{1}(\boldsymbol{\mathbf{X}}_{e,i}\in\mathcal{X}_{e}^{\phi})$, and $ \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\} $ for some $0<\phi<1$\footnote{Such truncation is basically used to improve the uniform convergence speed of kernel estimation. Similar method is applied in many research such as ichimura1993semiparametric and klein1993efficient. }. Since the above algorithm generalizes the conventional mini-batch gradient descent procedure to the semiparametric setup, I label the new algorithm the kernel-based mini-batch gradient descent algorithm (KMBGD). The algorithm is summarized in (ref).

\IncMargin{1em}

\DecMargin{1em}

Statistical Properties of KMBGD Estimator

In this section, I formally study the statistical properties of the proposed KMBGD estimator. Under some regularity conditions, I first show that as long as I update sufficiently many times, the KMBGD estimator is consistent. However, the convergence rate is slower than $1/\sqrt{n}$ if I choose $B\ll n$. Indeed, such rate is even slower than $1/\sqrt{B}$, which is the convergence rate of general mini-batch estimators forneron2022estimation. Then I will show that although KMBGD estimator itself converges at a slow rate, I can conduct averages across all the estimators produced during updates to accelerate the convergence rate. In particular, I show that if we properly choose subsample size, bandwidth prameter, order of kernel function, and number of iterations, the average estimator obtains $1/\sqrt{n}$-consistency.

Before I illustrate the main results, I first introduce some notations. Let $f_e\left(\mathbf{X}_e\right)$ and $f\left(\mathbf{X}\right)$ denote the joint density of $\mathbf{X}_e$ and $\mathbf{X}$\footnote{By assuming $\mathbf{X}_e$ has joint density function, we require that $\mathbf{X}_e$ is continuous, which facilitates our following discussion. However, I point out that my analysis can be trivially extended to the case where there are some discrete covariates, see KLTY. }. Define $z\left(\mathbf{X}_e, \boldsymbol{\beta}\right) = X_0 + \mathbf{X}^{\mathrm{T}}\boldsymbol{\beta}$. Let $f_{\mathbf{X}|z}\left(\left.\mathbf{X}\right|z,\boldsymbol{\beta}\right)$ be the conditional density of $\boldsymbol{\mathbf{X}}$ given $z\left(\mathbf{X}_e, \boldsymbol{\beta}\right) =z$ and $\boldsymbol{\beta}$. Define

The following technical assumptions are imposed.

Based on the above assumptions, now I formally study the statistical properties of the iterative estimator $\boldsymbol{\beta}_{k}$ based on iteration ((ref)) and ((ref)). I first introduce some further notations. Let $P$ denote the probability measure of the data set $\mathcal{D}_n$. Let $\mathbb{P}^{*}$ be the probability measure corresponding to random variables $\{\mathfrak{I}_{B,k}\}_{k=1}^{\infty}$ and $\mathbb{P}_{k}^{*}$ be probability measure corresponding to $\{\mathfrak{I}_{B,k^{\prime}}\}_{k^{\prime}\geq k}^{\infty}$ conditional on the observation of $\{\mathfrak{I}_{B,k^{\prime}}\}_{k^{\prime}=1}^{k-1}$ for $k\geq2$ and $\mathbb{P}_{1}^{*}=\mathbb{P}^{*}$. Let $\mathbb{E}^{*}$ and $\mathbb{E}_{k}^{*}$ be the expectation with respect to $\mathbb{P}^*$ and $\mathbb{P}_{k}^{*}$. Finally, let $\mathbb{P}$ be the probability measure of $\{\mathscr{D}_{n}, \mathfrak{I}_{B,1}, \mathfrak{I}_{B,2}, \cdots\}$, where $\mathscr{D}_{n}$ is the data set.

Recall that the Nadaraya-Watson kernel estimator for $\mathbb{E}\left(\left.y\right|X_{0}+\mathbf{X}^{\mathrm{T}}\boldsymbol{\beta}=z\right)$ based on the full data is given by $\widehat{G}\left(\left.z\right|\boldsymbol{\beta}\right)$ in ((ref)). For any $\boldsymbol{\beta}\in\mathcal{B}$, define $\Delta \boldsymbol{\beta} = \boldsymbol{\beta}- \boldsymbol{\beta}^{\star}$. I obviously have the following decomposition for the MBGD update ((ref)),

It's not difficult to see that if $\pi_{1,n,k} = \pi_{2,n,k} =0$, then ((ref)) degenerates to the full-sample-based KBGD algorithm. Indeed, $\pi_{1,n,k}$ describes the randomness caused by updating using only a subset of the data, whereas $\pi_{2,n,k}$ describes the randomness caused by performing nonparametric kernel estimation using only a subset of the data points. Essentially, $\pi_{1,n,k}$ is shared by all the mini-batch estimators, while $\pi_{2,n,k}$ is specific to the semiparametric setup I consider in this paper. I have the following lemma describing the properties of $\pi_{1,n,k}$ and $\pi_{2,n,k}$.

(ref) immediately yields the following result.

According to (ref), if I choose $B\ll n$ to improve computational speed, the upper bounded on the estimation error $\mathbb{E}^{*}\left(\left\Vert \Delta\boldsymbol{\beta}_{k}\right\Vert \right) $ will be of rate slower than $n^{-1/2}$ even when the order of the kernel function is large. The slower convergence rate is a common feature of all the mini-batch estimators. Indeed, the mini-batch estimators converge at the rate $1/\sqrt{B}$ at best, see, for example, Lemma 2 in forneron2022estimation. However, different from the conventional mini-batch estimator, my KMBGD estimators are guaranteed to converge no faster than $\sqrt{\log(n)/Bh_n^2}$. If I choose $B=1/\sqrt{n}$ and $h_n = n^{-1/6}$, then the convergence rate would be $\sqrt{\log(n)}n^{-1/12}$, which is much slower than $1/\sqrt{B}=n^{-1/4}$.

The slower convergence rate of the KMBGD estimator is mainly due to the fact that I use subsamples to construct the kernel estimator. In this case, the subsample-based gradient is no longer an unbiased estimator (conditional on the previous subsamples) of the full-sample-based gradient, that is, $\mathbb{E}^*(\pi_{2,n,k})\neq 0$. The bias makes the convergence rate of KMBGD estimator slower than $1/\sqrt{B}$. However, surprisingly, in the following I will show that if I appropriately choose the kernel function and bandwidth parameter, even with $B\ll n$, I can still obtain $1/\sqrt{n}$ by following polyak1992acceleration and conducting average across KMBGD estimators produced during iterations.

To formally show the above results, I first further decompose the KMBGD dynamics. To ease my following exposition, for any $z$ and $\boldsymbol{\beta}$ denote $A_{n,y}\left(z,\boldsymbol{\beta}\right)=\frac{1}{n}\sum_{i=1}^{n}K_{h_{n}}\left(z-X_{0,i}-\mathbf{X}_{i}^{\mathrm{T}}\boldsymbol{\beta}\right)y_{i}$, $A_{n,1}\left(z,\boldsymbol{\beta}\right)=\frac{1}{n}\sum_{i=1}^{n}K_{h_{n}}\left(z-X_{0,i}-\mathbf{X}_{i}^{\mathrm{T}}\boldsymbol{\beta}\right) $, $A_{n,y}\left(\left.z,\boldsymbol{\beta}\right|\mathfrak{I}_{B,k}\right)=\frac{1}{B}\sum_{i\in\mathfrak{I}_{B,k}}K_{h_{n}}\left(z-X_{0,i}-\mathbf{X}_{i}^{\mathrm{T}}\boldsymbol{\beta}\right)y_{i}$, and $A_{n,1}\left(\left.z,\boldsymbol{\beta}\right|\mathfrak{I}_{B,k}\right)=\frac{1}{B}\sum_{i\in\mathfrak{I}_{B,k}}K_{h_{n}}\left(z-X_{0,i}-\mathbf{X}_{i}^{\mathrm{T}}\boldsymbol{\beta}\right) $. I have the following lemma.

I now provide some discussion for (ref). Basically, if there are no noise terms $\varrho_{1,n,k}$, $\varrho_{2,n,k}$, and $\varrho_{3,n,k}$, then the dynamics of $\Delta \boldsymbol{\beta}_k$ simply degenerate to the full-sample-based KBGD algorithm in KLTY as implied in (ref) in Appendix. However, since I use subsamples to perform the update, additional noises due to subsampling are introduced into the update and these noises are captured by the above three terms. Basically, $\varrho_{1,n,k}$ describes the impacts of using subsamples instead of full sample to perform the update. Such error is shared by all the mini-batch-based methods. While the remaining two terms $\varrho_{2,n,k}$ and $\varrho_{3,n,k}$ describe the impacts of using subsamples instead of full sample to construct the Nadaraya-Watson kernel estimator, so are specific to my algorithm only. Simple calculation leads to $ \mathbb{E}^* \left(\varrho_{1,n,k}\right) = 0, $ $\mathbb{E}^{*}\left(\varrho_{2,n,k}\right) = O_p\left(1/Bh_n\right), $ and $\mathbb{E}^{*}\left(\varrho_{3,n,k}\right) = O_p\left(1/Bh_n\right) $ uniformly with respect to $k$. The above implies that for $k$ sufficiently large, the first-order difference between KBGD and KMBGD estimators almost constitute a martingale difference sequence. By “almost” I mean that the conditional expectation is of order $O_p(1/Bh_n)$, which can be made $n^{-1/2}$-trivial if I choose $B\gg n^{1/2}h_n^{-1}$.

(ref) implies that although the KMBGD estimator itself does not obtain $1/\sqrt{n}$-consistency due to noises caused by subsample-based kernel estimation and update, I can follow polyak1992acceleration to conduct average across the estimators produced during iterations to eliminate these noises. Similar to the conventional mini-batch gradient estimator, the resulting estimator will be $1/\sqrt{n}$-consistent as long as we choose $B$ that diverges at some rate. In particular, let $k^*$ be the number of burn-in iterations and $T$ be the number of follow-up iterations. The averaged KMBGD estimator (AKMBGD) is defined as follws,

I summarize the algorithm in (ref).

\IncMargin{1em}

\DecMargin{1em}

Now I provide the theoretical properties of the AKMBGD estimator.

(ref) is the key result of this paper. It demonstrates that even though I only use a random subsample whose size is substaintially smaller than the full sample size to conduct kernel estimation and perform update in each round of iteration, the average of estimators produced during iterations will be equivalent to the full-sample estimator up to some small order terms. The small order terms will be uniformly $1/\sqrt{n}$-trivial as long as I choose $B\gg \max\{\log^2(n)h_n^{-6},\sqrt{n}\log(n)h_n^{-2}\}$ and $T \gg nB^{-1}h_n^{-2}$. This implies that as long as I choose kernel function properly, my KMBGD estimator will be as efficient as the one based on the full sample, dispite the fact that I only use a much smaller subsample to perform the update in each round.

(ref) also suggests that the computational speed of each update can be improved by appropriately choosing the kernel function. In particular, since $h_n$ must satisfy $h_n \ll n^{-1/2D}$ according to the conditions required in the theorem, then $B\gg \max\{n^{3/D}\log^2(n), n^{1/2+1/D}\log (n)\}$ must hold, so the computational complexity will be of order at least $O(\max\{n^{6/D}\log^4(n), n^{1+2/D}\log^2 (n)\})$. Obviously, to improve the computational speed, I can choose a high-order kernel function. For example, if I choose a 8-th order kernel, the computational complexity is of order $O(n^{5/4}\log^2(n))$; if I choose a 12-th order kernel, the computational complexity is of order $O(n^{7/6}\log^2(n))$. If I can choose sufficiently large $D$, then the computational complexity is lower bounded by $n\log^2(n)$, which is almost the linear rate $O(n)$.

I finally discuss the total computational time of KBGD and KMBGD estimation. Suppose $k^*$ updates are necessary to eliminate the impacts of the initial guess, then the full-sample-based KBGD algorithm requires $O(k^*n^2)$ computational time in total, while the KMBGD algorithms requires $O(k^*B^2 + B^2T)$. Since (ref) requires that $T\gg nB^{-1}h_n^-2$, then the total computational time of KMBGD will be at least $O(k^*B^2 + nBh_n^{-2})$. If I choose $B\gg \sqrt{n}h_n^{-2}\log n$ and $h_n \ll n^{-1/2D}$, then $k^*B^2 + nBh_n^{-2}\gg k^*n^{1+2/D}\log^2(n) + n^{3/2+2/D}$. So the upper bound on the ratio between the total computational time of KBGD and KMBGD is of order \[ n^{1-2/D}\log^{-2}(n) + k^*n^{1/2 - 2/D}. \] Obviously, when $D\geq6$, the above ratio diverage at rate $n^{2/3} + k^*n^{1/6}$. More crucially, the above rate will be large when $k^*$, the number of burn-in updates, is large, which will often be the case when the number of covariates is large and $\underline{\varLambda}/\overline{\varLambda}$ is small,

Inference with Large $n$

In this section, I discuss the inference-related issues when the sample size $n$ is large. According to (ref), the AKMBGD estimator is asymptotically normally distributed, so inference on the true parameter $\boldsymbol{\beta}^{\star}$ can be conducted if I can consistently estimate the asymptotic covariance matrix $\Sigma_{\boldsymbol{\beta}}^{\phi}$. In their paper, KLTY provide a consistent estimator for the covariance matrix based on the full sample. However, to construct such estimator, I need to construct nonparametric estimators for conditional expectation $\mathbb{E}\left(\left. \mathbf{X}_{i}^{\phi}\right|z_i^{\star}\right)$ for each $i$, which may cost large amount of time when both $n$ and $p$ are large.

For parametric optimization, forneron2022estimation proposes a stochastic Newton-Raphson udpate and use the produced estimators for inference to alleviate the computational burden of statistical inference. But his method can not be applied in the current scenario even if I can approximate the “Hessian” matrix\footnote{Note that in our case, the “Hessian” refers to the matrix $\varLambda_{\phi} (\boldsymbol{\beta}^{\star})$, which is actually not symmetric.} accurately. This is because, apart from $\varrho_{1,n,k}$ that captures the distribution of $\boldsymbol{\xi}^{\phi}_n$, additional subsampling errors $\varrho_{2,n,k}$ and $\varrho_{3,n,k}$ are introduced because I use subsamples to construct the nonparametric estimator. Such additional errors are at least of the same order as $\varrho_{1,n,k}$, so they dampen the bootstrap-based inference.

To solve the above inference issue in the large $n$ scenario, this section provides a subsample-based estimator for the covariance matrix. Let $\{\mathfrak{I}_{B,r}\}_{r=1}^R$ be a sequence of random index sets defined in ((ref)). For each $1\leq r\leq R$, define \[ \widehat{\Sigma}_{\boldsymbol{\xi}}^{\phi, r}=\frac{1}{B}\sum_{i\in \mathfrak{I}_{B,r}}\left(\widehat{G}_{i}^r\left(1-\widehat{G}_{i}^r\right)\left(\mathbf{X}_{i}^{\phi}-\widehat{\mathbb{E}}^r\left(\left.\mathbf{X}_{i}^{\phi}\right|\widehat{z}_{i}\right)\right)\left(\mathbf{X}_{i}^{\phi}-\widehat{\mathbb{E}}^r\left(\left.\mathbf{X}_{i}^{\phi}\right|\widehat{z}_{i}\right)\right)^{\mathrm{T}}\right), \] and \[ \widehat{\varLambda}^r_{\phi}\left(\overline{\boldsymbol{\beta}}\right)=\frac{1}{B}\sum_{i\in\mathfrak{I}_{B,r}}\mathbf{X}_{i}^{\phi}\frac{\partial\widehat{G}\left(\left.z\left(\mathbf{X}_{e,i},\overline{\boldsymbol{\beta}}\right)\right|\overline{\boldsymbol{\beta}}, \mathfrak{I}_{B,r}, \overline{c}_f\right)}{\partial\boldsymbol{\beta}^{\mathrm{T}}}, \] where \[ \widehat{G}_{i}^r=\frac{\frac{1}{B}\sum_{j\in\mathfrak{I}_{B,r}}K_{h_{n}}\left(\widehat{z}_{i}-\widehat{z}_{j}\right)y_{j}}{\left\{\frac{1}{B}\sum_{j\in\mathfrak{I}_{B,r}}K_{h_{n}}\left(\widehat{z}_{i}-\widehat{z}_{j}\right)\right\}\lor \overline{c}_f},\ \widehat{\mathbb{E}}^r\left(\left.\mathbf{X}_{i}^{\phi}\right|\widehat{z}_{i}\right)=\frac{\frac{1}{B}\sum_{j\in\mathfrak{I}_{B,r}}K_{h_{n}}\left(\widehat{z}_{i}-\widehat{z}_{j}\right)\mathbf{X}_{j}^{\phi}}{\left\{\frac{1}{B}\sum_{j=\in\mathfrak{I}_{B,r}}K_{h_{n}}\left(\widehat{z}_{i}-\widehat{z}_{j}\right)\right\}\lor \overline{c}_f}, \] and $\widehat{z}_{i}=X_{0,i}+\mathbf{X}_{i}^{\mathrm{T}}\overline{\boldsymbol{\beta}}$. Also define

Then we have the following result.

Monte Carlo Experiments

This section conducts some Monte Carlo experiments to evaluate the finite-sample performance as well as the computational efficiency of the proposed KMBGD and AKMBGD estimators. Throughout this section, I consider the following data generating process

where $n$ is the sample size. For all $1\leq i\leq n$, $X_{0,i}\sim \mathcal{N}(0,1)$, $X_{1,i}\sim \text{Bernoulli}(1/2)$, $X_{2,i}\sim \text{Poisson}(2)$, and $X_{j,i}\sim (\chi^{2}(1)-1)/\sqrt{2}$ for $3\leq j \leq 9$. So I have a mixture of both continuous and discrete covariates. Moreover, $X_{j,i}$ is independent over $j$ for each $i$. $u_i$ is the random error with cumulative distribution function $G(u)$, which is independent of the covariates. $(X_{0,i}, \cdots, X_{9,i}, u_i)$ is iid over $i$. I set the true parameter vector as $\boldsymbol{\beta}^{\star}=\left(1,1,0.5,2,5,-0.5,-1,-2,-5\right)^{\mathrm{T}}.$ I consider four setups of error distrubtion: $\text{Cauchy}$, $t(4)$, $\chi^2(3)$, and $\mathcal{N}(0,1)$. Finally, in the following simulations, whenever I conduct the kernel estimation, I use sixth-order Epanechnikov kernel to construct the Nadaraya-Watson estimator, where the kernel function is given by $K(u) = \frac{525}{256}\left(1 - u^2\right)\left(1-6u^2 - \frac{33}{5}u^4\right)\mathds{1}\left(|u|\leq 1\right)$.

Finite-Sample Performance

In this subsection, I conduct some Monte Carlo experiments to study the finite sample performance of our AKMBGD estimator. I consider three setups of sample sizes: $n= 25000$, $n = 50000$, and $n = 100000$. I report the bias, root mean squared error (RMSE), and coverage rate of AKMBGD estimators for $\beta_{1}^{\star}$ to $\beta_{9}^{\star}$. Suppose that the simulation is repeated $R$ times, in the $r$-th round the estimator of $\beta_{j}^{\star}$ is denoted as $\widehat{\beta}_{j}^{r}$. Then the bias and RMSE of $\beta_{j}^{\star}$ is defined by \[ \text{Bias}=\left|\frac{1}{R}\sum_{r=1}^{R}\widehat{\beta}_{j}^{r}-\beta_{j}^{\star}\right|,\ \ \text{RMSE}=\sqrt{\frac{1}{R}\sum_{r=1}^{R}\left(\widehat{\beta}_{j}^{r}-\beta_{j}^{\star}\right)^{2}}. \] I consider nominal coverage rate $0.95$, so the actual coverage rate is given by \[ \mathrm{CR}=\frac{1}{R}\sum_{r=1}^{R}\mathds{1}\left(\widehat{\beta}_{j}^{r}-1.96\widehat{\sigma}_{j}^{r}\leq\beta_{j}^{\star}\leq\widehat{\beta}_{j}^{r}+1.96\widehat{\sigma}_{j}^{r}\right), \] where $\widehat{\sigma}_{j}^{r}$ is the subsample-based estimator of the variance of $\widehat{\beta}_{j}^{r}$.

The learning rate is chosen as $\gamma_k = 1$ for all $k$. The bandwidth used in the $k$-th round of update is $h_n = c_k\cdot h_n^{-1/10}$, where $c_k = \text{std}\left(z_{i,k}\right)$ and $z_{i,k} = X_{0,i}+\mathbf{X}_i^{\mathrm{T}}\boldsymbol{\beta}_k$. The initial guess is chosen as the Logit estimator. When constructing the AKMBGD estimator, I first run 2000 burn-in updates. Then the stopping rule is chosen as that in (ref) with $T = 10000$, $gap = 1000$, and $\varrho = 0.001$. The subsample size $B$ is chosen as 3000 for both estimation and inference. Finally, when conducting inference, i randomly draw 200 subsamples to construct the variance estimator.

The simulation results are reported in (ref). It can be seen that the AKMBGD estimators have small bias, whose RMSE decreases with sample size almost at rate $\sqrt{n}$. Moreover, the confidence interval constructed based on the subsample-based variance has actual coverage rate that is quite close to the nominal rate $0.95$. This demonstrates that the AKMBGD estimators and subsample-based variance estimator have great finite-sample performance.

Computational Efficiency

This subsection formally compares the computational efficiency of several gradient-based estimators for semiparametric montone index models. In particular, I compare KMBGD estimator with the KBGD and SBGD estimators proposed by khan2022.

I first compare the updating speed of each algorithm under different setups of sample sizes. In particular, for each algorithm, I keep updating 100 times and report the average running time of each single update. For kernel-based updates (KBGD and KMBGD), I consider two computation strategies: unparalleled and parallel computation. When using parallel computation, kernel estimators are simultaneously calculated over 6 cores. I consider six sample sizes: $n=2500, 5000, 10000, 20000, 500000,$ and $1000000$. For SBGD estimation, the sieve functions follow those used in khan2022. The order of sieves is chosen as $q=9$ when $n=2500$ and $5000$, $q=11$ when $n=10000$ and $20000$, and $q=31$ when $n=500000$ and $1000000$. The subsample size $B$ is chosen as $B=1000$ when $n\leq 20000$, $B=3000$ for $n=500000$, and $B = 5000$ for $n=1000000$. The simulation results are reported in (ref).

It can be seen that without parallel computation, the updating time of full-sample-based KBGD algorithm increases roughly at rate $n^2$, which is in linear with the previous discussion. In particular, when sample size is 2500, each single update requires 0.0475 seconds, which amounts to 21 updates within one second. However, such updating time increases to 0.2 seconds when sample size is 5000, which amounts to only 5 updates each second. When the sample size is 20000, without parallel computation, each single update of KBGD requires more than 3 seconds, indicating that 1000 updates may cost around 1 hour of computational time. For extremely large sample sizes $n=500000$ or $1000000$, KBGD is practically infeasible, so the computational time is not reported. It can also be seen that parallel computation may significantly decrease the updating time when $n$ is large ($n=10000, 20000$), but the updating time is still too long to be practically feasible.

I then look at the updating speed of SBGD and KMBGD. Apparently, when sample size is small or modest, SBGD exhibits excellent performance: when sample size is 2500, 5000, and 10000, each single update of SBGD requires only 0.0003, 0.0004, and 0.0006 seconds, which amounts to 3300, 2500, and 1600 updates within one second. Even when sample size is 20000, each update of SBGD requires only 0.0027 seconds, so 370 updates can be conducted within one second. This suggests that SBGD significantly outperforms KMBGD when the sample size $n$ is small or modest. However, when the sample size $n$ is extremely large, KMBGD starts dominating SBGD. In particular, when $n=500000$ and 1000000, the updating speed of KMBGD (with parallel computation) is roughly 4 and 5 times faster than that of SBGD.

Of course, the reduction of computational time of each single update of KMBGD compared with that of SBGD may come at the cost of longer total running time or large estimation error. To study whether it is the case, I then compare the total running time of SBGD and KMBGD. I also consider four setups of random error distributions as I did in (ref). I consider two extreme sample sizes: $n=500000$ and $n=1000000$. The subsample size $B = 3000$ when $n=500000$ and $B = 5000$ when $n=1000000$. The stopping rule for SBGD is $\max_{1\leq j\leq 9}|\beta_{j,k+1} - \beta_{j,k}|<10^{-6}$ and that for KMBGD is the same as before. For both updates, the initial guess is located at Logit estimator, and the maximum number of updates is 20000. For inference, I choose subsample size $B = 3000$ when $n=500000$ and $B = 6000$ when $n=1000000$. The number of subsamples is chosen as 200. Finally, I note here that for both estimation and inference, unparalleled computation is considered.

I report the RMSE and running time of both estimation and inference in (ref). As can be seen from the table, for all combinations of error distributions and sample sizes, the RMSE of SBGD and KMBGD are almost identical, indicating that updates based on subsamples do not result in loss of estimation accuracy. When looking at the running time, it's impressive to see that, the estimation time of KMBGD is substantially shorter compared with that of SBGD. When $n=500000$, KMBGD decreases the running time by roughly half, while when $n$ increases to 1000000, the reduction of estimation time is more significant: running time of KMBGD is only around one forth of that of SBGD. It is also interesting to see that, when the sample size increases and I use a larger subsample size, the running time of KMBGD even slightly decreases. This implies that although using a larger subsample size may make updating speed slightly slower, it makes convergence faster because the amount of noises in the update is decreased.

I finally look at the computational burden of inference based on different methods. As can be seen from (ref), the operational time of variance calculation of SBGD is over 3.2 hours without parallel computation when $n=500000$, and it rises to around 14 hours when $n=1000000$. This implies that even SBGD may have adequate computational efficiency in terms of estimation, it may still cost a large amount of time to conduct inference. When turning to the subsample-based infernece under KMBGD, it can be clearly seen that variance estimation only requires around 0.1 hours (10 min) when $n=500000$ and 0.4 hours (40 min) when $n=1000000$, which significantly improves the speed of inference. I also report in (ref) the true standard deviation and subsample-based estimator of the standard deviation of each estimator, which are close to each other. This implies that subsample-based inference improves the speed while does not suffer from much accuracy loss.

Empirical Illustration

In this section, I will illustrate the empirical applicability of the new subsample-based learning method by revisiting some empirical results in helpman2008estimating. In their paper, helpman2008estimating consider estimating the following model,

where $T_{ij}$ is an indicator of whether country $j$ exports to country $i$, $\xi_j^{\star}$ is the exporter fixed effect of the $j$-th country, $\zeta_i^{\star}$ is the importer fixed effect of the $i$-th country, $d_{ij}$ is the natural logarithm of the geographic distance between countries $i$ and $j$, and $\phi_{ij}$ is a vector of covariates that describe the variable country-pair fixed trade cost. The full sample contains a total of 248060 observations and 338 covariates, which features both large $n$ and $p$. The covariates contain 12 key variables including Distance, Land Border, Island, Landlock, Legal, Language, Colonial Ties, Currency Union, FTA, Religion, WTO (none) and WTO (both), and 158 exporter fixed effects, 158 importer fixed effects, and 10 year fixed effects.

When estimating ((ref)) based on the full sample, helpman2008estimating consider a parametric Probit setup, where $G$ is specified to be the CDF of standard normal distribution. In this section, I reestimate model ((ref)) without assuming the functional form of $G$ by applying the KMBGD algorithm. Such reestimation is well motivated because assuming normal distributed random shocks actually makes restrictive assumptions over the decreasing speed of the tails of the random shocks, which might be violated in some empirical applications. Misspecification of distribution of random shocks may dampen the estimation results as well as the subsequent inference, as we will see in the following analysis.

When conducting KMBGD estimation, I need to choose one covariate and normalize its coefficient to 1. To improve the numerical performance of the method, I choose to normalize the coefficient of the continuous variable Distance. According to khan2022, the covariate whose coefficient is normalized must have positive impacts on the conditional probability. Since a larger geographic distance is generally associated with higher trading costs, the covariate Distance has negative impacts on the conditional probability of the presence of trades between two countries\footnote{When I apply Logit or Probit to model ((ref)), the estimated coefficient of Distance is significantly negative.}. In this case, I use the negative value of (logarithm of) Distance instead of the original variable when performing iteration. So any covarite whose coefficient is estimated to be positive can be explained to have positive impacts on the conditional probability.

When estimating the model, I leave out as few fixed effects as possible to ensure that my covariate matrix is nonsingular. When conducting iteration for KMBGD, I choose learning rate $\delta_k = 1$ for all $k$ and subsample size $B = 1000$. When constructing kernel estimator, I choose sixth-order Epanechnikov kernel function, and the bandwidth $h_n$ is chosen as $h_n = c_k\cdot h_n^{-1/10}$, where $c_k = \text{std}\left(z_{i,k}\right)$ and $z_{i,k} = X_{0,i}+\mathbf{X}_i^{\mathrm{T}}\boldsymbol{\beta}_k$. The initial guess of the parameter is fixed at the Probit estimator. I update the estimator 500000 times and use the last 50000 updated estimators to construct the AKMBGD estimator.

Apart from KMBGD estimator, I also consider the full-sample-based SBGD estimator prposed in KLTY. To construct such estimator, I choose learning rate $\delta_k=1$ for all $k$ and the order of sieves $q=25$. The basis functions are the same as in KLTY. The initial guess is also fixed at the Probit estimator. The stopping rule is $\max_{1\leq j\leq p} \left|\beta_{j,k+1} - \beta_{j,k}\right|<10^{-6}$, where $\beta_{j,k+1}$ is the $j$-argument of $\boldsymbol{\beta}_k$ or the number of updates exceeds 500000. To further provide some comparisons between parametric and semiparametric estimation, I also consider parametric estimation based on Logit and Probit regression.

The estimation results are reported in (ref). I first compare the computational time of different methods. Obviously, parametric Probit and Logit estimation feature fast computation, which both take around 1 minute. On the other side, the semiparametric estimation based on KMBGD and SBGD take 8.0--9.0 hours, which are all computationally feasible. Comparitively, the subsample-based KMBGD is slightly faster in terms of estimation, and significantly outpeforms the SBGD method in terms of the operation time of inference.

Next I compare the estimation results of different estimation methods. I find that, first of all, the Logit estimator differs significantly from the Probit estimator for some coefficients. For example, the estimated coefficient of Island using Probit is 3.600 with standard deviation 0.100. So under Probit estimation, the 0.95 confidence interval for the coefficient of Island is $[3.404, 3.796]$, which does not include the Logit estimator $3.400$. This implies that if the random shock in the binary choice model actually has a Logistic distribution instead of standard normal distribution, then there is a high probability ($\geq 50\%$) that the confidence interval based on Probit does not include the unknown true parameter. Indeed, the semiparametric estimation results strongly favor such possibility. In particular, it can be seen that the KMBGD estimator is quite close to the Logit estimator. For example, for the coefficient of Island, the Logit estimator is $3.400$ and the KMBGD estimator is $3.395$, which almost coincide with each other. Similar patterns can also be seen from the estimation results of other coefficients. I further compare the SBGD estimator with both Probit and Logit estimators. I also find that comparatively, the SBGD estimator is closer to the Logit estimator. The above result highlights the potential of model misspecification of Probit estimation and motivates the use the semiparametric estimation.

I finally investigate convergence of KMBGD estimator. I plot the KMBGD estimation results (without WTO (both) and WTO (none)) of the first 9 covariates produced during 500000 iterations in (ref). It can be seen that different coefficients exhibit different converging behaviors. For example, for the coefficient of FTA, although the starting point of iteration (which is Probit estimator) deviates a lot from the final estimator, it converges very quickly and starts fluctuating around the AKMBGD estimator after roughly 100000 rounds of updates. While comparitively, the estimators of the coefficients of Island and Landlock converge slowly, which start fluctuating around the final estimators after roughly 300000 and 400000 rounds of updates, respectively.

Concluding Remarks

This paper investigates semiparametric estimation of monotone index models in a large-$n$ environment, where the number of observations is extremely large. I propose a novel subsample- and iteration-based estimation procedure. Essentially, starting from an initial guess of the parameter, in each round of iteration a subsample is randomly drawn and then used to update the parameter based on the gradient of some well-chosen loss function, where the unknown nonparametric component is replaced with its subsample-based kernel estimator. The proposed algorithm essentially generalizes the idea of mini-batch-based algorithms to the semiparametric setup. Compared with the KBGD algorithm proposed in KLTY, the computational speed of the new estimator substantially improves, so can be easily applied when the sample size $n$ is extremely large. I also show that further averaging across the estimators produced during iterations yields a $1/\sqrt{n}$ consistent and asymptotically normally distributed estimator.

As an empirical application of the new method, I revisit the Probit estimation of the presence of trade between countries in helpman2008estimating. Given the large sample size and number of covariates, the computational time of estimation and inference based on KMBGD algorithm is reasonable. I also find that compared with Probit specification, the semiparametric estimation results are more in favor of the Logistic distributed random shock in the binary choice model, which highlights the use of semiparametric estimation in the empirical applications.

Some issues in this paper remain to be addressed in the future studies. For example, similar to ichimura1993semiparametric, I show that a particular sequence of bandwidth satisfying some order conditions guarantees all the theorems. However, in the theorem the bandwidth is assumed to be unchanged across iterations. Obviously, as the updates proceed, the magnitude of the index value also changes, so a bandwidth adjusted to such change in index value in each round of iteration may lead to a better kernel estimator and improve the updating results. Similarly, other tuning parameters such as the learning rate $\delta $ and subsample size $B$ are all assumed to be given, while their optimal choices remain to be studied.

Another potential future research direction is to generalize the noval subsample-based updating techinque to the full-sample-based SBGD algorithm proposed in KLTY. Different from the kernel-based learning approach, the SBGD algorithm relies on the full sample to update the sieve coefficient in each iteration. So it is still unclear whether using subsamples to perform the update will also yield $1/\sqrt{n}$-consistent estimator. However, since the SBGD algorithm runs significantly faster than the KBGD algorithm, developing subsample-based SBGD algorithm may further improve the computational speed, which deserves further study.