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.
121,784 characters
Online Learning in Semiparametric Econometric Models
\title{Online Learning in Semiparametric Econometric Models}
\author{Xiaohong Chen\thanks{Dept. of Economics, Yale University.} \ \ \ Elie Tamer\thanks{Dept. of Economics, Harvard University.} \ \ \ Qingsong Yao\thanks{Dept. of Economics, LSU. }}
\maketitle
\begin{abstract}
\singlespacing
\noindent
Data in modern economic and financial applications often arrive as a stream, requiring models and inference to be updated in real time—yet most semiparametric methods remain batch-based and computationally impractical in large-scale streaming settings.
We develop an online learning framework for semiparametric monotone index models with an unknown monotone link function.
Our approach uses a two-phase learning paradigm. In a warm-start phase, we introduce a new online algorithm for the finite-dimensional parameter that is globally stable, yielding consistent estimation from arbitrary initialization.
In a subsequent rate-optimal phase, we update the finite-dimensional parameter using an orthogonalized score while learning the unknown link via an online sieve method; this phase achieves optimal convergence rates for both components.
The procedure processes only the most recent data batch, making it suitable when data cannot be stored (e.g., memory, privacy, or security constraints), and its resulting parameter trajectories enable online inference such as confidence regions - on parameters including policy-effect analysis with negligible additional computation. Monte Carlo experiments on both simulated and real data show adequate performace especially relative to full sample methods.
\end{abstract}
\noindent
\newpage
\section{Introduction}
Consider a generic monotone index model
\begin{align}
Y = F_0(x_0 + X^{\prime}\theta_0) + \varepsilon, \ \ \ \mathbb E(\varepsilon|x_0, X) = 0.\label{model}
\end{align}
where $Y$ is observed response, $F_0(\cdot )$ is monotonically increasing and is unknown, $(x_0, X^{\prime})^{\prime}$ is $(p + 1)\times 1$ observed regressor, $\theta_0$ is $p\times 1$ true parameter, and $\varepsilon$ is unobserved error term. For identification purpose, we normalize the coefficient of $x_0$ as 1 and absorb the scale into the unknown link $F_0$. Such class of models nests a wide range of econometric models such as linear regression, binary choice, censored, and hazard models \citep{han1987non,sherman1993limiting}. The researcher's primary interest here is to estimate the finite-dimensional parameter $\theta_0$ and uncover the shape of the infinite-dimensional parameter $F_0$. A large literature in econometrics developed methods for inference on this class of models. See for example \cite*{han1987non,powell1989semiparametric,ichimura1993semiparametric,klein1993efficient,neweyetal,ahn2018simple,fan2020rank,khanetal2021} and also \cite{powellhandbook} and references therein. Notably, all the above methods are designed by nature to be \textit{offline} in the sense that the training set is fixed prior to estimation, so the researcher conducts only one-time estimation based on such given sample and no repeated update of the estimation results is needed.
This paper considers a different learning task where data $\{(Y_i, x_{0,i}, X_i)\}_{i=1}^{\infty}$ arrive sequentially, so the sample expands continuously over time. Such learning task is therefore dynamic: in each period, the researcher produces updated estimates (and, when desired, conducts inference) for $\theta_0$ and $F_0$
using only the observations available up to that point. This setting falls under {\it online learning}, which has attracted growing attention as data are increasingly generated and updated in real time---for example, in financial trading or mobile/web applications. While the offline methods can still be applied in this online setting, they have important drawbacks in that they require re-estimating the model using the entire accumulated dataset whenever new observations arrive, which is computationally costly. More importantly, they require storing an ever-growing dataset and repeatedly accessing it in its entirety each time new data are incorporated.
Against this backdrop, the paper’s main contribution is to characterize the statistical properties of the semiparametric econometric models in
(\ref{model}) in this online setting, while explicitly allowing for unknown functions. The parameter of interest may be the finite-dimensional vector $\theta_0$ itself, or a functional of the model—such as an average or marginal effect—that typically depends on $\theta_0$ as well as on unknown objects like $F_0$. In doing so, the paper pushes the existing (and important) semiparametric econometrics toolkit toward accommodating online data environments.
\paragraph{Details and Literature}
To alleviate the computational and data storage burdens, online learning typically uses only the newest data point or a small number of most recent data points to update the estimators, which can be regarded as stochastic approximation (SA) originally proposed by \citet{robbins1951stochastic}. When a gradient is used for update, online learning can be viewed as stochastic gradient descent (SGD, \citet{bottou2007tradeoffs,bottou2010large,bottou2018optimization}), which solves the optimization problem in finite-dimensional spaces or infinite-dimensional spaces such as Hilbert space \citep{chen2002asymptotic,godichon2019lp}. A large body of literature studies the properties of SGD algorithms for loss function with strong convexity \citep{toulis2017asymptotic,jentzen2021strong,sebbouh2021almost}, local strong convexity \citep{godichon2019lp}, or nonconvexity \citep{ghadimi2013stochastic,ghadimi2016accelerated,mertikopoulos2020almost,sebbouh2021almost}. See also the references therein. There is also a small branch of literature focusing on the SGD/online estimation in the nonparametric setup, such as sieve estimation \citep{zhang2022sieve} or kernel density estimation/regression \citep{huang2013recursive}. More recently, SA methods have been applied to analyze econometric models such as GMM models \citep{chen2023sgmm,chen2025slim} under the context of large sample size or online learning.
While online learning has been extensively studied in the machine learning context, it has not been widely applied to semiparametric estimation, especially under the M-estimation setup. This is presumably due to the fact that, under the semiparametric setup, the loss function to be optimized depends both on the finite-dimensional parameter and the infinite-dimensional nonparametric component. Simply optimizing the loss with respect to both components may be ill-posed \citep{chen2007large}, while optimization with respect to the parametric component only is practically infeasible because the loss function and its gradient (with respect to the parametric component) depend on the unknown nonparametric component. One can instead replace $F_0$ with its kernel or sieve estimator $\hat F$ and use the plug-in loss function and its gradient to perform updates. The estimator $\hat F$ is usually a function of $\theta$ (and of the sieve coefficients in the sieve estimation setup; see \citet{shen1994convergence,shen1997methods,chen1998sieve}. This makes the plug-in loss function heavily nonlinear and nonconvex, and hence local optima creates serious problems \citep*{khanetal2021,khan2024inference}.
This paper develops a novel \textit{two-phase} online learning paradigm for models in (\ref{model}). We contribute to the literature on online learning and semiparametric estimation of monotone index models in several ways.
First, we propose a new online algorithm for learning $\theta_0$. The algorithm uses a new score function that is similar to the score of \citet{han1987non}'s MRC loss function. Under mild conditions, we show that the limiting Jacobian matrix associated with the score function is everywhere strictly positive definite, so our algorithm guarantees asymptotic global stability, meaning that the algorithm consistently estimates $\theta_0$ regardless of the choice of the starting point. We provide the almost sure convergence, the law of iterated logarithm, and the limiting distribution of the resulting estimator and its Polyak-Ruppert (PR) averages \citep{ruppert1988efficient,polyak1992acceleration}. Our first phase, ``warm-start phase'' has an attractive global stability property.
After the warm-start learning phase, the learning path is guaranteed to be in a small neighborhood of $\theta_0$. Starting within such small neighborhood, we initiate our second learning phase. The second learning phase simultaneously updates the estimators of both $\theta_0$ and $F_0$. For the update of $\theta_0$, we use an orthogonalized score function to remove the first-order impacts of the (estimated) nuisance parameter. We show that the resulting online estimator of $\theta_0$ is $1/\sqrt{N}$-consistent, where $N$ is the total number of accumulated data points. Law of iterated logarithm, limiting distribution, and functional central limit theorem are also provided for the resulting PR average estimator. On the other side, we propose to use the method of sieves to learn the nonparametric component $F_0$, where the order of the sieve functions increases at some rate as the update proceeds. We provide the asymptotic linear representation for the online sieve coefficient estimators, and the almost sure supreme-norm convergence rate. We show that the sup-norm convergence rate is comparable to that in \citet*{chen2015optimal} and \citet*{belloni2015some}. Notably, our results on the sup-norm rate of sieve online estimation error apply to general sieve-based nonparametric regression with generated regressors. In the second learning phase, the estimators of both $\theta_0$ and $F_0$ achieve the optimal convergence rate. As a result, such learning phase is labeled ``rate-optimal learning phase''.
As an online learning paradigm, the above-described two-phase learning procedure takes the most recent data points as inputs and produces the updated parameters as outputs. The updated parameters across updates create learning trajectories of $\theta_0$ and $F_0$, based on which online inference can be effectively conducted using random scaling \citep*{lee2022fast,chen2023sgmm,chen2025slim}. In particular, random scaling utilizes updated parameters to estimate the long-run variance of the trajectory, based on which confidence band for the true parameter can be constructed. Different from the
conventional semiparametric inference using estimated variance matrix as plug-in, inference based on random scaling takes only the trajectory as input so avoids intensive nonparametric estimation and requires almost no additional computational costs. Apart from online inference, we show that the learning trajectories of $\theta_0$ and $F_0$ can also be used to construct estimators of policy impacts, which can be formulated as functionals of both $\theta_0$ and $F_0$. Online policy evaluations and inference are then studied.
In section \ref{implementation} we provide detailed instructions via algorithms that show how one can implement our procedures for learning $\theta_0$. Also, Section \ref{sec:policy} extends these to procedures to obtain marginal effects. Finally, section \ref{sec:MC} implements an extensize Monte Carlo simulations and shows adequate performaance for our estimators. These are alway applied to real data from a trade application.
The semiparametric online learning paradigm developed in this paper can be extended further. For example, our algorithm can be adjusted to analyze sample selection effect \citep*{heckman1974,abhausmankhan,khan2024inference}. In general, if the observability of the data depends on an unobserved selection process, we can use the proposed algorithm to first learn the selection equation, based on which we can control for the selection bias and learn the outcome equation. Note that the selection and outcome equations are both semiparametric, and the learning of both equations are online. More broadly, this paper pushes the semiparametric analysis towards the online learning framework, which admits both data and model flexibility.
\paragraph{Notations} For any vector $A = (a_1, \cdots, a_n)^{\prime}$, we use $\Vert A\Vert = \sqrt{\sum_{i=1}^n a_i^2}$ to denote its Euclidean norm. For any matrix $A = (a_{ij})_{m\times n}$, we use $\Vert A \Vert = \sqrt{\sum_{i=1}^m\sum_{j=1}^n a_{ij}^2}$ to denote is Frobenius norm. For any function $f(z)$ with domain $\mathcal{Z}$, we use $\Vert f \Vert_{\infty} \equiv \sup_{z\in\mathcal{Z}} |f(z)|$ to denote the supreme norm of $f$. For any positive sequences $\{a_N\}_{N=1}^{\infty}$ and $\{b_N\}_{N=1}^{\infty}$, we write $a_N = o(b_N)$ if $\overline{\lim}_{N\rightarrow \infty} a_N/b_N= 0$, and $a_N = O(b_N)$ if $\overline{\lim}_{N\rightarrow \infty} a_N/b_N<\infty$. If both $a_N = O(b_N)$ and $b_N = O(a_N)$, we write $a_N\asymp
b_N$. For any symmetric matrix $A$, we use $\overline{\lambda}(A)$ and $\underline{\lambda}(A)$ to denote the largest and smallest eigenvalue of $A$.
\section{Learning Phase I: Warm-Start Learning} \label{section2}
This section introduces a new online algorithm for semiparametrically estimating $\theta_0$. The proposed algorithm features global stability, so can be used to quickly locate a small neighborhood of $\theta_0$ regardless of the starting point.
Suppose that at the beginning (the 0-th round) of the update, the initial guess is given by $\hat\theta_0$. In many cases $\hat\theta_0$ can be simply chosen as the origin point $\boldsymbol{0}_p$. For any $k\geq 1$, the researcher observes a batch of $B$ \textit{i.i.d.} realizations of $(x_0, X, Y)$ from model (\ref{model}) denoted as $(x_{0,i,k}, X_{i,k},Y_{i,k})$, $i = 1, 2, \cdots, B$. To ease exposition, we denote \begin{equation}\label{data}
W_k = \{(x_{0,i,k}, X_{i,k}, Y_{i,k})\}_{i=1}^B
\end{equation} as the data observed in the $k$-th round. In the rest of the paper, a maintained assumption over the observations is that
$\{W_k\}_{k=1}^{\infty}$
are \textit{i.i.d.}.
In the $k$-th round, the task for the researcher is to update $\hat\theta_{k-1}$ to $\hat\theta_{k}$ using $W_k$. Denote $z_{i,k}(\theta) \equiv x_{0,i,k} + X_{i,k}^{\prime}\theta$ for $i=1,\cdots, B$ and any $\theta$, we propose the following online updating algorithm:
\begin{equation}\label{first-step algorithm}
\hat \theta_{k} = \hat \theta_{k-1} + \frac{\gamma_k}{h_k\cdot B(B-1)}\cdot\sum_{i_1\neq i_2}^B \mathcal{K}\left(\frac{z_{i_1,k}(\hat\theta_{k-1}) - z_{i_2,k}(\hat\theta_{k-1})}{h_k}\right)\left(Y_{i_1,k} - Y_{i_2,k}\right)\left(X_{i_1,k} - X_{i_2,k}\right),
\end{equation}
where $\gamma_k>0$ is the learning rate in the $k$-th round of update, $\mathcal{K}$ is a kernel function and $h_k$ is the bandwidth parameter that depends on $k$.
Next we provide some comments on update (\ref{first-step algorithm}).
\begin{remark}
The updating algorithm (\ref{first-step algorithm}) is closely related to \citet{han1987non}'s MRC estimator defined as
\begin{equation}
\hat\theta_{MRC} = \arg\max_{\theta}\frac{1}{n^2}\sum_{1\leq i,j\leq n } I(Y_i> Y_j)I(x_{0,i}+X_i^{\prime}\theta> x_{0,j}+X_j^{\prime}\theta),
\end{equation}
where $I(\cdot)$ is indicator function, and $\{(x_{0,i}, X_i, Y_i)\}_{i=1}^n$ are observed i.i.d. data points. Note that the above MRC loss can be used for semiparametric estimation and does not involve the nuisance parameter $F_0$. However, the critical disadvantage of the MRC loss function is that it is not differentiable with respect to $\theta$, which makes optimization numerically difficult. Consider a smoothed version of the MRC loss function similar to \citet{horowitz1992smoothed} given by
\begin{equation}
L_n(\theta) = \frac{1}{n^2}\sum_{1\leq i,j\leq n} I(Y_i>Y_j)K\left(\frac{x_{0,i}+X_i^{\prime}\theta - x_{0,j} - X_j^{\prime}\theta}{h_n}\right),
\end{equation}
where $K$ is a smooth increasing function satisfying $\lim_{z\rightarrow -\infty} K(z) =0$ and $\lim_{z\rightarrow \infty} K(z) =1$. Note that the derivative of the smoothed MRC loss with respect to $\theta$ is
\[
\frac{\partial L_n(\theta)}{\partial \theta} = \frac{1}{n^2h_n}\sum_{1\leq i,j\leq n} I(Y_i>Y_j)\left(X_i - X_j\right)\mathcal{K}\left(\frac{x_{0,i} + X_i^{\prime}\theta - x_{0,j} - X_j^{\prime}\theta}{h_n}\right)
\]
where $\mathcal{K}(\cdot) = \partial K(\cdot)$. So our update is similar to the mini-batch gradient descent version of the smoothed MRC loss.
However our algorithm (\ref{first-step algorithm}) is still different from the MRC loss derivative in that we replace the indicator $I(Y_{i,k} > Y_{j,k})$ with the difference between the responses $Y_{i,k} -Y_{j,k}$. Note that when the responses are binary,
\begin{align*}
& \left(Y_{i,k} - Y_{j,k}\right)\left(X_{i,k} - X_{j,k}\right)
= I\left(Y_{i,k} > Y_{j,k}\right)\left(X_{i,k} - X_{j,k}\right)
+ I\left(Y_{j,k} >Y_{i,k}\right)\left(X_{j,k} - X_{i,k}\right).
\end{align*}
So our updating algorithm degenerates to the SGD update with smoothed MRC loss function. However, when the response is continuous, our algorithm is essentially different from the one directly derived from the MRC loss. In the following we show that such modification guarantees global contraction mapping of our algorithm.
Also note that the MRC derivative depends on the bandwidth parameter $h_n$, which is a function of sample size $n$. In the offline learning setup, the researcher acquires all the data points before estimation, so the sample size $n$ and the bandwidth parameter is well-defined. However, in the online learning setup, data points are received sequentially, so the sample size is essentially changing throughout time, leaving the bandwidth undefined. To accommodate such fact, we follow \citet*{huang2013recursive} and set the bandwidth in the $k$-th round of update to be $h_k$. In other words, the bandwidth of the one-time update depends on the number of data points collected up to that updating point.
\end{remark}
\begin{comment}
Even though both target a monotone index model, the iterative update in Equation (2.3) is distinct from the gradient decent offline approache used in \citep*{khan2024inference}. The latter's approach employs a \textit{pointwise} least-squares formulation, $\sum_{i} (Y_i - F(X_i'\theta))^2$, which requires simultaneous estimation of the unknown link function $F$ and the parameter $\theta$. While statistically efficient, this joint estimation results in a non-convex optimization problem that is sensitive to initialization and computationally intensive ($O(n^2)$ or $O(n)$ per step depending on the implementation).
In contrast, Equation (2.3) utilizes a \textit{pairwise} rank-based update derived from the Maximum Rank Correlation (MRC) objective, $\sum_{i \neq j} (Y_i - Y_j) \cdot \text{sgn}((X_i - X_j)'\theta)$. By focusing solely on the relative ordering of the data, this method eliminates the need to estimate $F$ directly, ensuring the objective is globally concave (after smoothing).
\end{comment}
In most of the cases, the online estimator $\hat\theta_k$ does not achieve the best convergence rate. So we also consider Polyak-Ruppert averages of $\hat\theta_k$ defined as
\begin{align}\label{PR average}
\bar{\theta}_N = \frac{1}{N}\sum_{k=1}^N\hat\theta_{k}.
\end{align}
Typically, $\bar{\theta}_N$ is the sample average of the first $N$ iterates obtained based on (\ref{first-step algorithm}). Note that evaluating $\bar{\theta}_N$ based on algorithm (\ref{PR average}) requires storing the full trajectory of updates. To save memory, we consider the following online update for evaluating $\bar{\theta}_N$
\begin{equation}\label{average_algorithm}
\bar{\theta}_N = \frac{N-1}{N}\bar{\theta}_{N-1}+ \frac{1}{N}\hat\theta_N.
\end{equation}
In this case, only the most recent PR average estimator needs to be stored in the memory.
\subsection{Preliminaries and Intuitive Discussion}\label{section2.1}
Define \begin{equation}\label{empirical score}\Phi_k\left(\theta, W_k\right) =\frac{1}{h_k\cdot B(B-1)}\cdot\sum_{i_1\neq i_2}^B \mathcal{K}\left(\frac{z_{i_1,k}(\theta) - z_{i_2,k}(\theta)}{h_k}\right)\left(Y_{i_1,k} - Y_{i_2,k}\right)\left(X_{i_1,k} - X_{i_2,k}\right),\end{equation}
\begin{equation}\label{theoretical score}
\Phi_k\left(\theta\right) = \mathbb E\left[\Phi_k\left(\theta, W_k\right)\right],
\end{equation}
and
\begin{equation}\label{limit score}
\Phi\left(\theta\right) = \lim_{k\rightarrow\infty}\Phi_k\left(\theta\right).
\end{equation}
Denote the joint density of $x_0$ and $X$ as $f(x_0,X)$. For any $\theta\in\mathbb{R}^p$, define $\Delta\theta \equiv \theta - \theta_0$. We impose the following conditions.
\begin{condition}\label{condition1}
(i) $F_0(\cdot)$ has up to $(s+1)$-th uniformly bounded derivative; (ii) $f(x_0,X)$ has up to $(s+1)$-th uniformly bounded partial derivative with respect to $x_0$; (iii) there exist nonnegative functions $g_0(X), \cdots, g_{s+1}(X)$ such that for $0\leq j \leq s+1$, there hold $\sup_{x_0\in\mathbb{R}}|\partial^jf(x_0, X)/\partial x_0^j|\leq g_j(X)$ and for all $j\geq 0$, $\int \left\Vert X_1 \right\Vert^4 g_j\left(X_1\right) dX_1<\infty$; moreover, $\mathbb{E}\Vert X\Vert^4<\infty$; (iv) define $\sigma^2(x_0, X) = \mathbb E(\varepsilon^2|x_0, X)$, there holds $\sup_{(x_0, X^{\prime})^{\prime}\in\mathbb{R}^{p+1}}\sigma^2(x_0, X)<\infty,$ moreover, $\sigma^{2}(x_0, X)$ has uniformly bounded derivative with respect to $x_0$.
\end{condition}
\begin{condition}\label{condition2}
The kernel function $\mathcal{K}$ satisfies: (i)
$\int_{-\infty}^{\infty} \mathcal{K}(t)dt =1, \int_{-\infty}^{\infty} t^{j}\mathcal{K}(t)dt = 0$ for $1\leq j\leq s-1,$ and $\int_{-\infty}^{\infty} |\mathcal{K}(t)t^{j}|dt <\infty$ for any $ j\geq s$; (ii) $\int_{-\infty}^{\infty}|\mathcal{K}^2(t)t^j|dt<\infty$ for any $j\geq 0$; (iii) $\int_{-\infty}^{\infty}|(\partial\mathcal{K}(t))^2t^j|dt<\infty$ for any $j\geq 0$.
\end{condition}
Given \autoref{condition1} and \autoref{condition2}, we have the following lemma that describes the properties of $\Phi_k(\theta,W_k), \Phi_k(\theta),$ and $\Phi(\theta)$.
\begin{lemma}\label{lemma1}
Let \autoref{condition1} and \autoref{condition2} hold, for any choice of $\{h_k\}_{k=1}^{\infty}$ with $h_k\downarrow 0$, we have that
\begin{enumerate}
\item There holds \begin{align*}
\Phi(\theta) & = \int \left[ F_0\left(z - X_1^{\prime}\Delta\theta\right)-F_0\left(z - X_2^{\prime}\Delta\theta\right)\right]\left(X_1 - X_2\right)\times\\
& \ \ \ \ \ \ \ \ \ \ \ \ \ \ f\left(z - X_1^{\prime}\theta, X_1\right)f\left(z - X_2^{\prime}\theta, X_2\right)dzdX_1dX_2,
\end{align*}
which is well defined for any $\theta\in\mathbb{R}^p$;
\item There exists a positive constant $C_{\Phi, 1}$ such that there holds \[\sup_{\theta\in\mathbb{R}^p}\left\Vert\Phi_k(\theta) - \Phi(\theta)\right\Vert\leq C_{\Phi, 1}h_k^s;\]
\item There holds
\[\mathbb E\left[\left(\Phi_k(\theta, W_k) - \Phi_k(\theta)\right)\left(\Phi_k(\theta, W_k) - \Phi_k(\theta)\right)^{\prime} \right] = \frac{2\mathcal{V}_{\mathcal{K}}(\theta)}{B(B-1)h_k} + O(1),\]
where
\begin{align*}
\mathcal{V}_{\mathcal{K}}(\theta) = & \int \mathcal{K}^2(t)dt \int \left(\sigma^2(z - X_1^{\prime}\theta, X_1) + \sigma^2(z - X_2^{\prime}\theta, X_2) + \left(F_0(z - X_1^{\prime}\Delta\theta) - F_0(z - X_2^{\prime}\Delta\theta)\right)^2\right) \times \\
& \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (X_1-X_2)(X_1-X_2)^{\prime} f(z- X_1^{\prime}\theta, X_1)f(z- X_2^{\prime}\theta, X_2)dzdX_1dX_2;
\end{align*}
So there exists a positive constant $C_{\Phi, 2}$ such that there holds \[\sup_{\theta\in\mathbb{R}^p}\mathbb E\left\Vert\Phi_k(\theta, W_k) - \Phi_k(\theta)\right\Vert^2\leq \frac{C_{\Phi, 2}}{B^2h_k};\]
\item There exists a positive constant $C_{\Phi, 3}$ such that \[\sup_{\theta\in\mathbb{R}^p} \mathbb E\left\Vert\partial_{\theta}\Phi_k(\theta, W_k) - \partial_{\theta}\Phi_k(\theta)\right\Vert^2\leq \frac{C_{\Phi, 3}}{B^2h_k^3}.\]
\end{enumerate}
\end{lemma}
Given \autoref{lemma1}, we can briefly discuss the intuition of our algorithm. Note that (\ref{first-step algorithm}) leads to
\begin{align}\label{intuition}
\hat\theta_k & = \hat\theta_{k-1} + \gamma_k \cdot \Phi_k(\hat\theta_{k-1},W_k)\nonumber \\
& = \underset{(I)}{\underbrace{\hat\theta_{k-1} + \gamma_k \cdot \Phi(\hat\theta_{k-1})}} - \underset{(II)}{\underbrace{ \gamma_k \cdot \left( \Phi(\hat\theta_{k-1}) - \Phi_k(\hat\theta_{k-1})\right)}} - \underset{(III)}{\underbrace{\gamma_k\cdot \left( \Phi_k(\hat\theta_{k-1}) - \Phi_k(\hat\theta_{k-1}, W_k)\right)}}.
\end{align}
In the above update, term (III) has zero expectation conditioned on observations up to period $k-1$, and term (II) will vanish as $k$ increases if we choose $h_k\downarrow 0$ according to \autoref{lemma1}. So we can regard the update (\ref{first-step algorithm}) as being mainly driven by term (I). Define
\begin{align}
\mathbb{H}(\theta, \tau) \equiv & \int \partial F_0(z - X_2^{\prime}\Delta\theta - \tau(X_1 - X_2)^{\prime}\Delta\theta)f(z - X_1^{\prime}\theta, X_1)f(z-X_2^{\prime}\theta, X_2)\nonumber \\
& \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \times\left(X_1 - X_2\right)\left(X_1 - X_2\right)^{\prime}dzdX_1dX_2,
\end{align}
then for any $\theta$, we obviously have that \begin{equation}\Phi(\theta) =- \left[\int_{0}^{1} \mathbb{H}(\theta,\tau)d\tau \right]\Delta\theta\end{equation} and \begin{equation}\underline{\lambda}\left(\int_{0}^1\mathbb{H}(\theta, \tau)d\tau\right)\geq 0.\end{equation}
If additional conditions are imposed on the data generating process so that $\underline{\lambda}(\int_{0}^1\mathbb{H}(\theta, \tau)d\tau)> 0$ for any $\theta$, (I) is a contraction mapping and the same applies to (\ref{first-step algorithm}) as $k\rightarrow \infty$. Consequently, our method guarantees global stability, meaning that we do not require any assumptions over the initial starting point of the update.
\begin{remark}
[Contraction Intuition] To better understand the global stability property presented above, we examine the deterministic component of the update rule. By ignoring the noise terms (II) and (III), the evolution of the estimator $\hat{\theta}_{k}$ is driven principally by the population score function $\Phi(\hat\theta_{k-1})$, so we have $\hat{\theta}_{k} \approx \hat{\theta}_{k-1} + \gamma_k \Phi(\hat{\theta}_{k-1})$. Recall that $\Delta \hat{\theta}_{k} = \hat{\theta}_{k} - \theta_0$, we obtain a recursive relationship for the estimation error
\begin{equation*}
\Delta \hat{\theta}_{k} \approx \left(\mathbb I_p - \gamma_k \bar{\mathbb{H}}_{k-1} \right) \Delta \hat{\theta}_{k-1},
\end{equation*}
where $\bar{\mathbb{H}}_{k-1} = \int_{0}^{1} \mathbb{H}(\hat{\theta}_{k-1}, \tau) \, d\tau$. By \autoref{lemma2} introduced in the next section, $\bar{\mathbb{H}}_{k-1}$ is strictly positive definite with lower-bounded eigenvalues over any compact ball of $\theta$. Consequently, for a sufficiently small learning rate $\gamma_k$, the eigenvalues of the operator $(\mathbb I_p - \gamma_k \bar{\mathbb{H}}_{k-1})$ lie strictly between 0 and 1. This spectral property ensures that the norm of the error vector contracts in every step ($\Vert \Delta \hat{\theta}_{k}\Vert \leq (1 - c\gamma_k) \Vert\Delta \hat{\theta}_{k-1}\Vert $ for some small $c>0$ with $c\gamma_k<1$), guaranteeing convergence to the true parameter $\theta_0$ regardless of the initialization as long as $\sum_{k=1}^{\infty}\gamma_k = \infty$.
\end{remark}
\subsection{Statistical Properties of $\hat\theta_k$ and $\bar\theta_N$}
This section formally studies the statistical properties of our proposed algorithm (\ref{first-step algorithm}) and (\ref{PR average}). We first impose some further conditions.
\begin{condition}\label{condition3}
There exist $-\infty<\underline{z}<\bar{z}<\infty$ and $r_X>0$ such that for any $z\in [\underline{z},\bar{z}]$ and $X$ with $\Vert X\Vert \leq r_X$, there holds $\partial F_0(z)\geq \underline{c}_F$ and $f(z, X)\geq \underline{c}_f$, where $\underline{c}_F$ and $\underline{c}_f$ are two positive constants.
\end{condition}
\begin{remark}
Note that $X$ can be normalized to have zero expectation, in which case \autoref{condition3} requires that $X$ has nonvanishing density around its expectation. Moreover, \autoref{condition3} also requires that $x_0$ has sufficiently large support, condition that is required for point identification.
\end{remark}
As we have pointed out before, for any $\theta$ and $\tau$, $\mathbb{H}(\theta,\tau)$ is positive semi-definite. The next lemma provides a more precise bound for the smallest and largest eigenvalues of $\mathbb{H}(\theta,\tau)$ based on the input $\theta$. The result is fundamental for the contraction mapping property of our algorithm.
\begin{lemma}\label{lemma2}
If \autoref{condition1} holds, then there exists $C_{\Phi, 4}>0$ such that
\[
\sup_{\theta\in\mathbb{R}^p}\sup_{\tau\in[0,1]}\overline{\lambda}\left(\mathbb{H}(\theta,\tau)\right)\leq C_{\Phi, 4}.
\]
If \autoref{condition3} holds, then
\[
\inf_{\tau\in[0,1]}\underline{\lambda}\left(\mathbb{H}(\theta,\tau)\right)\geq c_{\mathbb{H}} \cdot \left(\frac{\bar{z} - \underline{z}}{12(1 + \Vert\Delta\theta\Vert + \Vert\theta_0\Vert)}\wedge r_X\right)^{2p+2},
\]
where $c_{\mathbb{H}}$ is a positive constant that does not depend on the choice of $\theta$.
\end{lemma}
The lower bound of the smallest eigenvalue as a function of $\theta$ in \autoref{lemma2} is crucial because it guarantees that term (I) in (\ref{intuition}) is everywhere a strict contraction mapping. Note that such result also guarantees point identification of $\theta_0$.
Based on \autoref{lemma2}, we are ready to state the first theorem regarding the convergence of the online estimator $\hat\theta_k$.
\begin{theorem}\label{theorem1}
Let \autoref{condition1}--\autoref{condition3} hold. If further $\sup_{k\geq 1} \gamma_k C_{\Phi, 4}<1$, $\sum_{k=1}^{\infty}\gamma_k = \infty$, $\sum_{k=1}^{\infty}\gamma_k^2h_k^{-1}<\infty$ and $\sum_{k=1}^{\infty}h_k^{2s}<\infty$, we have that $\hat\theta_k \longrightarrow_{a.s.} \theta_0$.
\end{theorem}
\autoref{theorem1} states the a.s. convergence of our semiparametric online estimator. So for almost all paths of the data stream, our algorithm will lead to consistent estimator for the unknown parameter $\theta_0$ as long as the number of updates is sufficiently large.
Given \autoref{theorem1}, we can further show the convergence rate of $\Delta\hat\theta_k$ as well as its asymptotic distribution. To start with, we assume that the choices of learning rate $\gamma_k$ and bandwidth parameter $h_k$ satisfy the following condition.
\begin{condition}\label{condition4}
$\gamma_k = \gamma_0 k ^{-\alpha_{\gamma}}$ and $h_k = h_0 k ^{-\alpha_{h}}$ with $\gamma_0,h_0>0$, $\sup_{k\geq 1} \gamma_k C_{\Phi, 4}<1$, $0<\alpha_{\gamma}<1$, $\alpha_h>0$, and $2\alpha_{\gamma} - 3\alpha_h >1$.
\end{condition}
Under \autoref{condition4}, all the requirements on $\gamma_k$ and $h_k$ in \autoref{theorem1} are satisfied. Given the above conditions, we are able to illustrate the sharp rate of $\Vert \Delta\hat\theta_k\Vert$. We first introduce some additional notations. Let $\mathbb{H}_0$ be eigendecomposed as \begin{equation}\mathbb{H}_0 = \mathcal{P}_{\mathbb{H}_0}^{\prime}\varLambda_{\mathbb{H}_0}\mathcal{P}_{\mathbb{H}_0},\end{equation} where $\varLambda_{\mathbb{H}_0}$ is the diagonal matrix consisting of eigenvalues of $\mathbb{H}_0$, and $\mathcal{P}_{\mathbb{H}_0}\mathcal{P}_{\mathbb{H}_0}^{\prime} = \mathcal{P}_{\mathbb{H}_0}^{\prime}\mathcal{P}_{\mathbb{H}_0} = \mathbb{I}_p$. Define $
\mathbb{S}_N^2$ to be a $p\times p$ diagonal matrix with $j$-th diagonal argument being given by
\[
\left[\mathbb{S}_N^2\right](j,j) = \sum_{k=2}^N \frac{ [\mathcal{P}_{\mathbb{H}_0}\mathcal{V}_{\mathcal{K},0}\mathcal{P}_{\mathbb{H}_0}^{\prime}](j,j)\gamma_0^2h_0^{-1} k^{-2\alpha_{\gamma}+ \alpha_{h}}}{B(B-1)\left[\prod_{l=2}^k (1-\gamma_0\lambda_{\mathbb{H}_0,j}l^{-\alpha_{\gamma}})\right]^2}.
\]
where $[\mathcal{P}_{\mathbb{H}_0}\mathcal{V}_{\mathcal{K},0}\mathcal{P}_{\mathbb{H}_0}^{\prime}](j,j)$ is the $j$-th diagonal element of $\mathcal{P}_{\mathbb{H}_0}\mathcal{V}_{\mathcal{K},0}\mathcal{P}_{\mathbb{H}_0}^{\prime}$ with $\mathcal{V}_{\mathcal{K},0} = \mathcal{V}_{\mathcal{K}}(\theta_0)$, and $\lambda_{\mathbb{H}_0,j}$ is the $j$-th diagonal element of $\varLambda_{\mathbb{H}_0}$. Finally, define $\log(\log(\mathbb{S}_N^2))$ to be the diagonal matrix consisting of $\log(\log(\cdot))$ transformation of diagonal elements of $\mathbb{S}_N^2$.
We have the following theorem.
\begin{theorem}\label{theorem2}
Let \autoref{condition1}--\autoref{condition4} hold, we have that
\begin{equation}\label{lil_theta}
\underline{\overline{\lim}}_{N\rightarrow \infty} \left(2\mathbb{S}^2_N\log(\log(\mathbb{S}^2_N))\right)^{-\frac{1}{2}}\mathcal{P}_{\mathbb{H}_0}^{\prime} \prod_{k=2}^N (\mathbb{I}_p- \gamma_k \mathbb{H}_0)^{-1} \Delta\hat\theta_N = \pm \boldsymbol{1}_p, \ a.s.
\end{equation}
where $\boldsymbol{1}_p$ is a $p\times 1$ vector whose arguments are all 1,
and
\begin{equation}\label{normal_theta}
\left(\mathbb{S}_N^2\right)^{-\frac{1}{2}}\mathcal{P}_{\mathbb{H}_0}^{\prime} \prod_{k=2}^N (\mathbb{I}_p - \gamma_k \mathbb{H}_0)^{-1} \Delta\hat\theta_N \longrightarrow_d \mathbb N(0, \mathbb I_p).
\end{equation}
\end{theorem}
\autoref{theorem2} provides the sharp rate for our online estimator as well as its asymptotic distribution. Using \autoref{lemma4} provided in Section \ref{section2.3} of Appendix, we can verify that \[\left[\mathbb{S}_N^2\right](j) \asymp N^{-\alpha_{\gamma} + \alpha_h }\exp(2\gamma_0\lambda_{\mathbb{H}_0, j}N^{1-\alpha_{\gamma}}),\] and \[\left[\mathcal{P}_{\mathbb{H}_0}^{\prime} \prod_{k=2}^N (\mathbb{I}_p- \gamma_k \mathbb{H}_0)^{-1}\mathcal{P}_{
\mathbb{H}_0}\right](j,j)\asymp \exp(\gamma_0\lambda_{\mathbb{H}_0,j}N^{1-\alpha_{\gamma}}),\]
so (\ref{lil_theta}) implies that $\Vert \Delta\hat\theta_k\Vert = O(\sqrt{N^{-\alpha_{\gamma} + \alpha_{h} }\log(\log(N))})$ a.s.. Since \autoref{theorem2} is fundamental in proving the following theorems and is of interest on its own side, we provide detailed development in Section \ref{section2.3}.
The convergence rate of $\Delta\hat\theta_k$ depends on the choice of $\alpha_{\gamma}$ and hence the learning rate $\gamma_k$, and is not optimal in the sense that it can be further improved by PR average defined in (\ref{average_algorithm}). We have the following theorem.
\begin{theorem}\label{theorem3}Let \autoref{condition1}--\autoref{condition4} hold, we have that
\begin{equation}
\underline{\overline{\lim}}_{N\rightarrow \infty} \left( \frac{4\mathcal{V}_{\mathcal{K},0}N^{1+\alpha_h}\log(\log(N))}{B(B-1)h_0(1+\alpha_{h})}\right)^{-\frac{1}{2}}\mathbb{H}_0N\Delta\bar\theta_N = \pm \boldsymbol{1}_p, \ a.s. ,
\end{equation}
and
\begin{equation}
\left( \frac{2N^{1+\alpha_h}\mathcal{V}_{\mathcal{K},0}}{B(B-1)h_0(1+\alpha_{h})}\right)^{-\frac{1}{2}}\mathbb{H}_0N\Delta\bar\theta_N\longrightarrow_{d} \mathbb N(0, \mathbb I_p).
\end{equation}
\end{theorem}
Given the choice of $\alpha_{\gamma}$, the convergence rate of PR average is improved to $\sqrt{N^{-1 +\alpha_h}\log(\log(N))}$. When $\alpha_{\gamma}$ is chosen satisfying \autoref{condition4}, $\bar\theta_N $ always converges at a faster rate compared with $\hat\theta_N$.
However, both $\hat\theta_N$ and $\bar\theta_N$ converge at slower rates compared with the parametric online estimator and online average estimator, whose convergence rates are $\sqrt{N^{-\alpha_{\gamma}}\log(\log(N))}$ and $\sqrt{N^{-1}\log(\log(N))}$, respectively. This is basically due to the nature of the score function we use to perform the update in (\ref{first-step algorithm}). As we can see from \autoref{lemma1}, the variance of the score function $\Phi_k(\theta, W_k)$ is of order $k^{ \alpha_h }$, which diverges as $k$ increases. The diverging variance makes our update more volatile compared with parametric online learning, which results in a slower learning rate\footnote{For example, suppose we choose a small bandwidth decay to maintain smoothing, such as $\alpha_h = .1$. To satisfy the constraint $2\alpha_\gamma -3\alpha_h>1$, we need $2\alpha_\gamma >1.3$ or $\alpha_\gamma > .65.$ Assuming $\alpha_\gamma = .7$ yields a convergence rate for $\hat \theta_N$ of $N^{-1/3}$ while that of $\bar \theta_N$ of $N^{-.45}.$}.
Leveraging the global stability property, we recommend using the algorithm in this section to rapidly identify a small neighborhood around the true parameter. We use this procedure as the first stage of our online learning algorithm, which we refer to as the warm-start learning phase. Building on the output of this initial stage, Section \ref{section4} introduces a second learning phase that yields a $1/ \sqrt N$
-consistent estimator of $\theta_0.$
\section{Online Learning of $F_0$}\label{section3}
Before introducing the second learning phase, we briefly step aside to present preliminaries for online learning of the nonparametric component $F_0$. This learning procedure will be incorporated into the second phase alongside the estimation of the finite-dimensional parameter $\theta_0.$
\subsection{Problem Setup and Algorithm}
This section studies the question of inference on the unknown nonparametric component $F_0$. In addition to deriving conditions under which we can learn $F_0$ in an online regime, knowing $F_0$ allows also us to also get interesting parameters such as the marginal effect $\partial F_0(X^{\prime}\theta_0)\theta_0$ or the average marginal effect $\mathbb{E}(\partial F_0(X^{\prime}\theta_0))\theta_0$. Motivated by these examples, this section studies the online estimation of the unknown link function based on the \textit{method of sieves}.
Throughout the remainder of the paper, we will use a sequence of known basis functions $\{\psi_j(\cdot)\}_{j=1}^{\infty}$ to approximate the unknown function $F_0 $. The properties of the sieve functions will be detailed in the following conditions. For any positive integer $J$,
define the sieve vector consisting of the first $J$ sieve functions as \begin{equation}\Psi_J(z)= \left(
\psi_1(z), \psi_{2}(z),\cdots, \psi_J(z)\right)^{\prime}.\end{equation} Define $z_0 = x_0 + X^{\prime}\theta_0$ and $\Gamma_J = \mathbb{E}\left[\Psi_J(z_{0})\Psi_J(z_{0})^{\prime}\right]$. Define $\mathcal{S}_J$ as the sieve space spanned by the first $J$ sieve functions, that is, \begin{equation}\mathcal{S}_J = \left\{\sum_{j=1}^J b_j \psi_j(\cdot): b_1, \cdots, b_J\in \mathbb{R}\right\}.\end{equation} Following \citet{chen2015optimal}, we define the $L_2$-projection of function $F_0$ onto $\mathcal{S}_J$ as
\begin{align}\label{L2-projection}
\mathbb{P}_J(F_0)(z) = \Psi_J(z)^{\prime} \Gamma_{J}^{-1}\mathbb{E}\left[ \Psi_J( z_0)F_0(z_0)\right]\equiv \Psi_J(z)^{\prime}\mathcal{B}_{J,0},
\end{align}
where $\mathcal{B}_{J,0} \equiv \Gamma_{J}^{-1}\mathbb{E}\left[ \Psi_J( z_0)F_0(z_0)\right]$ is the pseudo true sieve coefficient vector. Then
for any function $g\in\mathcal{S}_J$, we have $\Vert \mathbb{P}_J(F_0) - F_0\Vert_{\infty} = \Vert \mathbb{P}_J(F_0 - g) - (F_0- g)\Vert_{\infty}$, so the sieve approximation error can be bounded by
\begin{equation}
\Vert \mathbb{P}_J(F_0) - F_0\Vert_{\infty}\leq \inf_{g\in\mathcal{S}_J}\Vert F_0 - g\Vert _{\infty}\left( 1+\frac{\Vert \mathbb{P}_J(F_0 - g)\Vert_{\infty}}{\Vert F_0 - g\Vert _{\infty}}\right).
\end{equation}
We impose the following condition regarding the supreme norm of operator $\mathbb{P}_J(\cdot)$.
\begin{condition}\label{condition5} (i)
There holds
\[
\overline{\lim}_{J\rightarrow \infty} \sup_{g\in \mathcal{S}_J: 0<\Vert F_0 - g \Vert _{\infty}<\infty} \Vert F_0 - g\Vert _{\infty}^{-1}\Vert \mathbb{P}_J(F_0 - g)\Vert_{\infty} <\infty;
\]
(ii) There exist positive constants $C$ such that for each $J$, there exists $g_J^*\in\mathcal{S}_J$ such that
$\Vert F_0 - g_J^*\Vert _{\infty}\leq CJ^{-s};$ (iii) $\Vert \partial F_0 - \partial \mathbb{P}_J(F_0)\Vert _{\infty}\leq CJ^{-s_1}$ for some $0<s_1\leq s$.
\end{condition}
\begin{remark}
\autoref{condition5}(i) bounds the supreme norm of the projection $\mathbb{P}_J$ uniformly with respect to $J$ for the unknown function $F_0$. Under such condition, we easily have that $\Vert \mathbb{P}_J(F_0) - F_0\Vert_{\infty}\leq C\inf_{g\in\mathcal{S}_J}\Vert F_0 - g\Vert _{\infty}$ for some $C$ and any $J$. Then if \autoref{condition5}(ii) further holds, we have that $\Vert \mathbb{P}_J(F_0) - F_0\Vert_{\infty}\leq CJ^{-s}$ for some $C$ and any $J$.
Note that the validity of \autoref{condition5}(i) depends on the choice of the basis functions. \citet{chen2015optimal} show that such condition holds for splines and wavelets. \autoref{condition5}(ii) holds when $F_0$ satisfies some smoothness conditions, see \citet{chen2007large}. Finally, note that the uniform approximation rate in \autoref{condition5}(ii) sometimes applies only to the case where the support of $F_0$ is bounded\footnote{For example, this is required by splines, wavelet, or Chebyshev polynomials. But we point out that for the case of unbounded support, Hermite polynomials can be used.}, which amounts to requiring that $(x_0, X)$ lies in a compact set. To allow for unbounded support of $(x_0, X)$, we may transform $z_0$ into a bounded set. For example, $2\pi^{-1}\arctan(z_0)\in (-1, 1)$ regardless of the support of $z_0$. Under such transformation, to approximate $F_0(z)$, it remains to approximate $G_0(z) = F_0(\tan(\pi z/2))$. In this case uniform bounded derivative of $G_0$ may imply stronger condition on $F_0$.
\end{remark}
In this section, we stick with the previous setup where in the $k$-th round of update, we observe a batch of $B$ i.i.d. observations $ W_k = \{(x_{0,i,k}, X_{i,k}, Y_{i,k})\}_{i=1}^B$.
To motivate our online sieve estimator for function $F_0$, consider the most special case where $F_0(z) = \sum_{j=1}^{J^*}\beta_{j,0}\psi_j(z)$ for some fixed $J^*$ and $\beta_{1,0},\cdots, \beta_{J^*, 0}\in \mathbb{R}$. In this case, $F_0(z)$ is an element of the finite-dimensional sieve space $\mathcal{S}_{J^*}$. Define $z_{0,i,k} = x_{0,i,k}+X^{\prime}_{i,k}\theta_0$ for $i= 1, 2, \cdots, B$, we have that
$
\mathbb E(Y_{i,k}|z_{0,i,k}) = \sum_{j=1}^{J^*}\beta_{j,0}\psi_j(z_{0,i,k})$ for each $i$.
This immediately leads to an OLS-type estimator for $\mathcal{B}_{J^*,0} = (\beta_{1,0}, \cdots, \beta_{J^*, 0})^{\prime}$, which can be formulated as an online estimator as follows
\[
\hat{\mathcal{B}}_k = \hat{\mathcal{B}}_{k-1} + \frac{\eta_k}{B}\sum_{i=1}^B \left(Y_{i,k} - \Psi_{J^*}(z_{0,i,k})^{\prime}\hat{\mathcal{B}}_{k-1}\right)\Psi_{J^*}(z_{0,i,k}),
\]
where $\eta_k>0$ is positive learning rate and $\hat{\mathcal{B}}_0$ is given. The above update obviously fits in the parametric online learning framework, and hence the properties of the estimator can be easily derived.
In the most general case, there may not be such a truncation integer $J^*$ that $F_0(z)$ can be exactly represented by a linear combination of the first $J^*$ sieve functions. Moreover, the parametric component $\theta_0$ in the model is also unknown, so $z_{0,i,k}$ indeed is impossible to construct. We propose the following solutions. For the first issue, we consider a sequence of non-decreasing integers that are indexed by $k$, $\{J_k\}_{k=1}^{\infty}$, and in the $k$-th round of update we use a total of $J_k$ basis functions for approximation. Since $J_k$ is non-decreasing with respect to $k$, we increase the effective number of basis functions as the updates proceed, which is in spirit equivalent to increasing the number of basis functions when more data points are available in the offline full sample sieve estimation. For the second issue, in the $k$-th round of update we can use an estimator of $\theta_0$ as the replacement of unknown parameter to perform the update. The estimator, denoted as $\check{\theta}_{k-1}$, can be constructed with flexibility as long as it converges fast enough; see conditions below. For example, we can use the PR average $\bar{\theta}_{k-1}$ proposed in the previous section, or the rate-optimal estimator proposed in the next section.
Different from the fixed-dimensional update, the dimension of the sieve coefficients increases as the update proceeds. When $k$ is such that $J_k>J_{k-1}$, the dimension of sieve space increases across consecutive updates, so the dimension of $\hat{\mathcal{B}}_k$ is strictly larger than that of $\hat{\mathcal{B}}_{k-1}$. In this case, we need to specify the starting point of the coefficient of the new sieve function $\psi_{J_k}(\cdot)$. We propose to choose the starting point of the new sieve coefficient as 0 for each of the new sieve function.
Define $\check{z}_{i,k} = x_{0,i,k} + X_{i,k}^{\prime}\check{\theta}_{k-1}$. Our update of the sieve coefficient is formally described as follows
\begin{align}\label{online_sieve}
\hat{\mathcal{B}}_{k} = \begin{cases}
\hat{\mathcal{B}}_{ k-1} + \frac{\eta_k}{B}\sum_{i=1}^B \left( Y_{i,k} - \Psi_{J_{k-1}}\left(\check{z}_{i,k}\right)^{\prime}\hat{\mathcal{B}}_{k-1} \right)\Psi_{J_{k-1}}\left(\check{z}_{i,k}\right), \ \ \ \ \ \ \ \ \text{if $J_{k} = J_{k-1}$} \\
(\hat{\mathcal{B}}_{k-1}^{\prime}, 0)^{\prime} + \frac{\eta_k}{B}\sum_{i=1}^B \left( Y_{i,k} - \Psi_{J_{k}}\left(\check{z}_{i,k}\right)^{\prime}(\hat{\mathcal{B}}_{k-1}^{\prime}, 0)^{\prime} \right)\Psi_{J_{k}}\left(\check{z}_{i,k}\right), \ \text{if $J_{k} = J_{k-1} + 1$}
\end{cases}
\end{align}
and the online estimator of $F_0(\cdot)$ is given by
\begin{equation}\label{online estimator of F}
\hat{F}_k (\cdot) = \Psi_{J_k}^{\prime}(\cdot)\hat{\mathcal{B}}_k.
\end{equation}
Based on (\ref{online estimator of F}), we can also define the Polyak-Ruppert average estimator for $F_0$. A direct construction is to average across all updates of $F_0$, which is simply given by
\begin{equation}\label{pr average of F_0 naive}
\bar{F}_N(\cdot) = \frac{1}{N}\sum_{k=1}^N \hat F_k(\cdot).
\end{equation}
While the above expression is difficult to implement in practice. Note that based on (\ref{pr average of F_0 naive}), to evaluate the estimated function at a specific value, we need to evaluate the estimators across past updates. This makes it inconvenient to evaluate $\bar{F}_N$ at many points or at some new points.
\paragraph{Online Construction of the Sieve.} For implementation purpose, we propose to formulate the PR average of $F_0$ by averaging across the sieve coefficients $\hat{\mathcal{B}}_k$. Note that this is different from the fixed-dimensional parameter update considered in the previous section. One of the key features of the online sieve estimation is the increasing dimensions of sieve functions as the update proceeds. So the sieve coefficient estimators obtained in different updates may have different dimensions, which can not be simply averaged across time. To properly define the PR average of the sieve coefficient, we introduce some additional notations. For any $k$ and any $ m\times 1$ vector $a$ such that $l\leq J_k$, define $[a]_k$ to be a $J_k\times 1$ vector such that \begin{equation} [a]_k = (a^{\prime}, 0, \cdots, 0)^{\prime}.\end{equation} Note that for any $m\times 1$ vectors $a$ and $b$ and arbitrary constants $c_1, c_2$, we have that $[c_1\cdot a + c_2\cdot b]_k = c_1\cdot [a]_k + c_2\cdot [b]_k$. Then the PR average of sieve coefficient is iteratively defined by $\bar{\mathcal{B}}_{0} = \hat{\mathcal{B}}_{0}$, and
\begin{align}
\bar{\mathcal{B}}_{k} = \frac{1}{k}\hat{\mathcal{B}}_k + \frac{k-1}{k}[\bar{\mathcal{B}}_{k-1}]_k.
\end{align}
For any $N$, we have
\begin{align}
\bar{\mathcal{B}}_N = \frac{1}{N}\sum_{k=1}^N [\hat{\mathcal{B}}_k]_N.
\end{align}
It's straightforward to verify that
\begin{align}\label{pr average of F_0}
\bar{F}_N (\cdot) = \Psi_{J_N}^{\prime}(\cdot)\bar{\mathcal{B}}_N
\end{align}
We point out that formulation (\ref{pr average of F_0}) will be used to study the statistical properties of the PR average estimator $\bar{F}_N(\cdot)$.
\subsection{Statistical Properties}
We impose the following conditions regarding the choices of the learning rate, sieve functions, and sieve orders.
\begin{condition}\label{condition6} (i) $\Vert \psi_j\Vert_{\infty}\leq C$ and $\Vert \partial \psi_j\Vert_{\infty}\leq Cj^{\alpha_{\psi}}$, where $\alpha_{\psi}$ is a positive constant; (ii) $0<\inf_{J}\underline{\lambda}(\Gamma_J)\leq \sup_{J}\bar{\lambda}(\Gamma_J)<\infty$.
\end{condition}
\begin{condition}\label{condition7}
$\eta_k$ and $J_k$ are chosen such that $\eta_k = \eta_0k^{-\alpha_{\eta}}$ and $J_k =[J_0k^{\alpha_{J}}]$, where $\eta_0, J_0, \alpha_{\eta}$ and $\alpha_{J}$ are all positive constants, and $J_0$ is a positive integer.
\end{condition}
\begin{condition}\label{new_theta_condition_for_sieve}
$\check{\theta}_k = O(k^{-\alpha_{\theta}}) \ a.s.$ up to some poly-log terms for some $\alpha_{\theta} >0$. Moreover, for each $k$, $\check{\theta}_k$ is $\sigma(W_1,\cdots, W_k)$-measurable.
\end{condition}
Define set
\begin{align*}
\mathcal{A}_{\alpha} & = \left\{\left(\alpha_{\eta}, \alpha_J, \alpha_{\psi}, s\right): 1>\alpha_{\eta}>2\alpha_J, \alpha_{\theta}>\alpha_J(1+ \alpha_{\psi}), s>1, \right.\\
&\left. \ \ \ \ \alpha_{\eta}+ (2s - 1)\alpha_{J}>1, \alpha_{\eta} + 2\alpha_{\theta} - \alpha_J(1+2\alpha_{\psi})> 1, 2\alpha_{\eta} > 1+ \alpha_J\right\}.
\end{align*}
For any $k\geq 1$, define $\Delta\hat{\mathcal{B}}_{k} = \hat{\mathcal{B}}_{k} - \mathcal{B}_{J_k, 0}$,
where $\mathcal{B}_{J_k, 0}$ was defined in (\ref{L2-projection}). We have the following result.
\begin{theorem}\label{theorem4}
Let \autoref{condition5}--\autoref{new_theta_condition_for_sieve} hold. If $(\alpha_{\theta}, \alpha_{\eta}, \alpha_J, \alpha_{\psi}, s)\in\mathcal{A}_{\alpha}$, $F_0$ has bounded derivative, and $\mathbb{E}\Vert X \Vert^2<\infty $, we have that \[\Vert \Delta\hat{\mathcal{B}}_{N}\Vert^2\longrightarrow_{a.s.} 0.\] Moreover, define $\underline{\alpha}_F = \min\{\alpha_{\eta}+ (2s - 1)\alpha_{J}-1, \alpha_{\eta} + 2\alpha_{\theta} - \alpha_J(1+2\alpha_{\psi}) - 1, 2\alpha_{\eta} - \alpha_J -1,(s -1)\alpha_J\}$, we have that for any $r>1$, there holds \[\Vert \Delta\hat{\mathcal{B}}_{N}\Vert = O(N^{ -\frac{\underline{\alpha}_F}{2}}\log^{\frac{r}{2}}(N)), \ a.s. \]
\end{theorem}
The rate in \autoref{theorem4} is rough. In the following, we will first refine the convergence rate of $\Vert \Delta\hat{\mathcal{B}}_N\Vert $. Then we will show that, similar to the finite-dimensional case in Section \ref{section2}, PR average will lead to optimal rate of our nonparametric estimator. We first introduce another technical condition.
\begin{condition} \label{condition8}
$\mathbb E|\varepsilon|^{\kappa} <\infty$ with $\kappa > 2 + 2/\alpha_{\eta}$.
\end{condition}
\begin{remark}
\autoref{condition8} further regulates the tail behavior of the error term $\varepsilon$. When $\alpha_{\eta}$ is chosen close to 1, the required order of the moments will be slightly higher than 4. Note that such condition is used to build a.s. consistency and supreme-norm rate for the sieve estimator. If we pursue probability rates, such condition can be readily weakened to $\kappa > 2$ as in \citet{chen2015optimal}.
\end{remark}
Based on the above conditions, we can provide a more refined convergence rate for $\Delta\hat{\mathcal{B}}_N$.
\begin{theorem}\label{theorem5}
Let all requirements in \autoref{theorem4} hold. Then if further $(\frac{1}{2}+\alpha_J)\vee (1 +2\alpha_J(1+\alpha_{\psi})-2\alpha_{\theta})<\alpha_{\eta}<(2\alpha_{\theta} - 2\alpha_J\alpha_{\psi})\wedge (\alpha_J(2s+1))$, we have that
\[
\Vert \Delta\hat{\mathcal{B}}_N \Vert = O\left( N^{-\frac{\alpha_{\eta} -\alpha_J}{2}}\log^{\frac{1}{2}}(N)\right), \ a.s.
\]
\end{theorem}
For each $N$, define \begin{align}
\Delta\bar{\mathcal{B}}_N = \bar{\mathcal{B}}_N - \mathcal{B}_{J_N, 0}.
\end{align}
We have the following result.
\begin{theorem}\label{theorem6}
Let all requirements in \autoref{theorem5} hold. Suppose further that $4\alpha_J < \alpha_{\eta} < 1- 4\alpha_J$, $\alpha_{\theta}>\max\{(1-s_1)\alpha_{J}+\frac{1}{2}, (1-s+\alpha_{\psi})\alpha_{J}+\frac{1}{2}, (\frac{3}{4} + \frac{\alpha_{\psi}}{2})\alpha_J+ \frac{1}{4}\}$, $F_0$ has bounded second-order derivative, $\mathbb{E}\Vert X\Vert^4<\infty$, and $\overline{\lim}_k\sup_{z\in\mathbb{R}}\Vert \mathbb{P}_{J_k}(\partial F_0(\cdot)\mu_{0}(\theta_0, \cdot )^{\prime})(z) - \partial F_0(z)\mu_{0}(\theta_0, z)^{\prime}\Vert_{\infty} =0$ with $\mu_0(\theta_0, z) = \mathbb{E}(X|z_0 = z)$ and $\sup_{z\in\mathbb{R}}\Vert \partial F_0(z)\mu_0(\theta_0, z)\Vert<\infty$, we have that
\begin{align*}
\Delta\bar{\mathcal{B}}_N & = \frac{1}{N}\sum_{k=1}^N\frac{1}{B}\sum_{i=1}^B \varepsilon_{i,k} \left[\Gamma_{J_{k}}^{-1}\Psi_{J_k}(z_{0,i,k})\right]_N - \frac{1}{N}\sum_{k=1}^N\mathbb E\left(\partial F_0(z_{0}) \left[\Gamma_{J_{k}}^{-1}\Psi_{J_k}(z_{0})\right]_N X^{\prime}\right)\Delta\check\theta_{k-1} \\
& \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ + \frac{1}{N}\sum_{k=1}^N \left([\mathcal{B}_{J_{k}, 0}]_N - \mathcal{B}_{J_N, 0}\right) + o\left(N^{-\frac{1+\alpha_J}{2}} \right), \ a.s.
\end{align*}
and
\[
\left\Vert \bar{F}_N(z) - F_0(z)\right\Vert_{\infty} = O\left(N^{-\frac{1 - \alpha_J}{2}}\log^{\frac{1}{2}}(N)+ N^{-\alpha_{J}s} + N^{- \alpha_{\theta}} \right), \ a.s.
\]
where the last term of the supreme-norm rate is up to some poly-log terms determined by the convergence rate of $\check\theta_k$.
\end{theorem}
\autoref{theorem6} is the key result of online sieve estimation with generated regressor. It first gives the first-order expansion of the average online estimator $\bar{\mathcal{B}}_N$, which contains four terms. The first term is comparable to variance term of the full-sample sieve estimator (e.g., \citet{chen2015optimal}), but since in the $k$-th round of update we only use $J_k$ sieve functions, the effective variance component in the $k$-th round is $\varepsilon_{i,k}\Gamma_{J_k}^{-1}\Psi_{J_k}(z_{0,i,k})$ instead of $\varepsilon_{i,k}\Gamma_{J_{N}}^{-1}\Psi_{J_N}(z_{0,i,k})$. The second term in the expansion describes the impacts of using generated regressor $\check z_{i,k}$ for sieve estimation. The third term describes the impacts on changing pseudo true sieve coefficients when the dimension of sieve increases.
When we consider the supreme-norm estimation error of $F_0$, we see that apart from the error caused by generated regressor $\check z_{i,k}$, the error rate is optimal in the sense that if we choose $\alpha_J = 1/(2s+1)$, then
\[
N^{-\frac{1 - \alpha_J}{2}}\log^{\frac{1}{2}}(N)+ N^{-\alpha_{J}s} \asymp N^{-\frac{s}{2s+1}}\log^{\frac{1}{2}}(N),
\]
and such error rate is optimal up to some poly-log terms \citep*{stone1982optimal,belloni2015some,chen2015optimal}. Note that when the plug-in estimator of $\theta$ converges fast enough, that is, $\alpha_{\theta} > \frac{s}{2s+1}$, we have that the plug-in error does not affect the supreme-norm rate.
\section{Learning Phase II: Rate-Optimal Learning}\label{section4}
This learning phase simultaneously updates estimators of $\theta_0$ and $F_0$, and produces estimators with optimal convergence rate. In Section \ref{section4.1}, we develop a Neyman-orthogonalized update algorithm, which takes the nonparametric estimator of $F_0$ as a plug-in and outputs trajectories of estimators of $\theta_0$. When the nonparametric estimation error meets some rate requirements, the orthogonalized update produces $1/\sqrt{N}$-consistent estimators for $\theta_0$. In Section \ref{Section 4.2}, we combine the orthogonalized update with the learning of $F_0$ based on online method of sieves proposed in Section \ref{section3}. This completes the second learning phase. We show that such combination delivers estimators of $\theta_0$ and $F_0$ with optimal convergence rate. Finally, in Section \ref{Section 4.3}, we discuss using trajectories to conduct online inference based on random scaling method.
\subsection{The Orthogonalized Update}\label{section4.1}
In Section \ref{section2}, we have shown that $\bar{\theta}_N$ is $1/\sqrt{N^{1 - \alpha_h }}$-consistent up to some poly-log terms. The convergence rate is obviously slower than the optimal rate $1/\sqrt{N}$. In this section we construct an algorithm that produces estimators that converge at $1/\sqrt{N}$ rate. According to the discussion following \autoref{theorem3}, the main cause of the slower convergence rate of $\Delta\hat\theta_N$ and $\Delta\bar{\theta}_N$ is the diverging variance of the score function $\Phi_k(\theta, W_k)$ used to perform the update. To improve the convergence rate, we resort to using a score function with bounded variance. In this section we consider the following score function
\begin{align}\label{score}
\varphi(\theta, x_0, X, Y, F_0) = \left( Y - F_0\left(x_0 + X^{\prime}\theta\right)\right)X,
\end{align}
and the associated update
\begin{align}\label{original update}
\tilde{\theta}_k = \tilde{\theta}_{k-1} + \frac{\xi_k}{B}\sum_{i=1}^B \left( Y_{i,k} - F_0\left(x_{0,i,k}+ X_{i,k}^{\prime}\tilde{\theta}_{k-1}\right)\right)X_{i,k},
\end{align}
where $\tilde\theta_0$ is the initial point and is given, and $\xi_k>0$ is the learning rate. The score function (\ref{score}) and corresponding update (\ref{original update}) are commonly used for monotone index models \citep*{agarwal2014least,khanetal2021,khan2024inference}. However, $F_0$ has to be estimated in the semiparametric setup. Note that we can simply use the online sieve estimator $\bar{F}_k$ proposed in Section \ref{section3}, but such estimator converges at a rate slower than $1/\sqrt{k}$ according to \autoref{theorem6}. As a result, the estimation error of the sieve estimator will remain and contaminate the update (\ref{original update}), which prevents us from getting $1/\sqrt{N}$-consistent estimator for $\theta_0$. This issue is especially serious when the plug-in estimator $\check{\theta}_k$ is not $1/\sqrt{N}$-consistent.
To deal with the above issue, we need to remove the first-order impacts of sieve estimation on the update. A straightforward way to so is to perform Neyman-orthogonalization for score function (\ref{score}), and use the orthogonalized score function for update. In particular, define $\mu_0(\theta, z) = \mathbb{E}( X| x_0 + X^{\prime}\theta = z)$, we consider the Neyman-orthogonalized score function given by
\begin{align}\label{Neyman score}
\tilde\varphi(\theta, x_0, X, Y, F_0, \mu_0) = \left(Y - F_0\left(x_0 + X^{\prime}\theta\right)\right)\left(X - \mu_{0}\left(\theta, x_0 + X^{\prime}\theta\right)\right).
\end{align}
Using (\ref{Neyman score}) will effectively remove the impact of sieve estimation of $F_0$ on the update of $\theta$. However, when we use such score function, two additional issues arise. The first issue is computational-wise. Note that the conditional expectation function $\mu_{0}(\theta, z)$ is unknown so has to be estimated. Such function has two inputs $z$ and $\theta$. Using the online sieve estimation proposed in Section \ref{section3}, we can easily provide an estimate for such function at any fixed parameter $\theta$. However, estimating the function that also takes $\theta$ as a varying input is computationally intensive, especially when the number of regressors is large (where tensors of sieve functions are required). The second issue is algorithm-wise. Note that for any fixed $\theta$, we have
\begin{align*}
\tilde\varphi(\theta, x_0, X, Y, F_0, \mu_0) & = \varepsilon\left(X - \mu_{0}\left(x_0 + X^{\prime}\theta\right)\right) \\
& - \int_{0}^1\partial F_0(z_0 +\tau X^{\prime}\Delta\theta)d\tau \left(X - \mu_{0}\left(\theta, x_0 + X^{\prime}\theta\right)\right)X^{\prime}\Delta\theta.
\end{align*}
where recall that $z_0 = x_0 + X^{\prime}\theta_0$. The first term on RHS is mean zero so only contributes to the variance of the update. The second term can be written as $M(\theta, x_0, X)\Delta\theta$ with \[M(\theta, x_0, X) = \int_{0}^1\partial F_0(z_0 +\tau X^{\prime}\Delta\theta)d\tau \left(X - \mu_{0}\left(\theta, x_0 + X^{\prime}\theta\right)\right)X^{\prime}.\] Note that the matrix $\mathbb EM(\theta,x_0, X)$ is not guaranteed to be positive-(semi)definite for arbitrary parameter $\theta$; indeed, for generally distributed $X$ and arbitrary $\theta$, $\mathbb EM(\theta, x_0, X)$ is even not guaranteed to be symmetric. This imposes challenges on the validity of the updating algorithm due to the lack of guaranteed contraction mapping.
To deal with the above two issues, we point out a key observation in the following equation
\begin{align}\label{key observation}
\mathbb EM(\theta_0, x_0, X) & = \mathbb E \left(\partial F_0(z_0)(X-\mu_{0}(\theta_0, z_0))X^{\prime}\right) \nonumber \\
& = \mathbb E \left(\partial F_0(z_0)(X-\mu_{0}(\theta_0, z_0))(X-\mu_{0}(\theta_0, z_0))^{\prime}\right)\equiv \mathcal{M}_0.
\end{align}
Under mild conditions, $\mathcal{M}_0$ is positive definite. This implies that the update is guaranteed to (on average) shrink $\theta$ towards $\theta_0$ if $\theta$ is close to $\theta_0$. Moreover, recall that the computational issue prevents us from estimating $\mu_0(\theta, z)$ for arbitrary $\theta$. As a result, if we can effectively ``gauge'' our update closely around $\theta_0$, we only need to estimate $\mu_{0}(\theta, z)$ at the true parameter $\theta_0$. Then the unknown function $\mu_{0}(\theta_0, z)$ becomes a univariate function of $z$, which is easy to estimate.
Motivated by the above discussion, we propose to construct a sequence of carefully designed ``gauge balls'' denoted by $\{\varTheta_k\}_{k=1}^{\infty}$, and use these neighborhoods to gauge our update so that the updated parameter will be closely centered around $\theta_0$. In particular, such sequence of gauge balls have shrinking diameters and they contain the true parameter $\theta_0$ almost surely. These properties are formally stated in the following condition.
\begin{condition}\label{condition9} $\varTheta_k\subseteq \mathbb{R}^p$ is convex, compact, and uniformly bounded, and there exists $d_k\downarrow 0$ such that
\[\lim\sup_{k\rightarrow\infty}d^{-1}_k\mathrm{diam}(\varTheta_k) \leq 1, \ a.s.\] and \[\lim\inf_{k\rightarrow \infty} \boldsymbol{1}\left(B(\theta_0, c_{\varTheta}d_k)\subseteq \varTheta_k\right) = 1, \ a.s. \] for some $0<c_{\varTheta}<1$. Finally, $\varTheta_k$ is $\sigma(W_1,\cdots, W_{k-1})$-measurable.
\end{condition}
The gauge ball can be simply constructed using the estimator $\bar{\theta}_N$ in the first phase. Given such sequence of gauge balls, we now describe how we gauge our update to make it centered around the true parameter $\theta_0$. For any convex and compact set $\varTheta\subseteq \mathbb{R}^p$ and $\theta\in \mathbb{R}^p$, define projection on to set $\varTheta$
\begin{align}
\Pi_{\varTheta} \theta = \arg\min_{\vartheta\in\varTheta} \left\Vert \vartheta - \theta \right\Vert.
\end{align}
For any positive integer $k$, define $\Pi_k(\cdot) = \Pi_{\varTheta_k}(\cdot)$ as the projection to set $\varTheta_k$. We modify the update (\ref{original update}) as follows
\begin{align}\label{gauged update}
\tilde{\theta}_k = \Pi_{k}\left[\tilde{\theta}_{k-1} + \frac{\xi_k}{B} \sum_{i=1}^B \left(Y_{i,k} - \check{F}_{k-1}\left(x_{0,i,k} + X_{i,k}^{\prime}\tilde{\theta}_{k-1}\right) \right)\left(X_{i,k} - \check{\mu}_{k-1}\left(x_{0,i,k} + X_{i,k}^{\prime}\check\theta_{k-1}\right)\right)\right].
\end{align}
where $\check{F}_{k-1}(\cdot)$ and $\check{\mu}_{k-1}\left(\cdot\right)$ are estimators of $F_0(\cdot)$ and $\mu_{0}(\theta_0, \cdot)$, and $\check{\theta}_{k-1}$ is the estimator of $\theta_0$. We note that $\check{F}_{k-1}(\cdot)$, $\check{\mu}_{k-1}\left(\cdot\right)$, and $\check{\theta}_{k-1}$ can be constructed flexibly.
To study the asymptotic behavior of the proposed estimator $\tilde\theta_k$, we impose the following conditions.
\begin{condition}\label{condition10} (i) $0<\lambda(\mathcal{M}_0)<\infty$; (ii) $\mathbb{E}\Vert X\Vert^6<\infty$, $\mathbb{E}\Vert X\Vert^4\Vert \mu_{0}(\theta_0, z_0)\Vert^2<\infty $, and $\sup_{x_0, X}\mathbb{E}(\varepsilon^2|x_0, X)<\infty$; (iii) $\mu_{0}(\cdot, \theta_0)$ has uniformly bounded first derivative, $F_0$ has uniformly bounded first and second derivative; (iv) $\check{F}_{k-1}(\cdot)$, $\check{\mu}_{k-1}(\cdot)$, $\check{\theta}_{k-1}$ are $\sigma(W_1, \cdots, W_{k-1})$-measurable. There exist positive constants $\alpha_{\xi}, \alpha_{d}, \alpha_{F}, \alpha_{\partial F},$ $\alpha_{\mu}$, and $\alpha_{\theta}$ such that $
\xi_k = \xi_0k^{-\alpha_{\xi}}$, $
d_k = d_0 k^{-\alpha_d}$, and there a.s. hold $
\Vert \check{F}_k - F_0 \Vert_{\infty} = O\left(k^{-\alpha_{ F}}\right)$, $
\Vert \partial \check{F}_k - \partial F_0 \Vert_{\infty} = O\left(k^{-\alpha_{\partial F}}\right)$, $
\Vert \check{\mu}_k - \mu_{0}(\theta_0,\cdot) \Vert = O\left(k^{-\alpha_{\mu}}\right)$, and $\Vert \Delta\check{\theta}_k\Vert = O(k^{-\alpha_{\theta}})$ up to some poly-log terms.
Moreover, $\frac{1}{2}<\alpha_{\xi}<1$, $2\alpha_{\partial F} +\alpha_{\xi}>1$, $2\alpha_{\mu}\wedge \alpha_{\theta} + \alpha_{\xi}>1$, and $2\alpha_F + 2\alpha_{\mu}\wedge\alpha_{\theta} > 1$. Finally, define $\omega = \alpha_{\xi} + \min\{\alpha_{\xi}, 2\alpha_{F}+\min\{2\alpha_{\mu}, 2\alpha_{\theta}\}\}$,
there hold
$\omega > 1 + 2\alpha_d$.
\end{condition}
Based on the above conditions, we now provide the statistical properties of of $\tilde{\theta}_k$. Define
\begin{equation}
\bar{\tilde\theta}_N = \frac{1}{N} \sum_{k=1}^N \tilde\theta_k,
\end{equation}
and
\begin{equation} \mathcal{V}_1 = \mathbb{E}\left[\varepsilon^2 \left(X - \mu_{0}(\theta_0, z_0)\right) \left(X - \mu_{0}(\theta_0, z_0)\right)^{\prime}\right],\end{equation} we have the following theorem.
\begin{theorem}\label{theorem7} If \autoref{condition9} and \autoref{condition10} hold, we have that
\begin{equation}
\bar{\tilde\theta}_N - \mathcal{M}_0^{-1}\frac{1}{N}\sum_{k=1}^N \frac{1}{B}\sum_{i=1}^B\left( X_{i,k} - \mu_{0}(\theta_0, z_{0,i,k})\right)\varepsilon_{i,k} = o\left(N^{-\frac{1}{2}}\right), \ a.s.
\end{equation}
Consequently, we have that
\begin{equation}
\underline{\overline{\lim}}_{N\rightarrow \infty}\left(B^{-1}\mathcal{M}_0^{-1}\mathcal{V}_1\mathcal{M}_0^{-1}\right)^{-\frac{1}{2}} \frac{\sqrt N\bar{\tilde\theta}_N}{\sqrt{ \log(\log(N))}} = \pm \boldsymbol{1}_{p}, \ a.s.,
\end{equation}
\begin{equation}
\left(B^{-1}\mathcal{M}_0^{-1}\mathcal{V}_1\mathcal{M}_0^{-1}\right)^{-\frac{1}{2}}\sqrt N\bar{\tilde\theta}_N \longrightarrow_d N(0, \mathbb{I}_p),
\end{equation}
and
\begin{equation}\label{FCLT for N-O theta}
\left(B^{-1}\mathcal{M}_0^{-1}\mathcal{V}_1\mathcal{M}_0^{-1}\right)^{-\frac{1}{2}}\sqrt {[Nr]}\bar{\tilde\theta}_{[Nr]} \Longrightarrow \mathbb{W}_p(r), \ r\in [0, 1].
\end{equation}
\end{theorem}
\subsection{ Rate-Optimal Learning}\label{Section 4.2}
Given the Neyman-orthogonalized update proposed in the previous section, this section describes the rate-optimal learning phase in detail and its properties. Suppose that in period $k-1$, we are given $\tilde\theta_{k-1}$, $\bar{\tilde \theta}_{k-1}$, and the PR averages $\hat F_{k-1}$, and $\bar F_{k-1}$. Then we perform the following updates.
\begin{enumerate}
\item Choose $\check{\theta}_{k-1} = \tilde\theta_{k-1}$, and $\check{F}_{k-1} = \bar F_{k-1} $, and update $\tilde\theta_{k-1}$ to $\tilde\theta_k$ by (\ref{gauged update});
\item Update $\hat F_{k-1}$ to $\hat F_{k}$ by (\ref{online_sieve});
\item Update the PR averages to $\bar{\tilde \theta}_{k}$ and $\bar F_k$.
\end{enumerate}
We have the following theorem for the second learning phase.
\begin{theorem}\label{rate-optimal}
Let \autoref{condition5}--\autoref{condition9}, \autoref{condition10}(i)--(iii), and all the requirements in \autoref{theorem6} hold. If \autoref{condition10}(iv) holds with $\alpha_{F} = \min\{\frac{1-\alpha_J}{2}, s\alpha_J, \alpha_{d}\}$, $\alpha_{\partial F} = \min\{\frac{1-\alpha_J(1 + 2\alpha_{\psi})}{2}, s_1\alpha_J, \alpha_{d} - \alpha_J\alpha_{\psi}\}$, and $\alpha_{\theta} = \alpha_d$, then we have that \autoref{theorem7} holds, and the expansion of $\Delta\bar{\mathcal{B}}_N$ in \autoref{theorem6} holds. Consequently, the rate of $\Vert \bar F_N - F_0\Vert_{\infty}$ can be improved to $O(N^{-\frac{1-\alpha_J}{2}}\log^{\frac{1}{2}}(N) + N^{-\alpha_Js})$. So both updates are rate-optimal.
\end{theorem}
\subsection{Online Inference based on Random Scaling}\label{Section 4.3}
The inference in semiparametric estimation with plug-in estimated variance-covariance matrix estimator is often computationally exhaustive due to intensive nonparametric estimation uisng full sample. This issue becomes important when the sample size is too large for storage. Fortunately, as a consequence of our functional central limit theorem (FCLT) in \autoref{theorem7} (\ref{FCLT for N-O theta}), we can apply a computationally super efficient algorithm (see Algorithm \ref{alg:rs_inference}) to conduct online inference for $\theta_0$. This algorithm is based on the random scaling method proposed by \citet*{lee2022fast}. The random scaling only takes the updated PR averages $\bar{\tilde{\theta}}_1, \bar{\tilde{\theta}}_2, \cdots$ as input, and automatically produces the confidence intervals for $\theta_0$, so the inference requires almost no additional computational time. \citet {lee2022fast} initially proposed this random scaling method for their online parametric least squares estimators, and is subsequently used in \citet*{chen2023sgmm,chen2025slim} for online generalized method of moment estimators. Notice that the theoretical justification of the random scaling just relies on the validity of a FCLT, which we establish for our online estimator in \autoref{theorem7}. We summarize the inference for our online semiparametric estimator via random scaling in Algorithm \ref{alg:rs_inference} in Section \ref{implementation} for easy reference.
\section{Practical Implementations}\label{implementation}
This section provides detailed discussion on the practical implementations of the algorithms proposed in the previous sections.
The warm-start phase estimation proceeds with kernel-based online learning proposed in Section \ref{section2}. Due to the global contraction mapping properties indicated in \autoref{lemma2} and \autoref{theorem1}, we can use any initial guess as the starting point of such learning phase. An option can be the origin point $\boldsymbol{0}_p$. When conducting first learning phase, it's recommended to use Epanechnikov kernel function. The order of the kernel function should be chosen depending on the order of smoothness of the unknown function. In simulations and empirical applications, we use 6-th order Epanechnikov kernel function. Given the choice of the kernel function, it's recommended to choose a large bandwidth to improve finite sample performance. For example, when we choose $s$-th order kernel, the bandwidth can be chosen as $h_k = c_k\cdot k^{-\frac{1}{s-1}}$, where $c_k$ is some positive constant depending on $k$. A rule-of-thumb choice of $c_k$ is the estimator of standard deviation of $z_{i,k}(\hat\theta_{k-1}) - z_{j,k}(\widehat\theta_{k-1})$ with $i\neq j$. In particular, in the $k$-th update, define $\hat z_{i,k} = x_{0,i}+ X_{i}^{\prime}\hat\theta_{k-1}$ and $\hat\sigma_{z, k}$ as the sample standard deviation of $\hat z_{i, k}$, then $c_k$ should be chosen as \begin{equation}\label{rule-of-thumb-bandwidth}
c_k = 1.414 \hat\sigma_{z, k }.
\end{equation}
In practice, the warm-start phase learning will be performed only for a small number of rounds to quickly locate the neighborhood of the true parameter. For this purpose, the learning rate $\gamma_k$ can be chosen as a small constant such as 0.05 or 0.01 to speed up the convergence. We recommend starting the PR average update of $\theta$ when $\hat\theta_k$ starts fluctuating around some point. While we update the PR average of $\theta$, we can also produce warm-start estimators for $F_0(\cdot)$ and $\mu_{0}(\theta_0, \cdot)$, which will be used as the warm-start of the second learning phase. Let $N_0$ be the round from which we start conducting PR average update for $\theta$, and $N_1$ be the total number of updates in warm-start learning phase. Let $J_0$ be the initial number of sieve functions. Then the initial sieve coefficients of $F_0(\cdot)$ and $\mu_{0}(\theta_0, \cdot)$ are estimated by
\begin{equation}\label{Y_sieve}
\hat{\mathcal{B}}_{Y,N_1} = \left(\sum_{k=N_0}^{N_1}\sum_{i=1}^B\Psi(J_0, \hat z_{i,k})\Psi(J_0, \hat z_{i,k})^{\prime}\right)^{-1}\left(\sum_{k=N_0}^{N_1}\sum_{i=1}^B\Psi(J_0, \hat z_{i,k})Y_{i,k}\right),
\end{equation}
\begin{equation}\label{mu_sieve}
\hat{\mathcal{B}}_{X,N_1}= \left(\sum_{k=N_0}^{N_1}\sum_{i=1}^BX_{i,k}\Psi(J_0, \hat z_{i,k})^{\prime}\right)\left(\sum_{k=N_0}^{N_1}\sum_{i=1}^B\Psi(J_0, \hat z_{i,k})\Psi(J_0, \hat z_{i,k})^{\prime}\right)^{-1}.
\end{equation}
Note that (\ref{Y_sieve}) and (\ref{mu_sieve}) can be computed based on online update similar to the PR average update. The first learning phase is then summarized in Phase I in Algorithm \autoref{alg:no_online_update}.
\begin{algorithm}[t!]
\caption{Online Learning of $\theta_0$ and $F_0$}
\label{alg:no_online_update}
\begin{algorithmic}[1]
\Require
Sequential observations $W_1, \cdots, W_{N}$, number of initial updates $N_0$, number of warm-up updates $N_1$, kernel function $K(\cdot)$, bandwidth parameter $\alpha_h$, number of sieves $J_0$, sieve vector $\Psi(J, z)$, constant learning rate $\gamma_0$ for warm-start phase learning, and sequence of learning rates $\xi_k$ and $\eta_k$ for rate-optimal phase learning
\Ensure
Online estimators $\bar{\tilde\theta}_N$ and $\bar{F}_N(\cdot)$
\State \textbf{Learning Phase I: Warm-Start Phase }
\State Initialization: $\hat\theta_0 \leftarrow \boldsymbol{0}_p$, $\bar\theta_{N_0 - 1} \leftarrow \boldsymbol{0}_p$, $M_1 \leftarrow 0$, $M_2 \leftarrow 0$, $M_3 \leftarrow 0$;
\For{$k = 1, 2, \dots, N_1$}
\State $h_k \leftarrow c_k\cdot k^{-\alpha_h}$ with $c_k$ being define in (\ref{rule-of-thumb-bandwidth});
\State Update $\hat\theta_{k-1}$ to $\hat\theta_k$ using (\ref{first-step algorithm});
\If{$k\geq N_0$}
\State $\bar{\theta}_{k}\leftarrow \bar{\theta}_{k-1 }\cdot (1 - \frac{1}{k - N_0 + 1}) + \hat\theta_{k}\cdot \frac{1}{k - N_0 + 1}$;
\State $M_1 \leftarrow M_1 + \sum_{i=1}^B \Psi(J_0, \hat z_{i,k})\Psi(J_0, \hat z_{i,k})^{\prime}$, where $\hat z_{i,k} = x_{0,i,k} + X_{i,k}^{\prime}\hat\theta_k$;
\State $M_2 \leftarrow M_2 + \sum_{i=1}^B \Psi(J_0, \hat z_{i,k})y_{i,k}$;
\State $M_3 \leftarrow M_3 + \sum_{i=1}^B X_{i,k}\Psi(J_0, \hat z_{i,k}) ^{\prime}$;
\EndIf
\EndFor
\State \textbf{Learning Phase II: Rate-Optimal Phase}
\State Initialization: $\tilde\theta_{N_1} \leftarrow \bar\theta_{N_1}$, $\bar{\tilde\theta}_{N_1} \leftarrow \bar\theta_{N_1}$, $\hat{\mathcal{B}}_{Y,N_1} \leftarrow M_1^{-1}M_2$, $\hat{\mathcal{B}}_{X,N_1} \leftarrow M_3M_1^{-1}$;
\For{$k = N_1 + 1, \dots, N$}
\State Update $\tilde\theta_{k-1}$ to $\tilde\theta_k$ using (\ref{gauged update}) with gauge ball $\varTheta_k = \mathbb{R}^p$, and $\check{F}_{k-1}(\cdot)$ and $\check\mu_{k-1}(\cdot)$ being chosen in (\ref{nonparametric update});
\State Update sieve coefficients $\hat {\mathcal{B}}_{Y, k}$ and $\hat {\mathcal{B}}_{X, k} $ using (\ref{online estimator of F}) and (\ref{sieve_for_mu});
\State Update PR average estimators using (\ref{pr update});
\EndFor
\State \Return $\bar{\tilde\theta}_N$ and $\bar{F}_N(\cdot) = \Psi(J_N, \cdot)^{\prime}\bar{\mathcal{B}}_{Y,N}$
\end{algorithmic}
\vspace{3pt}
\makebox[\linewidth][l]{
\footnotesize
\parbox{\linewidth}{
\noindent
\textbf{Note:} The learning rate in rate-optimal learning phase should adjust to the number of warm-start phase learning. For example, $\xi_k$ should be chosen as $\xi_k = \xi_0\cdot (k - N_1)^{-\alpha_{\xi}}$ for some $\xi_0$ and $\alpha_{\xi}$.
}
}
\normalsize
\end{algorithm}
After $N_1$ rounds of warm-start phase learning, $\bar\theta_{N_1}$ will fall into a small neighborhood of the true parameter, and in this case, we can start the second phase of learning and perform updates for $\theta_0$ and $F_0$ simultaneously. Throughout discussion we will assume $J_k$, $\xi_k$ and $\eta_k$ are all given; we will discuss the choices of these tuning parameters in the end. For initialization, we set $\tilde\theta_{N_1} = \bar\theta_{N_1}$.
We also initialize the PR average estimators
$\bar{\tilde\theta}_{N_1} = \tilde\theta_{N_1}$, $\bar{\mathcal{B}}_{Y,N_1} = \hat{\mathcal{B}}_{Y,N_1}$, $ \bar{\mathcal{B}}_{X,N_1} = \hat{\mathcal{B}}_{X,N_1}
$.
Given all the initializations, we start the online learning procedure. In the $k$-th round with newly-arrived data $W_k$, we first update $\tilde\theta_{k-1}$ to $\tilde\theta_{k}$ using (\ref{gauged update}). Some details should be noted here. When we choose a large $N_1$, $\hat\theta_{N_1}$ will fall into a small neighborhood of $\theta_0$ almost surely, in which case the projection is no loner needed. So in practice, we can choose the the gauge ball $\varTheta_k$ to be $\mathbb{R}^p$ to get rid of projection. For the choice of $\check F_{k-1}(\cdot)$ and $\check \mu_{k-1}(\cdot)$ we use
\begin{equation}\label{nonparametric update}
\check{F}_{k-1} (\cdot) = \Psi(J_{k}, \cdot)^{\prime}[\bar{\mathcal{B}}_{Y,k-1}]_k, \ \ \check{\mu}_{k-1}(\cdot) = [\bar{\mathcal{B}}_{X,k-1}]_k^{\prime}\Psi(J_{k},\cdot),
\end{equation}
where $J_k$ is the number of sieve functions in the $k$-th update, and $[\bar{\mathcal{B}}_{X,k-1}]_k$ is understood to be adding a zero column the right side of matrix $\bar{\mathcal{B}}_{X,k-1}$ if the dimension of $\bar{\mathcal{B}}_{X,k-1}$ and $\Psi(J_k, \cdot)$ does not match. We note that the use of PR average estimators for the above two nonparametric estimators is to make the first-order error decay at the best rate.
After the update of $\tilde\theta_{k}$, we update $\hat{\mathcal{B}}_{Y,k-1}$ to $\hat{\mathcal{B}}_{Y,k}$ using (\ref{online estimator of F}), where $\check{\theta}_{k-1}$ is chosen as $\tilde\theta_{k-1}$. We also update $\hat{\mathcal{B}}_{X,k-1}$ to $\hat{\mathcal{B}}_{X,k}$ by
\begin{equation}\label{sieve_for_mu}
\hat {\mathcal{B}}_{X, k} = \hat {\mathcal{B}}_{X, k-1} + \frac{\eta_k}{B}\sum_{i=1}^B\left(X_{i, k} - [\hat{\mathcal{B}}_{X, k-1}]_k\Psi(J_k, \tilde z_{i,k} )\right)\Psi(J_k, \tilde z_{i,k} )^{\prime};
\end{equation}
where $\tilde z_{i,k} = x_{0,i,k}+X_{i,k}^{\prime}\tilde\theta_{k-1}$. We note that in the update of sieve coefficients, we use the unaveraged estimator $\tilde\theta_{k-1}$ instead of the PR average estimator. According to \autoref{theorem6}, such choice leads to the best convergence rate of PR average estimator of sieve coefficients.
After the above updates, we finally update the PR average estimator by
\begin{align}\label{pr update}
& \bar{\tilde{\theta}}_{k} = \bar{\tilde{\theta}}_{k-1 }\left(1 - \frac{1}{k - N_1}\right) + \frac{\tilde \theta_{k}}{k - N_1},\\
& \bar{\mathcal{B}}_{Y, k} = [\bar{\mathcal{B}}_{Y, k-1}]_k\left(1 - \frac{1}{k - N_1}\right) + \frac{\hat{\mathcal{B}}_{Y, k}}{k - N_1},\\
& \bar{\mathcal{B}}_{X, k} = [\bar{\mathcal{B}}_{X, k-1}]_k\left(1 - \frac{1}{k - N_1}\right) + \frac{\hat{\mathcal{B}}_{X, k}}{k - N_1}.
\end{align}
This concludes the $k$-th update.
Based on the PR average estimators of $\theta_0$ in the second phase, we can conduct online inference based on random scaling \citep{lee2022fast}. The algorithm is provided in Algorithm \autoref{alg:rs_inference}.
We finally discuss the choice of tuning parameters $J_k$, $\xi_k$ and $\eta_k$. In most of the empirical applications, we can choose a large $J_0$ such as 30, and then fix $J_k = J_0$ for all $k$. For $\xi_k$, we choose $\xi_k = \xi_0(k - N_1)^{-\alpha_{\xi}}$. We can choose $\xi_0 = 1$ and we recommend choosing a small $\xi$ to speed up the convergence. In simulation and applications we choose $\alpha_{\xi} = 0.51$. Finally, for $\eta_k$, we recommend choosing $\eta_k = \eta_0(k-N_1)^{\alpha_{\eta}}$ with $\eta_0 = 1$. When $J_k$ is fixed throughout updates, we can choose $\alpha_{\eta} = 0.61$.
\begin{algorithm}[t]
\caption{Online Inference of $\theta_0$ based on Random Scaling of \citet{lee2022fast}}
\label{alg:rs_inference}
\begin{algorithmic}[1]
\Require
Sequential updates of PR average $\bar{\tilde\theta}_{N_1+ 1}, \cdots, \bar{\tilde\theta}_N$
\Ensure
Sequential vectors of confidence bands $\nu_{N_1 + 1},\cdots, \nu_N$ (nominal coverage rate 0.95)
\State Initialization: $A\leftarrow 0$, $b\leftarrow 0$, $c\leftarrow 0$
\For{$k = N_1 + 1, \dots, N$}
\State Update $A, b, c$
\[
A\leftarrow A + (k - N_1)^2\bar{\tilde \theta}_k\bar{\tilde \theta}_k^{\prime},
\]
\[
b\leftarrow b + (k - N_1)^2\bar{\tilde \theta}_k,
\]
\[
c\leftarrow c + (k - N_1)^2,
\]
\State $V_k\leftarrow \left(A - \bar{\tilde \theta}_k\cdot b^{\prime} - b\cdot \bar{\tilde \theta}_k^{\prime} + c\cdot \bar{\tilde \theta}_k\bar{\tilde \theta}_k^{\prime} \right)/k^2$;
\State Construct vector of confidence band for $\bar{\tilde\theta}_k$
\[
\nu_k \leftarrow 6.747\cdot \sqrt{\text{diag}(V_k)/(k - N_1)}
\]
\EndFor
\State \Return $\nu_1,\cdots, \nu_N$
\end{algorithmic}
\vspace{3pt}
\makebox[\linewidth][l]{
\footnotesize
\parbox{\linewidth}{
\noindent
\textbf{Note:} For any matrix $V$, $\text{diag}(V)$ returns a vector consisting of the diagonal elements of $V$.
}
}
\normalsize
\end{algorithm}
\begin{comment}
\paragraph{Phase 1: Initialization (Warm Start)}
Before processing the data stream, we establish a reliable starting point using a small offline sample.
\begin{enumerate}
\item {Collect Warm-Up Sample:} Gather the first $N_{warm}$ observations.
\item {Initial Estimation ($\hat{\theta}_0$):} Compute the static estimator using the Stochastic Batched Gradient Descent (SBGD) algorithm (Khan et al., 2021).
\item \text{Calibrate Hyperparameters:}
\begin{itemize}
\item Compute $Z_i = X_i'\hat{\theta}_0$.
\item Calculate $c_z = \text{std}(Z_i - Z_j)$ from pairs in the warm-up sample.
\item Set initial bandwidth constant based on $c_z$.
\end{itemize}
\item {Initialize Tracks:} Set $\hat{\theta}_{robust, 0} = \hat{\theta}_0$ and $\tilde{\theta}_{efficient, 0} = \hat{\theta}_0$.
\end{enumerate}
\paragraph{Phase 2: The Online Loop}
For each new data batch $k = 1, 2, \dots$ arriving in the stream:
\begin{algorithm}[t]
\caption{Online Update Procedure at Step $k$}
\begin{algorithmic}[1]
\State \textbf{Input:} New data batch $W_k$, previous estimates $\hat{\theta}_{k-1}$, $\tilde{\theta}_{k-1}$.
\Statex \textit{// Step A: Robust Update (The Guardrail)}
\State Update bandwidth $h_k = c_z \cdot k^{-0.1}$ and rate $\gamma_k \propto k^{-0.66}$.
\State Compute Kernel Gradient $\hat{S}_k(\hat{\theta}_{k-1})$ using Eq. (2.3).
\State Update robust estimator: $\hat{\theta}_k \leftarrow \hat{\theta}_{k-1} + \gamma_k \hat{S}_k$.
\State Normalize: $\hat{\theta}_k \leftarrow \hat{\theta}_k / ||\hat{\theta}_k||$.
\Statex \textit{// Step B: Construct Safety Zone (Gauge Ball)}
\State Define feasible set $\Theta_k = \{ \theta : |\theta - \hat{\theta}_k| \le k^{-\alpha_\theta} \}$.
\Statex \textit{// Step C: Efficient Update (The Target)}
\State Compute Orthogonal Gradient update (Newton-step) for $\tilde{\theta}_{k}$.
\State \textbf{Projection:} If $\tilde{\theta}_{k} \notin \Theta_k$, project it back onto the boundary of $\Theta_k$.
\Statex \textit{// Step D: Link Function Update}
\State Update Sieve coefficients $\hat{\beta}_k$ using generated regressor $X'\hat{\theta}_{k-1}$.
\State \textbf{Output:} Current efficient estimate $\tilde{\theta}_k$ and link function $\hat{F}_k$.
\end{algorithmic}
\end{algorithm}
\end{comment}
\section{Online Learning of Policy Interventions}\label{sec:policy}
So far we have obtained the online estimators for $\theta_0$ and $F_0$. Inference based on random scaling is also provided. Given such estimators, we can provide online estimation and inference for functionals. In this section, we consider an example of functionals of particular interest, the \textit{weighted marginal effect} (WME) defined as
\begin{equation}\label{wme}
\tau_{\rho, 0} = \mathbb{E} \left(\partial F_0(x_0+ X^{\prime}\theta_0)\rho(x_0, X)^{\prime}\theta_0\right),
\end{equation}
where $\rho$ is assumed to be a bounded $p\times 1$ vector of functions that measure policy interventions over individuals with different features. For example, if the first argument of $X$ is the individual wealth level, then a policy targeted at low-income group can be described by function $\rho(x_0, X) = \boldsymbol{1}(x_1 \leq c_{\mathrm{wealth}})$, where $c_{\mathrm{wealth}}$ is a policy-specific threshold of income level for the classification of low-income group.
The target of this section is to construct a sequence of estimators $\bar{\tau}_{1,\rho}, \bar{\tau}_{2,\rho},\cdots$ such that (i) $\bar\tau_{k,\rho}$ can be constructed by sequentially using the first $k$ batches of data with minimal computational efforts, and (ii) the estimation sequence $\bar{\tau}_{1,\rho}, \bar{\tau}_{2,\rho},\cdots$ can be directly used for inference using random scaling, so no plug-in variance needs to be calculated.
Define $\rho_{i, k} = \rho(x_{0,i,k}, X_{i, k})$. We define the following iterative averages: $ \bar{\partial\Psi}_{\rho, 0} = 0$,
\begin{equation}
\bar{\partial\Psi}_{\rho, k} = \frac{k - 1}{k}[\bar{\Psi}_{k-1}]_k + \frac{1}{kB}\sum_{i=1}^{B}\left[\partial \Psi_{J_{k-1}}\left(x_{0,i,k}+X_{i, k}^{\prime}\tilde\theta_k\right)\right]_k\rho_{i,k}^{\prime},
\end{equation}
where for any $m\times p$ matrix $A$, $[A]_N$ is understood to be $[A]_N = (A^{\prime}, 0)^{\prime}$, where the last matrix is $J_N\times p$. We have that
\[
\bar{\partial\Psi}_{\rho, N} = \frac{1}{N}\sum_{k=1}^N \frac{1}{B}\sum_{i=1}^B \left[\partial \Psi_{J_{k-1}}\left(x_{0,i,k}+X_{i, k}^{\prime}\tilde\theta_k\right)\right]_N\rho_{i,k}^{\prime}.
\]
Then we define
\begin{equation}
\bar\tau_{\rho, N} = \bar{\mathcal{B}}_{N}^{\prime}\bar{\partial \Psi}_{\rho, N}\cdot\frac{1}{N}\sum_{k=1}^N \tilde{\theta}_k.
\end{equation}
Note that $\bar{\tau}_{\rho, N}$ can be constructed sequentially.
\begin{theorem}\label{theorem8}
Let all conditions in \autoref{rate-optimal} hold and $\Vert \rho\Vert_{\infty}<\infty$, then we have that
\begin{align}
&\bar\tau_{\rho, N} - \tau_{\rho,0}\nonumber \\
& = \frac{\mathbb{E}(\rho(X)^{\prime}\theta_0)}{N}\sum_{k=1}^N\frac{1}{B}\sum_{i=1}^B \varepsilon_{i,k} + \frac{1}{N}\sum_{k=1}^N \left(\frac{1}{B}\sum_{i=1}^B \partial F_0\left(z_{0,i,k}\right)\rho_{i,k}^{\prime} \theta_0 - \mathbb{E}[\partial F_0\left(z_0\right)\rho(x_0,X)^{\prime} \theta_0]\right) \nonumber \\
& + \left(\mathbb{E} \left(\partial F_0(z_0)\rho(x, X)^{\prime}\right) + \mathbb{E}\left(\partial F_0(z_{0})X^{\prime}\right)+ \mathbb{E}\left[\partial^2 F_0\left(z_0\right)\rho(x,X)^{\prime} \theta_0 X^{\prime}\right]\right)\bar{\tilde\theta}_N + o\left(\frac{1}{\sqrt{N}}\right), \ a.s.
\end{align}
\end{theorem}
\begin{comment}
\subsection{Full-Sample Asymptotic normality of sieve GLS From Chen (2007, subsection 4.2.2)}
Recall that for all the models belonging to the first subclass of the
conditional moment restrictions:
$E\{\rho (Z,\alpha
_{o})|X\}=0$, where $\rho (Z,\alpha )-\rho (Z,\alpha _{o})$ does not depend
on endogenous variables $Y$, we can estimate $\alpha _{o}=(\theta
_{o},F_{o})\in B\times \mathbb{H}$ by the sieve GLS procedure:
\begin{equation*}
\hat{\alpha }_{n}=(\hat{\theta }_{n},\hat{F}_{n})=\arg
\min_{(\theta ,F)\in B\times \mathbb{H}_{n}}\frac{1}{n}\sum_{i=1}^{n}\rho
(Z_{i},\theta ,F)^{\prime }\Sigma (X_{i})^{-1}\rho (Z_{i},\theta ,F),
\end{equation*}
where $\Sigma (X_{i})$ is a positive definite weighting matrix. The optimal weighting matrix $\Sigma _{o}(X_{i})\equiv
Var\{\rho (Z_{i},\alpha_{o})|X_{i}\}$. Applying to our paper's model, we have:
\[
\rho (Z,\alpha _{o})=Y-F_o (x_0 + X'\theta_o )=\varepsilon~,~~\Sigma _{o}(X_{i})=Var\{\varepsilon_i|X_{i}\}
\]
and
\[
\frac{\partial\rho (Z_{i},\alpha _{o})}{\partial \alpha }[\alpha -\alpha_{o}]
=-f_o (x_0 + X'\theta_o )X'[\theta-\theta_o] -[F(x_0 + X'\theta_o )-F_o (x_0 + X'\theta_o )]
\]
where $f_o$ is the derivative of $F_o$. Let
\[
[F(x_0 + X'\theta_o )-F_o (x_0 + X'\theta_o )]=-w(x_0 + X'\theta_o )[\theta-\theta_o]
\]
with $w=(w_{1},...,w_{d_{\theta }})$. Let
\begin{equation*}
D_{w_{j}}(X)=
\frac{\partial \rho (Z,\alpha _{o})}{\partial \theta _{j}}-\frac{\partial
\rho (Z,\alpha _{o})}{\partial F}[w_{j}],
\end{equation*}
\[
\frac{\partial\rho (Z_{i},\alpha _{o})}{\partial \alpha }[\alpha -\alpha_{o}] :=D_{w}(X)[\theta-\theta_o]
=(w(x_0 + X'\theta_o )-f_o (x_0 + X'\theta_o )X')\times [\theta-\theta_o]
\]
with $D_{w}(X)=\left(
D_{w_{1}}(X),...,D_{w_{d_{\theta }}}(X)\right)$ being a $1\times d_{\theta }-$ matrix valued measurable
function of $x_0 + X'\theta_o $.
We now apply Proposition 4.4 of Chen (2007) or Ai and Chen (2003) to obtain root-$n$ asymptotic normality
of the sieve GLS estimator $\hat{\theta }_{n}$. Define the norm $\Vert
\alpha -\alpha_{o}\Vert ^{2}=E\{(\frac{\partial \rho (Z_{i},\alpha _{o})}{
\partial \alpha }[\alpha -\alpha_{o}])\Sigma (X_{i})^{-1}(\frac{\partial
\rho (Z_{i},\alpha _{o})}{\partial \alpha }[\alpha -\alpha_{o}])\}$. For our paper it becomes
\[
\Vert
\alpha -\alpha_{o}\Vert ^{2}=[\theta-\theta_o]'E\{D_w(X)^{\prime }\Sigma (X)^{-1}D_w(X)\}[\theta-\theta_o]
\]
Let $w^{\ast }=(w_{1}^{\ast },...,w_{d_{\theta }}^{\ast })$,
where for $j=1,...,d_{\theta }$, $w_{j}^{\ast }$ solves
\begin{equation*}
E\{D_{w_{j}^{\ast }}(X)^{\prime }\Sigma (X)^{-1}D_{w_{j}^{\ast
}}(X)\}=\inf_{w_{j}}E\{D_{w_{j}}(X)^{\prime }\Sigma (X)^{-1}D_{w_{j}}(X)\}.
\end{equation*}
Let $D_{w^{\ast }}(X)=(w^{\ast }(x_0 + X'\theta_o )-f_o (x_0 + X'\theta_o )X')$.
Let $v_{\theta }^{\ast }=\left( E\{D_{w^{\ast }}(X)^{\prime }\Sigma
(X)^{-1}D_{w^{\ast }}(X)\}\right) ^{-1}\lambda $, $v_{F}^{\ast }=-w^{\ast
}v_{\theta }^{\ast }$ and $v^{\ast }=(v_{\theta }^{\ast },v_{F}^{\ast })$.
\smallskip
\begin{assumption}
(i) $\theta _{o}\in int(B)$; (ii) $E[D_{w^{\ast }}(X)^{\prime }\Sigma
(X)^{-1}D_{w^{\ast }}(X)]$ is positive definite; (iii) there is a $\pi
_{n}v^{\ast }\in \mathcal{A}_{n}$ such that $||\pi _{n}v^{\ast }-v^{\ast
}||\times \Vert \hat{\alpha}_{n}-\alpha _{o}\Vert =o_{P}(n^{-\frac{1}{2}})$.
\end{assumption}
\begin{proposition}
Let $\hat{\alpha}_n $ be the sieve GLS estimate. Suppose Assumption 4.1 - 4.2 hold. Then $n^{\frac{1}{2}}(\hat{\theta }_{n}-\theta
_{o})\stackrel{d}{\longrightarrow }\mathcal{N}(0,V_{1}^{-1}V_{2}V_{1}^{-1})$
where
\[
V_{1}=E[D_{w^{*}}(X)^{\prime }\Sigma (X)^{-1}D_{w^{*}}(X)],\text{ }
V_{2}=E[D_{w^{*}}(X)^{\prime }\Sigma (X)^{-1}\Sigma _{o}(X)\Sigma
(X)^{-1}D_{w^{*}}(X)].
\]
Semiparametric efficient variance corresponds to the case $\Sigma (X)=\Sigma _{o}(X)$, which becomes $n^{\frac{1}{2}}(\hat{\theta }_{n}-\theta
_{o})\stackrel{d}{\longrightarrow }\mathcal{N}(0,V_{o}^{-1})$ with
\[
V_{o}=\inf_wE[D_{w}(X)^{\prime }\Sigma_o (X)^{-1}D_{w}(X)]=E[D_{w_o}(X)^{\prime }\Sigma_o (X)^{-1}D_{w_o}(X)].
\]
Suppose that $\Sigma _{o}(X)=\sigma^2 (x_0 + X'\theta_o )$, then the optimal $w_o$ has a closed form solution:
\[
w_o (x_0 + X'\theta_o ) =f_o (x_0 + X'\theta_o ) E[X' |x_0 + X'\theta_o ]
\]
and
\[
V_{o}=E\left[\frac{f^2_o (x_0 + X'\theta_o )}{\sigma^2 (x_0 + X'\theta_o )}\left(X- E[X |x_0 + X'\theta_o ]\right)\left(X'-E[X'|x_0 + X'\theta_o ]\right)\right].
\]
\end{proposition}
\end{comment}
\section{Monte Carlo Experiments}\label{sec:MC}
\subsection{Simulated Data}
This section reports Monte Carlo simulation results for our algorithm. We consider the binary-choice data-generating process
\begin{align}\label{monte_carlo_model}
Y = \mathbbm{1}\!\left(x_{0} + X^{\prime}\theta_0 - u \ge 0\right),
\end{align}
where $Y$ is the observed response, $(x_0, X')'$ is the regressor vector with $X=(x_1,\dots,x_p)'$, and $u$ is an unobserved shock independent of $(x_0,X)$. The coefficient on $x_0$ is normalized to one.
The true parameter $\theta_0\in\mathbb{R}^p$ is set to
\[
\theta_0=(\theta_{1,0}',\theta_{2,0}',\theta_{3,0}')',
\]
where $\theta_{1,0}\in\mathbb{R}^{0.2p}$ decreases linearly from $1$ to $0$, $\theta_{2,0}\in\mathbb{R}^{0.2p}$ increases linearly from $-1$ to $0$, and $\theta_{3,0}=\mathbf{0}_{0.6p}$. The regressors $x_0,x_1,\ldots,x_p$ are i.i.d. and independent of $u$. We consider two regressor distributions: (i) $x_j\sim N(0,1)$ and (ii) $x_j\sim t(6)/\sqrt{1.5}$ (scaled to have unit variance). We also consider two shock distributions: (i) $u\sim t(2)$ and (ii) a skewed-normal--type distribution generated as
\[
u=\frac{v_1+|v_2|}{\sqrt{2}},\qquad v_1,v_2\ \text{i.i.d. }N(0,1).
\]
We examine two batch sizes, $B\in\{20,40\}$, and fix the total number of online updates at $N=2\times 10^7$. For each configuration, we run 100 independent Monte Carlo replications.
Algorithm \ref{alg:no_online_update} summarizes the online update procedure. Across the 100 replications, we report bias, root mean squared error (RMSE), confidence interval coverage rate (CR), and average confidence interval length. Results are presented in \autoref{table1}--\autoref{table4}. Overall, the online estimator performs well across all designs: coverage at the final update is close to the nominal 0.95 level. Holding the number of updates fixed, larger batch sizes typically yield smaller RMSE and shorter confidence intervals, consistent with a larger effective sample size per update.
Tables 1--4 show that the online estimator delivers accurate point estimates and well-calibrated inference across a range of heavy-tailed and non-Gaussian designs (all reported Bias/RMSE/Length entries are scaled by $10^{-4}$). In the baseline ``hard'' setting with \emph{normal regressors and Cauchy errors} (Table~1), increasing the mini-batch size from $B=20$ to $B=40$ uniformly improves precision: for example, the \emph{average RMSE} falls from $1.930$ to $1.337$ when $p=50$, and from $2.230$ to $1.564$ when $p=100$, while the \emph{average CI length} drops from $10.000$ to $6.699$ ($p=50$) and from $11.000$ to $7.695$ ($p=100$). Coverage is close to nominal throughout, with the \emph{average} CR around $0.944$--$0.948$ across the four $(p,B)$ cells in Table~1.
Comparing across tables highlights how distributional features affect finite-sample difficulty while preserving the same qualitative patterns. Moving from \emph{Cauchy} to \emph{skewed} errors (Tables~2 and~4) typically yields smaller RMSEs and shorter confidence intervals than the corresponding Cauchy-error designs: with normal regressors (Table~2) the average RMSE ranges from $0.787$ to $1.243$ and the average CI length from $3.874$ to $6.220$, both noticeably below Table~1. In this skewed-error setting, the component $\omega_0(1)$ can show mild under-coverage at $B=20$ (about $0.89$--$0.91$), but improves with $B=40$ (about $0.95$--$0.97$), and \emph{average} CR stays near $0.945$ overall. With \emph{heavy-tailed $t(6)$} regressors, the same batch-size and dimension patterns persist: under Cauchy errors (Table~3), average RMSE declines from $1.916$ to $1.372$ ($p=50$) and from $2.191$ to $1.568$ ($p=100$) when increasing $B$, with average CR still around $0.949$--$0.952$. Under skewed errors (Table~4), average RMSE and CI length are again smaller (e.g., Avg RMSE $0.777$ at $(p,B)=(50,40)$ and $0.859$ at $(p,B)=(100,40)$; Avg length $3.911$--$6.217$), with average coverage tightly clustered around $0.946$--$0.951$.
\begin{table}[tbp]
\centering
\caption{Simulation results: Normal Regressors and Cauchy Errors}
\label{table1}
\small
\setlength{\tabcolsep}{6pt}
\begin{tabular}{l c *{4}{c}}
\toprule
& & \multicolumn{2}{c}{$p = 50$} & \multicolumn{2}{c}{$p = 100$} \\
\cmidrule(lr){3-4} \cmidrule(lr){5-6}
&
& $B = 20$ & $B = 40$
& $B = 20$ & $B = 40$ \\
\midrule
\multirow{2}{*}{Kernel RMSE}
& $\theta_0(1)$ & 261.00 & 94.000 & 382.00 & 135.00 \\
& Avg & 116.00 & 50.000 & 161.00 & 62.000 \\
\midrule
\multirow{2}{*}{Bias}
& $\theta_0(1)$ & 0.2384 & 0.2281 & 0.8600 & 0.6481 \\
& Avg & 0.1665 & 0.1034 & 0.2952 & 0.1405 \\
\cmidrule(lr){2-6}
\multirow{2}{*}{RMSE}
& $\theta_0(1)$ & 2.4700 & 1.5330 & 3.2810 & 2.2140\\
& Avg & 1.9300 & 1.3370 & 2.2300 & 1.5640 \\
\midrule
\multirow{2}{*}{CR}
& $\theta_0(1)$ &0.9600 & 0.9800 & 0.9600 &0.9500 \\
& Avg & 0.9456 & 0.9480 &0.9458 & 0.9434 \\
\cmidrule(lr){2-6}
\multirow{2}{*}{Length}
& $\theta_0(1)$ & 13.000 & 9.0510 & 15.000 & 9.9040 \\
& Avg & 10.000 & 6.6990 & 11.000 & 7.6950 \\
\bottomrule
\end{tabular}
\vspace{4pt}
\makebox[\linewidth][l]{
\footnotesize
\parbox{\linewidth}{
\noindent
NOTE: The numbers in Kernel RMSE, Bias, RMSE, and Length are all in magnitude of $10^{-4}$ (For example, the true Kernel RMSE of $\theta_0(1)$ with $p=50$ and $B = 20$ is 0.0261). DGP is $Y = \boldsymbol{1} \left(x_{0} + X^{\prime}\theta_0 - u\geq 0\right)$ with $\theta_0 = (\theta_{1,0}^{\prime},\theta_{2,0}^{\prime}, \theta_{3,0}^{\prime})^{\prime}$, where $\theta_{1,0} = \mathrm{linspace}(1, 0, 0.2p)$, $\theta_{2,0} = -\mathrm{linspace}(1, 0, 0.2p)$, and $\theta_{3,0} = \boldsymbol{0}_{0.6p}$. The Table numbers are Monte Carlo averages over 100 independent simulation replications. All online estimators/algorithms started with zero, with $N_0 = 9\cdot 10^5$, $N_1 = 10^6$, and $N = 4\cdot 10^7$, see Algorithm \ref{alg:no_online_update} for details on $N_0, N_1$ and $N$. $B$ denotes the mini-batch size. The nominal coverage rate is set at 0.95. The Confidence Interval (CR)
is constructed via random scaling in Algorithm \autoref{alg:rs_inference}. $\theta_0(1)$ refers to the first argument of $\theta_0$, whose value is 1. Avg reports the simple average of the summary statistics across all parameters. For example, let $\mathrm{CR}_j$ denote the coverage rate of $\theta_0(j)$ across 100 repetitions, then $\mathrm{Avg} = \frac{1}{p}\sum_{j=1}^p \mathrm{CR}_j$. The above also applies to \autoref{table2}---\autoref{table4}
}
}
\normalsize
\end{table}
\begin{table}[tbp]
\centering
\caption{Simulation results: Normal Regressor and Skewed Normal Errors}
\label{table2}
\small
\setlength{\tabcolsep}{6pt}
\begin{tabular}{l c *{4}{c}}
\toprule
& & \multicolumn{2}{c}{$p = 50$} & \multicolumn{2}{c}{$p = 100$} \\
\cmidrule(lr){3-4} \cmidrule(lr){5-6}
&
& $B = 20$ & $B = 40$
& $B = 20$ & $B = 40$ \\
\midrule
\multirow{2}{*}{Kernel RMSE}
& $\theta_0(1)$ & 117.00 & 54.000 & 188.00 & 61.000 \\
& Avg & 61.000 & 28.000 & 81.000 & 33.000 \\
\midrule
\multirow{2}{*}{Bias}
& $\theta_0(1)$ & 0.0024 & 0.1029 & 0.1812 & 0.0092 \\
& Avg & 0.0807 & 0.0635 & 0.1483 & 0.0063\\
\cmidrule(lr){2-6}
\multirow{2}{*}{RMSE}
& $\theta_0(1)$ & 1.4360 & 1.1020 & 1.7550 &1.1170 \\
& Avg & 1.1020 & 0.7870 & 1.2430 & 0.8530 \\
\midrule
\multirow{2}{*}{CR}
& $\theta_0(1)$ & 0.8900 & 0.9500 & 0.9100 & 0.9700 \\
& Avg & 0.9456 & 0.9424 & 0.9467 & 0.9463 \\
\cmidrule(lr){2-6}
\multirow{2}{*}{Length}
& $\theta_0(1)$ & 6.8000 & 4.9400 & 8.4050 & 6.0200 \\
& Avg & 5.3960 & 3.8740 & 6.2200 & 4.3230\\
\bottomrule
\end{tabular}
\end{table}
\begin{table}[tbp]
\centering
\caption{Simulation results: $t(6)$ Regressors and Cauchy Errors}
\label{table3}
\small
\setlength{\tabcolsep}{6pt}
\begin{tabular}{l c *{4}{c}}
\toprule
& & \multicolumn{2}{c}{$p = 50$} & \multicolumn{2}{c}{$p = 100$} \\
\cmidrule(lr){3-4} \cmidrule(lr){5-6}
&
& $B = 20$ & $B = 40$
& $B = 20$ & $B = 40$ \\
\midrule
\multirow{2}{*}{Kernel RMSE}
& $\theta_0(1)$ & 284.00 & 9.5000 & 428.00 & 135.00 \\
& Avg & 123.00 & 5.1000 & 168.00 & 64.000 \\
\midrule
\multirow{2}{*}{Bias}
& $\theta_0(1)$ & 0.7205 & 0.0547 & 1.5200 & 0.4600 \\
& Avg & 0.2180 & 0.0757 & 0.2840 & 0.1943\\
\cmidrule(lr){2-6}
\multirow{2}{*}{RMSE}
& $\theta_0(1)$ & 2.5970 & 2.0960 & 3.5060 & 2.3910\\
& Avg & 1.9160 & 1.3720 & 2.1910 & 1.5680\\
\midrule
\multirow{2}{*}{CR}
& $\theta_0(1)$ & 0.9400 & 0.9000 & 0.9200 &0.9300 \\
& Avg & 0.9524 &0.9490 & 0.9507 & 0.9488\\
\cmidrule(lr){2-6}
\multirow{2}{*}{Length}
& $\theta_0(1)$ & 13.000 & 9.2580 & 17.000 &11.000 \\
& Avg & 10.000 & 6.8280 & 11.000 & 8.0000\\
\bottomrule
\end{tabular}
\end{table}
\begin{table}[tbp]
\centering
\caption{Simulation results: $t(6)$ Regressors and Skewed Normal Errors}
\label{table4}
\small
\setlength{\tabcolsep}{6pt}
\begin{tabular}{l c *{4}{c}}
\toprule
& & \multicolumn{2}{c}{$p = 50$} & \multicolumn{2}{c}{$p = 100$} \\
\cmidrule(lr){3-4} \cmidrule(lr){5-6}
&
& $B = 20$ & $B = 40$
& $B = 20$ & $B = 40$ \\
\midrule
\multirow{2}{*}{Kernel RMSE}
& $\theta_0(1)$ & 150.00 & 53.000 & 197.00 & 67.000 \\
& Avg & 67.000 & 29.000 & 84.000 & 34.000\\
\midrule
\multirow{2}{*}{Bias}
& $\theta_0(1)$ & 0.1663 & 0.0359 & 0.4241 & 0.0031 \\
& Avg & 0.0821 & 0.0592 & 0.1264 &0.0737 \\
\cmidrule(lr){2-6}
\multirow{2}{*}{RMSE}
& $\theta_0(1)$ & 1.5290 &1.0150 & 1.8730 &1.1110 \\
& Avg & 1.0920 & 0.7770 & 1.2490 & 0.8590 \\
\midrule
\multirow{2}{*}{CR}
& $\theta_0(1)$ & 0.9600 & 0.9100 & 0.9500 &0.9500 \\
& Avg & 0.9510 & 0.9458 & 0.9461 & 0.9476 \\
\cmidrule(lr){2-6}
\multirow{2}{*}{Length}
& $\theta_0(1)$ & 8.2230 & 5.0790 & 8.6050 &5.6940 \\
& Avg & 5.5850 & 3.9110 & 6.2170 & 4.3240\\
\bottomrule
\end{tabular}
\end{table}
\begin{table}[tbp]
\centering
\caption{Comparison with Full Sample Estimation: $n=2.5\times 10^{6}$ (Times in Hours)}
\label{table5}
\small
\setlength{\tabcolsep}{7pt}
\begin{tabular}{l l cc cc}
\toprule
& & \multicolumn{2}{c}{$p = 50$} & \multicolumn{2}{c}{$p = 100$} \\
\cmidrule(lr){3-4} \cmidrule(lr){5-6}
& Method & RMSE & Time & RMSE & Time \\
\midrule
& Full Sample & 0.0252 & 0.2300& 0.0404 & 0.6602 \\
\addlinespace[3pt]
\midrule
\addlinespace[2pt]
\multirow{3}{*}{Online Learning: $B = 20$}
& Phase I & 0.0786 & 0.0032& 0.1310 & 0.0052 \\
& Phase II (One Pass) & 0.0283 & 0.0010 & 0.0616 & 0.0017 \\
& Phase II (20 Passes) & 0.0253 &0.0200& 0.0405 & 0.0337 \\
\addlinespace[4pt]
\multirow{3}{*}{Online Learning: $B = 40$}
& Phase I & 0.0418 & 0.0055 & 0.0587 & 0.0094 \\
& Phase II (One Pass) & 0.0282 &0.0008 & 0.0472 & 0.0014 \\
& Phase II (20 Passes) & 0.0253 & 0.0161 & 0.0404 & 0.0286 \\
\bottomrule
\end{tabular}
\vspace{4pt}
\makebox[\linewidth][l]{
\footnotesize
\parbox{\linewidth}{
\noindent
NOTE: Regressors are standard normally distributed and the error term follows a Cauchy distribution.
$n$ denotes the full sample size. The sieve order is fixed at 30.
Full sample estimation uses the SBGD algorithm \citet{khanetal2021}.
Phase I learning conducts kernel updates with fixed step size (0.01) with $N_0 = 4\times 10^5$ and $N_1 = 5\times 10^5$. Online updating with one pass performs $n/B$ updates, each using a batch of size $B$.
Multiple-pass updating processes the full sample 20 times with random reshuffling after each pass.
Initial points for both full sample estimation and online learning are both fixed at origin point $\boldsymbol{0}_p$. For multiple passes, see Algorithm \ref{alg:multiple pass}.
}
}
\normalsize
\end{table}
\begin{table}[tbp]
\centering
\caption{Comparison with Full Sample Estimation: $n=5\times 10^{6}$ (Times in Hours)}
\label{table5}
\small
\setlength{\tabcolsep}{7pt}
\begin{tabular}{l l cc cc}
\toprule
& & \multicolumn{2}{c}{$p = 50$} & \multicolumn{2}{c}{$p = 100$} \\
\cmidrule(lr){3-4} \cmidrule(lr){5-6}
& Method & RMSE & Time & RMSE & Time \\
\midrule
& Full Sample & 0.0178 & 0.4503 & 0.0254 & 1.3371 \\
\addlinespace[3pt]
\midrule
\addlinespace[2pt]
\multirow{3}{*}{Online Learning: $B = 20$}
& Phase I & 0.0785 & 0.0034 & 0.1376 & 0.0069 \\
& Phase II (One Pass) & 0.0190 & 0.0021 & 0.0442 & 0.0037 \\
& Phase II (20 Passes) & 0.0179 & 0.0422 & 0.0259 & 0.0721 \\
\addlinespace[4pt]
\multirow{3}{*}{Online Learning: $B = 40$}
& Phase I & 0.0345 & 0.0060 & 0.0508 & 0.0104 \\
& Phase II (One Pass) & 0.0185 & 0.0017 & 0.0364 & 0.0031 \\
& Phase II (20 Passes) & 0.0178 & 0.0350 & 0.0259 & 0.0626 \\
\bottomrule
\end{tabular}
\end{table}
\begin{table}[tbp]
\centering
\caption{Average Marginal Effect: Normal Regressors and Cauchy Errors}
\label{functional1}
\small
\setlength{\tabcolsep}{6pt}
\begin{tabular}{c c *{4}{c}}
\toprule
& & \multicolumn{2}{c}{$p = 50$} & \multicolumn{2}{c}{$p = 100$} \\
\cmidrule(lr){3-4} \cmidrule(lr){5-6}
&
& $B = 20$ & $B = 40$
& $B = 20$ & $B = 40$ \\
\midrule
\multirow{2}{*}{Bias}
& $\mathbb{E}\partial F_0(x_{0}+X^{\prime}\theta_0)/\partial x_1 \cdot \theta_{0,1}$ & 0.0443 & 0.0367 & 0.0747 & 0.0225 \\
& Avg & 0.0253 & 0.0133 & 0.0245 & 0.0145 \\
\multirow{2}{*}{RMSE}
& $\mathbb{E}\partial F_0(x_{0}+X^{\prime}\theta_0)/\partial x_1 \cdot \theta_{0,1}$ & 0.1487 & 0.1039 & 0.1586 & 0.1014 \\
& Avg & 0.1574 & 0.1099 & 0.1485 & 0.1055 \\
\multirow{2}{*}{CR}
& $\mathbb{E}\partial F_0(x_{0}+X^{\prime}\theta_0)/\partial x_1 \cdot \theta_{0,1}$ & 0.9700 & 0.9300 & 0.9500 & 0.9400 \\
& Avg & 0.9454 & 0.9474 & 0.9451 & 0.9434 \\
\multirow{2}{*}{Length}
& $\mathbb{E}\partial F_0(x_{0}+X^{\prime}\theta_0)/\partial x_1 \cdot \theta_{0,1}$ & 0.7653& 0.5318 & 0.7676 & 0.5343 \\
& Avg & 0.7795 & 0.5501 & 0.7424 & 0.5232 \\
\bottomrule
\end{tabular}
\vspace{4pt}
\makebox[\linewidth][l]{
\footnotesize
\parbox{\linewidth}{
\noindent
NOTE: Bias, RMSE, and Length are all in $10^{-4}$.
}
}
\normalsize
\end{table}
\begin{comment}
\begin{table}[tbp]
\centering
\caption{Simulation results: $t(8)$ Regressors and Cauchy Random Errors}
\label{table3}
\small
\setlength{\tabcolsep}{6pt}
\begin{tabular}{c c *{4}{c}}
\toprule
& & \multicolumn{2}{c}{$p = 25$} & \multicolumn{2}{c}{$p = 50$} \\
\cmidrule(lr){3-4} \cmidrule(lr){5-6}
&
& $B = 20$ & $B = 40$
& $B = 20$ & $B = 40$ \\
\midrule
\multirow{2}{*}{Bias}
& $\theta_0(1)$ & 0.0001 &0.0001 &0.0006 & 0.0002 \\
& Avg &0.0001 &0.0001 &0.0005 & 0.0002\\
\multirow{2}{*}{RMSE}
& $\theta_0(1)$ &0.0003 &0.0003 &0.0011 &0.0008 \\
& Avg &0.0004 &0.0003&0.0011 &0.0008 \\
\multirow{2}{*}{CR}
& $\theta_0(1)$ & 0.9600& 0.9000&0.9700 &0.9700 \\
& Avg &0.9424 &0.9364 &0.9556 &0.9656 \\
\multirow{2}{*}{Length}
& $\theta_0(1)$ & 0.0018&0.0013 &0.0048 & 0.0038\\
& Avg &0.0019 &0.0014 &0.0050 & 0.0038\\
\bottomrule
\end{tabular}
\vspace{4pt}
\makebox[\linewidth][l]{
\footnotesize
\parbox{\linewidth}{
\noindent
}
}
\normalsize
\end{table}
\begin{table}[tbp]
\centering
\caption{Simulation results: $t(8)$ Regressors and $\chi^2(1)$ Random Errors}
\label{table4}
\small
\setlength{\tabcolsep}{6pt}
\begin{tabular}{c c *{4}{c}}
\toprule
& & \multicolumn{2}{c}{$p = 25$} & \multicolumn{2}{c}{$p = 50$} \\
\cmidrule(lr){3-4} \cmidrule(lr){5-6}
&
& $B = 20$ & $B = 40$
& $B = 20$ & $B = 40$ \\
\midrule
\multirow{2}{*}{Bias}
& $\theta_0(1)$ & 0.0000 & 0.0000 & 0.0001 & 0.0000 \\
& Avg & 0.0000 & 0.0000 & 0.0002 &0.0000\\
\multirow{2}{*}{RMSE}
& $\theta_0(1)$ & 0.0002 & 0.0002 &0.0004 &0.0004\\
& Avg & 0.0002 & 0.0002 & 0.0004 &0.0004\\
\multirow{2}{*}{CR}
& $\theta_0(1)$ & 0.9300 & 0.9600 & 0.9400 &0.9400\\
& Avg & 0.9256 & 0.9344 & 0.9444 &0.9664\\
\multirow{2}{*}{Length}
& $\theta_0(1)$ &0.0011 & 0.0008 &0.0021 &0.0019\\
& Avg & 0.0010 & 0.0008 &0.0021 &0.0019\\
\bottomrule
\end{tabular}
\vspace{4pt}
\makebox[\linewidth][l]{
\footnotesize
\parbox{\linewidth}{
\noindent
}
}
\normalsize
\end{table}
\begin{table}[tbp]
\centering
\caption{Average Marginal Effect: Normal Regressors and Cauchy Errors}
\label{functional1}
\small
\setlength{\tabcolsep}{6pt}
\begin{tabular}{c c *{4}{c}}
\toprule
& & \multicolumn{2}{c}{$p = 25$} & \multicolumn{2}{c}{$p = 50$} \\
\cmidrule(lr){3-4} \cmidrule(lr){5-6}
&
& $B = 20$ & $B = 40$
& $B = 20$ & $B = 40$ \\
\midrule
\multirow{2}{*}{Bias}
& $\theta_0(1)$ & 0.0049 & 0.0099 & 0.0286 & 0.0113 \\
& Avg & 0.0120& 0.0106 & 0.0124 & 0.0114 \\
\multirow{2}{*}{RMSE}
& $\theta_0(1)$ & 0.1659 & 0.0997 & 0.1489 & 0.1138\\
& Avg & 0.1726 & 0.1218 & 0.1522 & 0.1127\\
\multirow{2}{*}{CR}
& $\theta_0(1)$ & 0.9600 & 0.9600 & 0.9400 & 0.9600\\
& Avg & 0.9452 & 0.9404 & 0.9512 & 0.9388\\
\multirow{2}{*}{Length}
& $\theta_0(1)$ & 0.8477 & 0.5836 & 0.7680 &0.5259 \\
& Avg & 0.8489 & 0.5931 & 0.7750 & 0.5406\\
\bottomrule
\end{tabular}
\vspace{4pt}
\makebox[\linewidth][l]{
\footnotesize
\parbox{\linewidth}{
\noindent
NOTE: Bias, RMSE, and Length are all in $10^{-4}$.
}
}
\normalsize
\end{table}
\begin{table}[tbp]
\centering
\caption{Average Marginal Effect: Normal Regressors and $\exp(1)$ Errors}
\label{functional1}
\small
\setlength{\tabcolsep}{6pt}
\begin{tabular}{c c *{4}{c}}
\toprule
& & \multicolumn{2}{c}{$p = 25$} & \multicolumn{2}{c}{$p = 50$} \\
\cmidrule(lr){3-4} \cmidrule(lr){5-6}
&
& $B = 20$ & $B = 40$
& $B = 20$ & $B = 40$ \\
\midrule
\multirow{2}{*}{Bias}
& $\theta_0(1)$ & 0.0219&0.0018 &0.0096 & 0.0004 \\
& Avg &0.0096 & 0.0085& 0.0074& 0.0077 \\
\multirow{2}{*}{RMSE}
& $\theta_0(1)$ &0.1425 &0.1034 &0.1249 & 0.0724 \\
& Avg &0.1311 &0.0913&0.1119 & 0.0795 \\
\multirow{2}{*}{CR}
& $\theta_0(1)$ & 0.9300& 0.8900& 0.9200& 0.9400 \\
& Avg & 0.9540&0.9516 &0.9510 & 0.9508 \\
\multirow{2}{*}{Length}
& $\theta_0(1)$ &0.6764 &0.4359&0.6033 & 0.4027 \\
& Avg &0.6692 &0.4632 & 0.5851& 0.4080 \\
\bottomrule
\end{tabular}
\vspace{4pt}
\makebox[\linewidth][l]{
\footnotesize
\parbox{\linewidth}{
\noindent
NOTE: Bias, RMSE, and Length are all in $10^{-4}$.
}
}
\normalsize
\end{table}
\begin{table}[tbp]
\centering
\caption{Marginal Effect Curve: Normal Regressors and Cauchy Errors}
\label{functional1}
\small
\setlength{\tabcolsep}{6pt}
\begin{tabular}{l *{4}{c}}
\toprule
& \multicolumn{2}{c}{$p = 25$} & \multicolumn{2}{c}{$p = 50$} \\
\cmidrule(lr){2-3} \cmidrule(lr){4-5}
& $B = 20$ & $B = 40$
& $B = 20$ & $B = 40$ \\
\midrule
Bias & 0.3633 & 0.1389 & 0.6811 &0.3783 \\
RMSE & 1.5909 & 0.9958 & 1.8393 & 1.3115 \\
CR & 0.9295 & 0.9580 & 0.9395 & 0.9300 \\
Length & 6.4902 & 0.5059 & 7.8522 & 5.1855 \\
\bottomrule
\end{tabular}
\vspace{4pt}
\makebox[\linewidth][l]{
\footnotesize
\parbox{\linewidth}{
\noindent
NOTE: Bias, RMSE, and Length are all in $10^{-4}$.
}
}
\normalsize
\end{table}
\begin{table}[tbp]
\centering
\caption{Comparison with Full Sample Estimation: $n=2\times 10^{6}$ (Times in Hour)}
\label{table5}
\small
\setlength{\tabcolsep}{6pt}
\begin{tabular}{l *{4}{c}}
\toprule
& \multicolumn{2}{c}{$p = 25$} & \multicolumn{2}{c}{$p = 50$} \\
\cmidrule(lr){2-3} \cmidrule(lr){4-5}
& RMSE & Time & RMSE & Time \\
\midrule
Warm Start & 0.1354 & --- & 0.2117 &--- \\
Full Sample & 0.0221 & 0.3483 & 0.0411 & 0.7423 \\
Online Kernel Update: $B = 20$, One Pass & 0.0881 &0.0030 & 0.1432&0.0034 \\
Online Kernel Update: $B = 20$, Multiple Pass &0.0643 &0.1566 & 0.1070&0.1768 \\
Online Orthogonalized Update: $B = 20$, One Pass & 0.0386 & 0.0002 & 0.0979 & 0.0004 \\
Online Orthogonalized Update: $B = 20$, Multiple Passes & 0.0218 & 0.0213 & 0.0424 & 0.0353 \\
Online Kernel Update: $B = 40$, One Pass &0.0919&0.0060 & 0.1527&0.0067 \\
Online Kernel Update: $B = 40$, Multiple Pass &0.0693 & 0.3090& 0.1091& 0.3344\\
Online Orthogonalized Update: $B = 40$, One Pass & 0.0464 & 0.0002 & 0.0981 & 0.0003 \\
Online Orthogonalized Update: $B = 40$, Multiple Passes & 0.0215 & 0.0172 & 0.0455 & 0.0295 \\
\bottomrule
\end{tabular}
\vspace{4pt}
\makebox[\linewidth][l]{
\footnotesize
\parbox{\linewidth}{
\noindent
NOTE: Regressors are standard normally distributed, and error term is Cauchy distributed. $n$ is full sample size. Order of sieve functions is fixed at 20. Full sample estimation uses SBGD algorithm. Online update with one pass updates the parameter with $n/B$ times, each update uses a batch of size $B$. Online update with multiple passes learns the full sample 50 times, in which random reshuffling was used after each one pass of the data. Finally, initial points are generated by SBGD algorithm based on an independent small sample of size 50,000. Simulations are independently repeated 20 times. The same setup applies to \autoref{table6}.
}
}
\normalsize
\end{table}
\begin{table}[tbp]
\centering
\caption{Comparison with Full Sample Estimation: $n=4\times 10^{6}$ (Times in Hour)}
\label{table6}
\small
\setlength{\tabcolsep}{6pt}
\begin{tabular}{l *{4}{c}}
\toprule
& \multicolumn{2}{c}{$p = 25$} & \multicolumn{2}{c}{$p = 50$} \\
\cmidrule(lr){2-3} \cmidrule(lr){4-5}
& RMSE & Time & RMSE & Time \\
\midrule
Warm Start & 0.1354 & ---- & 0.2117 & ---\\
Full Sample &0.0153 & 0.6535 & 0.0265 & 2.9554 \\
Online Kernel Update: $B = 20$, One Pass & 0.0818 &0.0062 &0.1292 &0.0089 \\
Online Kernel Update: $B = 20$, Multiple Passes & 0.0585& 0.3286&0.1023 &0.5580 \\
Online Orthogonalized Update: $B = 20$, One Pass &0.0294 &0.0017 &0.0793 & 0.0034 \\
Online Orthogonalized Update: $B = 20$, Multiple Passes &0.0151 & 0.0617 & 0.0273 & 0.1991\\
Online Kernel Update: $B = 40$, One Pass & 0.0856 &0.0124 & 0.1390& 0.0147\\
Online Kernel Update: $B = 40$, Multiple Passes &0.0649 &0.6206 &0.1052 & 0.8109\\
Online Orthogonalized Update: $B = 40$, One Pass & 0.0367 & 0.0010 & 0.0846 & 0.0007 \\
Online Orthogonalized Update: $B = 40$, Multiple Passes & 0.0150 &0.0460 & 0.0288 & 0.2484 \\
\bottomrule
\end{tabular}
\end{table}
\end{comment}
\subsection{Data Experiments}
\begin{figure}[t!]
\includegraphics[
width=1.0\linewidth,
trim=1.8cm 4.5cm 1.5cm 4cm,
clip
]{empirical1.pdf}
\caption{Empirical Results: Online Learning Updates and PR Averages}
\label{fig:empirical}
\bigskip{}
\footnotesize {NOTE: Batch size $B = 20$. We choose $N_0 = 4.5\times 10^6$, $N_1 = 5\times 10^6$, and $N = 2\cdot 10^7$. So the first learning phase (warm-start phase) takes 5 million updates, where the PR average is calculated at the last $5\cdot 10^5$ updates. The second learning phase takes 1.5 millions updates. The black solid line corresponds to PR average estimator, red solid line corresponds to $\hat\theta_k$ in the first learning phase and $\tilde\theta_k$ in the second phase. Both estimators are constructed using Algorithm \ref{alg:no_online_update}, where the starting point is origin point $\boldsymbol{0}_p$. For warm-start learning phase, 4-th order Epanechnikov kernel function is used with $h_k = c_k \cdot (B\cdot k)^{-1/7}$. For second learning phase, we choose Chebyshev polynomials, sieve order $J_k = J_0 = 30$, $\xi_k = (k - N_1)^{-0.51}$, and $\eta_k = (k - N_1)^{-0.61}$. }
\end{figure}
\begin{figure}[t!]
\includegraphics[
width=1.0\linewidth,
trim=1.8cm 4.0cm 1.5cm 4cm,
clip
]{empirical2.pdf}
\caption{Empirical Results: Confidence Band based on Random Scaling}
\label{fig:empirical}
\bigskip{}
\end{figure}
This section illustrates the empirical applicability of our proposed semiparametric online learning method by analyzing the data set of \citet*{helpman2008estimating}. The original data contain 248,060 observations with 333 covariates including 10 regressors, 157 exporter country dummies, 157 importer country dummies, and 9 year dummies\footnote{The data contain 158 countries and 10 years. When construct exporter country, importer country, and year dummies, we leave out one country and one year to ensure identification of the model.}, which features high dimensionality. In applications, the trading data can be collected at various frequency such as monthly, quarterly, or annually. To mimic such availability of stream data, we repeatedly draw random batches of the data points from the original data set with replacement, and take each random draw as the new data points. Following \citet{helpman2008estimating}, we use the data stream to learn the following model describing bilateral trading relationship \begin{align}
P\left(\left.T_{jl}=1\right|\text{observed variables}\right) = F_0\left(\theta_{I,0} + \theta_{E,0} + d_{jl} + X^{\prime}_{jl}\theta_0
\right),\label{empirical_model}
\end{align}
where $T_{jl}$ is an indicator of whether country $l$ exports to country $j$, $\theta_{E,0}$ is the exporter fixed effect of the $l$-th country, $\theta_{I,0}$ is the importer fixed effect of the $j$-th country, $d_{jl}$ is the natural logarithm of the geographic distance between countries $j$ and $l$, and $X_{jl}$ is a vector of covariates that describe the variable country-pair fixed trade cost.
\citet{helpman2008estimating} consider estimating (\ref{empirical_model}) using the original data set based on Probit specification, which amounts to specifying $F_0$ to be the CDF of standard normal distribution. In our application of semiparametric online estimation, we leave $F_0$ unspecified, which allows for non-normal random shocks in the trading relationships.
\paragraph{Empirical learning trajectories and confidence bands.}
Figures~1--2 summarize how the online estimates evolve over the update index for a set of representative bilateral-trade covariates (Land Border, Island, Landlock, Legal, Language, Colonial Ties, Currency Union, FTA, and Religion). In Figure~1, the red path reports the sequential updated parameter estimates, while the black path reports the corresponding Polyak--Ruppert (PR) averages. The updated estimates move quickly early on and can be visibly noisy, especially before and around the transition from the warm-start phase; after the PR averaging begins (the black series appears once averaging starts), the PR-average trajectories become stable and the red iterates fluctuate around them. Across the panels, the limiting levels are negative for \emph{Land Border} (around $-0.10$) and positive for the remaining indicators, with noticeably larger magnitudes for \emph{Island} (about $2.4$) and \emph{Landlock} (about $3.0$) than for the other covariates.
Figure~2 focuses on the PR-average trajectories in the post--warm-start region and overlays 95\% confidence bands constructed via random scaling. The bands are widest near the phase transition and then tighten as the number of updates increases, remaining centered on the same stabilized PR-average paths. Visually, this yields relatively tight uncertainty envelopes around the eventual PR-average levels for each coefficient, with especially narrow bands for the smaller-magnitude trade-cost indicators (e.g., Legal and Colonial Ties) once the algorithm has accumulated many updates.
\section{Conclusion}
This paper develops a two-phase online learning paradigm for semiparametric monotone index models that is designed for streaming environments where continuously re-estimating on the full sample (or even storing it) is infeasible.
The warm-start phase uses a globally stable update rule that consistently learns the finite-dimensional index parameter from arbitrary initialization, and then transitions to a rate-optimal phase that jointly updates the index parameter using a Neyman-orthogonalized score and the unknown monotone link using an online sieve procedure.
In this second phase, both components attain optimal convergence rates, delivering a practical semiparametric toolkit for real-time estimation.
A key byproduct of the procedure is a sequence of parameter updates—learning trajectories—that can be used for online inference via random scaling with essentially no additional nonparametric estimation burden.
The same trajectories also support online estimation and inference for policy-relevant functionals that depend on both the parametric and nonparametric components.
Monte Carlo experiments show good finite-sample performance with coverage rates close to nominal,
and an empirical illustration using the Helpman et al. (2008) trade dataset demonstrates feasibility in a high-dimensional setting while leaving the link function unspecified.
Finally, the framework naturally invites extensions—e.g., to online semiparametric models with sample selection—highlighting a broader agenda of bringing semiparametric econometrics into online learning environments.
\newpage
\bibliography{overall_1}
\newpage