Online Learning in Semiparametric Econometric Models

Online Learning in Semiparametric Econometric Models

Introduction

Consider a generic monotone index model

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 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 han1987non,powell1989semiparametric,ichimura1993semiparametric,klein1993efficient,neweyetal,ahn2018simple,fan2020rank,khanetal2021 and also powellhandbook and references therein. Notably, all the above methods are designed by nature to be 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)) 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 robbins1951stochastic. When a gradient is used for update, online learning can be viewed as stochastic gradient descent (SGD, bottou2007tradeoffs,bottou2010large,bottou2018optimization), which solves the optimization problem in finite-dimensional spaces or infinite-dimensional spaces such as Hilbert space chen2002asymptotic,godichon2019lp. A large body of literature studies the properties of SGD algorithms for loss function with strong convexity toulis2017asymptotic,jentzen2021strong,sebbouh2021almost, local strong convexity godichon2019lp, or nonconvexity 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 zhang2022sieve or kernel density estimation/regression huang2013recursive. More recently, SA methods have been applied to analyze econometric models such as GMM models 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 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 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 two-phase online learning paradigm for models in ((ref)). 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 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 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) we provide detailed instructions via algorithms that show how one can implement our procedures for learning $\theta_0$. Also, Section (ref) extends these to procedures to obtain marginal effects. Finally, section (ref) 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$.

Learning Phase I: Warm-Start Learning

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$ i.i.d. realizations of $(x_0, X, Y)$ from model ((ref)) denoted as $(x_{0,i,k}, X_{i,k},Y_{i,k})$, $i = 1, 2, \cdots, B$. To ease exposition, we denote

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 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:

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)).

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

Typically, $\bar{\theta}_N$ is the sample average of the first $N$ iterates obtained based on ((ref)). Note that evaluating $\bar{\theta}_N$ based on algorithm ((ref)) requires storing the full trajectory of updates. To save memory, we consider the following online update for evaluating $\bar{\theta}_N$

In this case, only the most recent PR average estimator needs to be stored in the memory.

Preliminaries and Intuitive Discussion

Define

and

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.

Given (ref) and (ref), we have the following lemma that describes the properties of $\Phi_k(\theta,W_k), \Phi_k(\theta),$ and $\Phi(\theta)$.

Given (ref), we can briefly discuss the intuition of our algorithm. Note that ((ref)) leads to

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 (ref). So we can regard the update ((ref)) as being mainly driven by term (I). Define

then for any $\theta$, we obviously have that

and

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)) 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.

Statistical Properties of $\hat\theta_k$ and $\bar\theta_N$

This section formally studies the statistical properties of our proposed algorithm ((ref)) and ((ref)). We first impose some further conditions.

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.

The lower bound of the smallest eigenvalue as a function of $\theta$ in (ref) is crucial because it guarantees that term (I) in ((ref)) is everywhere a strict contraction mapping. Note that such result also guarantees point identification of $\theta_0$.

Based on (ref), we are ready to state the first theorem regarding the convergence of the online estimator $\hat\theta_k$.

(ref) 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 (ref), 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.

Under (ref), all the requirements on $\gamma_k$ and $h_k$ in (ref) 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

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.

(ref) provides the sharp rate for our online estimator as well as its asymptotic distribution. Using (ref) provided in Section (ref) 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)) implies that $\Vert \Delta\hat\theta_k\Vert = O(\sqrt{N^{-\alpha_{\gamma} + \alpha_{h} }\log(\log(N))})$ a.s.. Since (ref) is fundamental in proving the following theorems and is of interest on its own side, we provide detailed development in Section (ref).

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)). We have the following 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 (ref), $\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)). As we can see from (ref), 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) introduces a second learning phase that yields a $1/ \sqrt N$ -consistent estimator of $\theta_0.$

Online Learning of $F_0$

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.$

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 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

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,

Following chen2015optimal, we define the $L_2$-projection of function $F_0$ onto $\mathcal{S}_J$ as

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

We impose the following condition regarding the supreme norm of operator $\mathbb{P}_J(\cdot)$.

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

and the online estimator of $F_0(\cdot)$ is given by

Based on ((ref)), 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

While the above expression is difficult to implement in practice. Note that based on ((ref)), 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

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

For any $N$, we have

It's straightforward to verify that

We point out that formulation ((ref)) will be used to study the statistical properties of the PR average estimator $\bar{F}_N(\cdot)$.

Statistical Properties

We impose the following conditions regarding the choices of the learning rate, sieve functions, and sieve orders.

Define set

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)). We have the following result.

The rate in (ref) 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), PR average will lead to optimal rate of our nonparametric estimator. We first introduce another technical condition.

Based on the above conditions, we can provide a more refined convergence rate for $\Delta\hat{\mathcal{B}}_N$.

For each $N$, define

We have the following result.

(ref) 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., 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.

Learning Phase II: Rate-Optimal Learning

This learning phase simultaneously updates estimators of $\theta_0$ and $F_0$, and produces estimators with optimal convergence rate. In Section (ref), 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), we combine the orthogonalized update with the learning of $F_0$ based on online method of sieves proposed in Section (ref). 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), we discuss using trajectories to conduct online inference based on random scaling method.

The Orthogonalized Update

In Section (ref), 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 (ref), 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

and the associated update

where $\tilde\theta_0$ is the initial point and is given, and $\xi_k>0$ is the learning rate. The score function ((ref)) and corresponding update ((ref)) 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), but such estimator converges at a rate slower than $1/\sqrt{k}$ according to (ref). As a result, the estimation error of the sieve estimator will remain and contaminate the update ((ref)), 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)), 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

Using ((ref)) 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), 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

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

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.

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$

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)) as follows

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.

Based on the above conditions, we now provide the statistical properties of of $\tilde{\theta}_k$. Define

and

we have the following theorem.

Rate-Optimal Learning

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.

We have the following theorem for the second learning phase.

Online Inference based on Random Scaling

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 (ref) ((ref)), we can apply a computationally super efficient algorithm (see Algorithm (ref)) 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. 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 (ref). We summarize the inference for our online semiparametric estimator via random scaling in Algorithm (ref) in Section (ref) for easy reference.

Practical Implementations

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). Due to the global contraction mapping properties indicated in (ref) and (ref), 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

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

Note that ((ref)) and ((ref)) 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 (ref).

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)). 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

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)), 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

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 (ref), 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

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 lee2022fast. The algorithm is provided in Algorithm (ref).

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$.

Online Learning of Policy Interventions

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 weighted marginal effect (WME) defined as

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$,

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

Note that $\bar{\tau}_{\rho, N}$ can be constructed sequentially.

Monte Carlo Experiments

Simulated Data

This section reports Monte Carlo simulation results for our algorithm. We consider the binary-choice data-generating process

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) 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 (ref)--(ref). 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 normal regressors and Cauchy errors (Table 1), increasing the mini-batch size from $B=20$ to $B=40$ uniformly improves precision: for example, the average RMSE falls from $1.930$ to $1.337$ when $p=50$, and from $2.230$ to $1.564$ when $p=100$, while the 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 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 Cauchy to 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 average CR stays near $0.945$ overall. With 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$.

Data Experiments

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 helpman2008estimating, we use the data stream to learn the following model describing bilateral trading relationship

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.

helpman2008estimating consider estimating ((ref)) 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 Land Border (around $-0.10$) and positive for the remaining indicators, with noticeably larger magnitudes for Island (about $2.4$) and 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.

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.