Stochastic Gradient Descent in Continuous Time

Introduction

This paper develops a statistical learning algorithm for continuous-time models, which are common in science, engineering, and finance. We study its theoretical convergence properties as well as its computational performance in a number of benchmark problems. Although the method is broadly applicable, this paper mainly focuses on applications in finance. Given a continuous stream of data, stochastic gradient descent in continuous time (SGDCT) can estimate unknown parameters or functions in stochastic differential equation (SDE) models for stocks, bonds, interest rates, and financial derivatives. The statistical learning algorithm can also be used for the optimization of high-dimensional continuous-time models, such as American options. High-dimensional American options have been a longstanding computational challenge in finance. SGDCT is able to accurately solve American options even in 100 dimensions.

Batch optimization for the statistical estimation of continuous-time models can be impractical for large datasets where observations occur over a long period of time. Batch optimization takes a sequence of descent steps for the model error for the entire observed data path. Since each descent step is for the model error for the entire observed data path, batch optimization is slow (sometimes impractically slow) for long periods of time or models which are computationally costly to evaluate (e.g., partial differential equations). Typical existing approaches in the financial statistics literature use batch optimization.

SGDCT provides a computationally efficient method for statistical learning over long time periods and for complex models. SGDCT continuously follows a (noisy) descent direction along the path of the observation; this results in much more rapid convergence. Parameters are updated online in continuous time, with the parameter updates $\theta_{t}$ satisfying a stochastic differential equation. We prove that $\lim_{t\rightarrow\infty}\nabla\bar{g}(\theta_{t})=0$ where $\bar{g}$ is a natural objective function for the estimation of the continuous-time dynamics.

The stochastic gradient descent update in continuous time follows the SDE:

where $\nabla_{\theta}f(X_{t};\theta_{t})$ is matrix valued and $\alpha_{t}$ is the learning rate. The parameter update (1.2) can be used for both statistical estimation given previously observed data as well as online learning (i.e., statistical estimation in real-time as data becomes available). SGDCT will still converge if $\sigma\sigma^{\top}$ in (1.2) is replaced by the identity matrix $I$.

We assume that $X_{t}$ is sufficiently ergodic (to be concretely specified later in the paper) and that it has some well-behaved $\pi(dx)$ as its unique invariant measure. As a general notation, if $h(x,\theta)$ is a generic $L^{1}(\pi)$ function, then we define its average over $\pi(dx)$ to be

The gradient $\nabla_{\theta}g(X_{t},\theta)$ cannot be evaluated since $f^{\ast}(x)$ is unknown. However, $dX_{t}=f^{\ast}(X_{t})dt+\sigma dW_{t}$ is a noisy estimate of $f^{\ast}(x)dt$, which leads to the algorithm (1.2). SGDCT follows a noisy descent direction along a continuous stream of data produced by $X_{t}$.

Heuristically, it is expected that $\theta_{t}$ will tend towards the minimum of the function $\bar{g}(\theta)=\int_{\mathcal{X}}g(x,\theta)\pi(dx)$. The data $X_{t}$ will be correlated over time, which complicates the mathematical analysis. This differs from the standard discrete-time version of stochastic gradient descent where the the data is usually considered to be i.i.d. at every step.

In this paper we show that if $\alpha_{t}$ is appropriately chosen then $\nabla\bar{g}(\theta_{t})\rightarrow 0$ as $t\rightarrow\infty$ with probability 1 (see Theorem 2.4). Results like this have been previously derived for stochastic gradient descent in discrete time; see and . proves convergence in the absence of the $X$ term. proves convergence of stochastic gradient descent in discrete time with the $X$ process but requires stronger conditions than .

Although stochastic gradient descent for discrete time has been extensively studied, stochastic gradient descent in continuous time has received relatively little attention. We refer readers to and for a thorough review of the very large literature on stochastic gradient descent. There are also many algorithms which modify traditional stochastic gradient descent (stochastic gradient descent with momentum, Adagrad, RMSprop, etc.). For a review of these variants of stochastic gradient descent, see . We mention below the prior work which is most relevant to our paper.

Our approach and assumptions required for convergence are most similar to , who prove convergence of discrete-time stochastic gradient descent in the absence of the $X$ process. The presence of the $X$ process is essential for considering a wide range of problems in continuous time, and showing convergence with its presence is considerably more difficult. The $X$ term introduces correlation across times, and this correlation does not disappear as time tends to infinity. This makes it challenging to prove convergence in the continuous-time case. In order to prove convergence, we use an appropriate Poisson equation associated with $X$ to describe the evolution of the parameters for large times.

proves, in a setting different than ours, convergence in $L^{2}$ of projected stochastic gradient descent in discrete time for convex functions. In projected gradient descent, the parameters are projected back into an a priori chosen compact set. Therefore, the algorithm cannot hope to reach the minimum if the minimum is located outside of the chosen compact set. Of course, the compact set can be chosen to be very large for practical purposes. Our paper considers unconstrained stochastic gradient descent in continuous time and proves the almost sure convergence $\nabla\bar{g}(\theta_{t})\rightarrow 0$ as $t\rightarrow\infty$ taking into account the $X$ component as well. We do not assume any stability conditions on $X$ (except that it is ergodic with a unique invariant measure).

Another approach for proving convergence of discrete-time stochastic gradient descent is to show that the algorithm converges to the solution of an ODE which itself converges to a limiting point. This is the approach of . See also . This method, sometimes called the “ODE method”, requires the assumption that the iterates (i.e., the model parameters which are being learned) remain in a bounded set with probability one. It is unclear whether the ODE method of proof can be successfully used to show convergence for a continuous-time stochastic gradient descent scheme. In this paper we follow a potentially more straightforward method of proof by analyzing the speed of convergence to equilibrium with an appropriately chosen Poisson type of equation.

studies continuous-time stochastic mirror descent in a setting different than ours. In the framework of , the objective function is known. In this paper, we consider the statistical estimation of the unknown dynamics of a random process (i.e. the $X$ process satisfying (1.1)).

Statisticians and financial engineers have actively studied parameter estimation of SDEs, although typically not with statistical learning or machine learning approaches. The likelihood function will usually be calculated from the entire observed path of $X$ (i.e., batch optimization) and then maximized to find the maximum likelihood estimator (MLE). Unlike in this paper, the actual optimization procedure to maximize the likelihood function is often not analyzed.

Some relevant publications in the financial statistics literature include , , , and . derives the likelihood function for continuously observed $X$. The MLE can be calculated via batch optimization. and consider the case where $X$ is discretely observed and calculate MLEs via a batch optimization approach. estimates parameters by a Bayesian approach. Readers are referred to for thorough reviews of classical statistical inference methods for stochastic differential equations.

2 Applications of SGDCT

Continuous-time models are especially common in finance. Given a continuous stream of data, the stochastic gradient descent algorithm can be used to estimate unknown parameters or functions in SDE models for stocks, bonds, interest rates, and financial derivatives. Numerical analysis of SGDCT for two common financial models is included in Sections 5.1, 5.2, and 5.5. The first is the well-known Ornstein-Uhlenbeck (OU) process (for examples in finance, see , , , and ). The second is the multidimensional CIR process which is a common model for interest rates (for examples in finance, see , , , , and ).

Scientific and engineering models are also typically in continuous-time. There are often coefficients or functions in these models which are uncertain or unknown; stochastic gradient descent can be used to learn these model parameters from data. In Section 5, we study the numerical performance for two example applications: Burger’s equation and the classic reinforcement learning problem of balancing a pole on a moving cart. Burger’s equation is a widely used nonlinear partial differential equation which is important to fluid mechanics, acoustics, and aerodynamics.

A natural question is why use SGDCT versus a straightforward approach which (1) discretizes the continuous-time dynamics and then (2) applies traditional stochastic gradient descent. For some of the same reasons that scientific models have been largely developed in continuous time, it can be advantageous to develop continuous-time statistical learning for continuous-time models.

SGDCT allows for the application of numerical schemes of choice to the theoretically correct statistical learning equation for continuous-time models. This can lead to more accurate and more computationally efficient parameter updates. Numerical schemes are always applied to continuous-time dynamics and different numerical schemes may have different properties for different continuous-time models. A priori performing a discretization to the system dynamics and then applying a traditional discrete-time stochastic gradient descent scheme can result in a loss of accuracy. For example, there is no guarantee that (1) using a higher-order accurate scheme to discretize the system dynamics and then (2) applying traditional stochastic gradient descent will produce a statistical learning scheme which is higher-order accurate in time. Hence, it makes sense to first develop the continuous-time statistical learning equation, and then apply the higher-order accurate numerical scheme.

Besides model estimation, SGDCT can be used to solve continuous-time optimization problems, such as American options. We combine SGDCT with a deep neural network to solve American options in up to $100$ dimensions (see Section 6). An alternative approach would be to discretize the dynamics and then use the Q-learning algorithm (traditional stochastic gradient descent applied to an approximation of the discrete HJB equation). However, Q-learning is biased while SGDCT is unbiased. Furthermore, in SDE models with Brownian motions, the Q-learning algorithm can blow up as the time step size $\Delta$ becomes small; see Section 6 for details.

The convergence issue with Q-learning highlights the importance of studying continuous-time algorithms for continuous-time models. It is of interest to show that (1) a discrete-time scheme converges to an appropriate continuous-time scheme as $\Delta\rightarrow 0$ and (2) the continuous-time scheme converges to the correct estimate as $t\rightarrow\infty$. These are important questions since any discrete scheme for a continuous-time model incurs some error proportional to $\Delta$, and therefore $\Delta$ must be decreased to reduce error. It is also important to note that in some cases, such as Q-learning, computationally expensive terms in the discrete algorithm (such as expectations over high-dimensional spaces) may become much simpler expressions in the continuous-time scheme (differential operators).

3 Organization of Paper

The paper is organized into five main sections. Section 2 presents the assumption and the main theorem. In Section 3 we prove the main result of this paper for the convergence of continuous-time stochastic gradient descent. The extension of the stochastic gradient descent algorithm to the case of a variable diffusion coefficient function is described in Section 4. Section 5 provides numerical analysis of SGDCT for model estimation in several applications. Section 6 discusses SGDCT for solving continuous-time optimization problems, particularly focusing on American options.

Assumptions and Main Result

Before presenting the main result of this paper, Theorem 2.4, let us elaborate on the standing assumptions. In regards to the learning rate $\alpha_{t}$ the standing assumption is

Assume that $\int_{0}^{\infty}\alpha_{t}dt=\infty$, $\int_{0}^{\infty}\alpha^{2}_{t}dt<\infty$, $\int_{0}^{\infty}|\alpha_{s}^{\prime}|ds<\infty$ and that there is a $p>0$ such that $\lim_{t\rightarrow\infty}\alpha_{t}^{2}t^{1/2+2p}$=0.

A standard choice for $\alpha_{t}$ that satisfies Condition 2.1 is $\alpha_{t}=\frac{1}{C+t}$ for some constant $0<C<\infty$. Notice that the condition $\int_{0}^{\infty}|\alpha_{s}^{\prime}|ds<\infty$ follows immediately from the other two restrictions for the learning rate if it is chosen to be a monotonic function of $t$.

Let us next discuss the assumptions that we impose on $\sigma$, $f^{\ast}(x)$ and $f(x,\theta)$. Condition 2.2 guarantees uniqueness and existence of an invariant measure for the $X$ process.

We assume that $\sigma\sigma^{\top}$ is non-degenerate bounded diffusion matrix and $\lim_{|x|\rightarrow\infty}f^{\ast}(x)\cdot x=-\infty$

Moreover, there exists $K>0$ and $q>0$ such that

The function $f^{\ast}(x)$ is $C^{2+\alpha}(\mathcal{X})$ with $\alpha\in(0,1)$. Namely, it has two derivatives in $x$, with all partial derivatives being Hölder continuous, with exponent $\alpha$, with respect to $x$.

Condition 2.3 allows one to control the ergodic behavior of the X process. As will be seen from the proof of the main convergence result Theorem 2.4, one needs to control terms of the form $\int_{0}^{t}\alpha_{t}(\nabla\bar{g}(\theta_{s})-g(X_{s},\theta_{s}))ds$. Due to ergodicity of the $X$ process one expects that such terms are small in magnitude and go to zero as $t\rightarrow\infty$. However, the speed at which they go to zero is what matters here. We treat such terms by rewriting them equivalently using appropriate Poisson type partial differential equations (PDE). Condition 2.3 guarantees that these Poisson equations have unique solutions that do not grow faster than polynomially in the $x$ variable (see Theorem A.1 in Appendix A).

The main result of this paper is Theorem 2.4.

Assume that Conditions 2.1, 2.2 and 2.3 hold. Then we have that

Proof of Theorem 2.4

We proceed in a spirit similar to that of . However, apart from continuous versus discrete dynamics, one of the main challenges of the proof here is the presence of the ergodic $X$ process. Let us consider an arbitrarily given $\kappa>0$ and $\lambda=\lambda(\kappa)>0$ to be chosen. Then set $\sigma_{0}=0$ and consider the cycles of random times

The purpose of these random times is to control the periods of time where $\|\nabla\bar{g}(\theta_{\cdot})\|$ is close to zero and away from zero. Let us next define the random time intervals $J_{k}=[\sigma_{k-1},\tau_{k})$ and $I_{k}=[\tau_{k},\sigma_{k})$. Notice that for every $t\in J_{k}$ we have $\|\nabla\bar{g}(\theta_{t})\|<\kappa$.

Let us next consider some $\eta>0$ sufficiently small to be chosen later on and set $\sigma_{k,\eta}=\sigma_{k}+\eta$. Lemma 3.1 is crucial for the proof of Theorem 2.4.

Assume that Conditions 2.1, 2.2 and 2.3 hold. Let us set

The idea is to use Theorem A.1 in order to get an equivalent expression for the term $\Gamma_{k,\eta}$ that we seek to control.

Let us consider the function $G(x,\theta)=\nabla_{\theta}g(x,\theta)-\nabla_{\theta}\bar{g}(\theta)$. Notice that by definition and due to Condition 2.3, the function $G(x,\theta)$ satisfies the centering condition (A.1) of Theorem A.1 componentwise. So, the Poisson equation (A.2) will have a unique smooth solution, denoted by $v(x,\theta)$ that grows at most polynomially in $x$. Let us apply Itô formula to the vector valued function $u(t,x,\theta)=\alpha_{t}v(x,\theta)$. Doing so, we get for $i=1,\cdots,n$

where $\mathcal{L}_{x}$ and $\mathcal{L}_{\theta}$ denote the infinitesimal generators for processes $X$ and $\theta$ respectively.

Recall now that $v(x,\theta)$ is the solution to the given Poisson equation and that $u(s,x,\theta)=\alpha_{s}v(x,\theta)$. Using these facts and rearranging the previous Itô formula, we get in vector notation

The next step is to treat each term on the right hand side of (3) separately. For this purpose, let us first set

By Theorem A.1 and Proposition 2 of there is some $0<K<\infty$ (that may change from line to line below) and $0<q<\infty$ such that for $t$ large enough

By Condition 2.1 let us consider $p>0$ such that $\lim_{t\rightarrow\infty}\alpha_{t}^{2}t^{1/2+2p}=0$ and for any $\delta\in(0,p)$ define the event $A_{t,\delta}=\left\{J_{t}^{(1)}\geq t^{\delta-p}\right\}$. Then we have for $t$ large enough such that $\alpha_{t}^{2}t^{1/2+2p}\leq 1$

Therefore, by Borel-Cantelli lemma we have that for every $\delta\in(0,p)$ there is a finite positive random variable $d(\omega)$ and some $n_{0}<\infty$ such that for every $n\geq n_{0}$ one has

Thus for $t\in[2^{n},2^{n+1})$ and $n\geq n_{0}$ one has for some finite constant $K<\infty$

The latter display then guarantees that for $t\geq 2^{n_{0}}$ we have with probability one

By the bounds of Theorem A.1 we see that there are constants $0<K<\infty$ (that may change from line to line) and $0<q<\infty$ such that

The first inequality follows by Theorem A.1, the second inequality follows by Proposition 1 in and the third inequality follows by Condition 2.1.

The latter display implies that there is a finite random variable $\bar{J}_{\infty,0}^{(2)}$ such that

The last term that we need to consider is the martingale term

Notice that the Burkholder-Davis-Gundy inequality and the bounds of Theorem A.1 (doing calculations similar to the ones for the term $J_{t,0}^{(2)}$) give us that for some finite constant $K<\infty$, we have

Thus, by Doob’s martingale convergence theorem there is a square integrable random variable $\bar{J}_{\infty,0}^{(3)}$ such that

Let us now go back to (3). Using the terms $J_{t}^{(1)}$, $J_{t,0}^{(2)}$ and $J_{t,0}^{(3)}$ we can write

The last display together with (3.2), (3.3) and (3.4) imply the statement of the lemma. ∎

Assume that Conditions 2.1, 2.2 and 2.3 hold. Choose $\lambda>0$ such that for a given $\kappa>0$, one has $3\lambda+\frac{\lambda}{4\kappa}=\frac{1}{2L_{\nabla\bar{g}}}$, where $L_{\nabla\bar{g}}$ is the Lipschitz constant of $\nabla\bar{g}$. For $k$ large enough and for $\eta>0$ small enough (potentially random depending on $k$), one has $\int_{\tau_{k}}^{\sigma_{k,\eta}}\alpha_{s}ds>\lambda$. In addition we also have $\frac{\lambda}{2}\leq\int_{\tau_{k}}^{\sigma_{k}}\alpha_{s}ds\leq\lambda$ with probability one.

We proceed with an argument via contradiction. In particular let us assume that $\int_{\tau_{k}}^{\sigma_{k,\eta}}\alpha_{s}ds\leq\lambda$ and let us choose arbitrarily some $\epsilon>0$ such that $\epsilon\leq\lambda/8$.

Let us now make some remarks that are independent of the sign of $\int_{\tau_{k}}^{\sigma_{k,\eta}}\alpha_{s}ds-\lambda$. Due to the summability condition $\int_{0}^{\infty}\alpha^{2}_{t}dt<\infty$, $\frac{\kappa}{\|\nabla\bar{g}(\theta_{\tau_{k}})\|}\leq 1$ and Conditions 2.1 and 2.3, we have that

Hence, the martingale convergence theorem applies to the martingale $\int_{0}^{t}\alpha_{s}\frac{\kappa}{\|\nabla\bar{g}(\theta_{\tau_{k}})\|}\nabla_{\theta}f(X_{s},\theta_{s})\sigma^{-1}dW_{s}$. This means that there exists a square integrable random variable $M$ such that $\int_{0}^{t}\alpha_{s}\frac{\kappa}{\|\nabla\bar{g}(\theta_{\tau_{k}})\|}\nabla_{\theta}f(X_{s},\theta_{s})\sigma^{-1}dW_{s}\rightarrow M$ both almost surely and in $L^{2}$. This means that for the given $\epsilon>0$ there is $k$ large enough such that $\left\|\int_{\tau_{k}}^{\sigma_{k,\eta}}\alpha_{s}\frac{\kappa}{\|\nabla\bar{g}(\theta_{\tau_{k}})\|}\nabla_{\theta}f(X_{s},\theta_{s})\sigma^{-1}dW_{s}\right\|<\epsilon$ almost surely.

Let us also assume that for the given $k$, $\eta$ is so small such that for any $s\in[\tau_{k},\sigma_{k,\eta}]$ one has $\|\nabla\bar{g}(\theta_{s})\|\leq 3\|\nabla\bar{g}(\theta_{\tau_{k}})\|$.

Let us next bound appropriately the Euclidean norm of the vector-valued random variable

By Lemma 3.1 we have that for the same $0<\epsilon<\lambda/8$ that was chosen before there is $k$ large enough such that almost surely

Hence, using also the fact that $\frac{\kappa}{\|\nabla\bar{g}(\theta_{\tau_{k}})\|}\leq 1$ we obtain

The latter then implies that we should have

The latter statement will then imply that

But then we would necessarily have that $\int_{\tau_{k}}^{\sigma_{k,\eta}}\alpha_{s}ds>\lambda$, since otherwise $\sigma_{k,\eta}\in[\tau_{k},\sigma_{k}]$ which is impossible.

Next we move on to prove the second statement of the lemma. By definition we have $\int_{\tau_{k}}^{\sigma_{k}}\alpha_{s}ds\leq\lambda$. So it remains to show that $\frac{\lambda}{2}\leq\int_{\tau_{k}}^{\sigma_{k}}\alpha_{s}ds$. Since we know that $\int_{\tau_{k}}^{\sigma_{k,\eta}}\alpha_{s}ds>\lambda$ and because for $k$ large enough and $\eta$ small enough one should have $\int_{\sigma_{k}}^{\sigma_{k,\eta}}\alpha_{s}ds\leq\lambda/2$, we obtain that

Lemma 3.3 shows that the function $\bar{g}$ and its first two derivatives are uniformly bounded in $\theta$.

Assume Conditions 2.1, 2.2 and 2.3. For any $q>0$, there is a constant $K$ such that

In addition we also have that there is a constant $C<\infty$ such that $\sum_{i=0}^{2}\|\nabla^{i}_{\theta}\bar{g}(\theta)\|\leq C$.

By Theorem 1 in , the density $\mu$ of the measure $\pi$ admits, for any $p$, a constant $C_{p}$ such that $\mu(x)\leq\frac{C_{p}}{1+|x|^{p}}$. Choosing $p$ large enough that $\int_{\mathcal{X}}\frac{1+|x|^{q}}{1+|x|^{p}}dy<\infty$, we then obtain

concluding the proof of the first statement of the lemma. Let us now focus on the second part of the lemma. We only prove the claim for $i=0$, since due to the bounds in Condition 2.3, the proof for $i=1,2$ is the same. By Condition 2.3 and by the first part of the lemma, we have that there exist constants $0<q,K,C<\infty$ such that

Our next goal is to show that if the index $k$ is large enough, then $\bar{g}$ decreases, in the sense of Lemma 3.4.

Assume Conditions 2.1, 2.2 and 2.3. Suppose that there are an infinite number of intervals $I_{k}=[\tau_{k},\sigma_{k})$. There is a fixed constant $\gamma=\gamma(\kappa)>0$ such that for $k$ large enough, one has

Let’s first consider $\Theta_{1,k}$. Notice that for all $s\in[\tau_{k},\sigma_{k}]$ one has $\frac{\|\nabla\bar{g}(\theta_{\tau_{k}})\|}{2}\leq\|\nabla\bar{g}(\theta_{s})\|\leq 2\|\nabla\bar{g}(\theta_{\tau_{k}})\|$. Hence, for sufficiently large $k$, we have the upper bound:

since Lemma 3.1 proved that $\int_{\tau_{k}}^{\sigma_{k}}\alpha_{s}ds\geq\frac{\lambda}{2}$ for sufficiently large $k$.

We next address $\Theta_{2,k}$ and show that it becomes small as $k\rightarrow\infty$. First notice that we can trivially write

By Condition 2.3 and Itô isometry we have

where $R_{s}$ is defined via (3.5). Hence, by Doob’s martingale convergence theorem there is a square integrable random variable $M$ such that $\int_{0}^{t}\alpha_{s}\left<\frac{\nabla\bar{g}(\theta_{s})}{R_{s}},\nabla_{\theta}f(X_{s},\theta_{s})dW_{s}\right>\rightarrow M$ both almost surely and in $L^{2}$. The latter statement implies that for a given $\epsilon>0$ there is $k$ large enough such that almost surely

where we have used Condition 2.3 and Lemma 3.3. Bound (3.7) implies that

is finite almost surely, which in turn implies that there is a finite random variable $\Theta_{3}^{\infty}$ such that

with probability one. Since $\Theta_{3}^{\infty}$ is finite, $\int_{\tau_{k}}^{\sigma_{k}}\frac{\alpha_{s}^{2}}{2}\text{tr}\left[(\nabla_{\theta}f(X_{s},\theta_{s})\sigma^{-1})(\nabla_{\theta}f(X_{s},\theta_{s})\sigma^{-1})^{\top}\nabla_{\theta}\nabla_{\theta}\bar{g}(\theta_{s})\right]ds\rightarrow 0$ as $k\rightarrow\infty$ with probability one.

Finally, we address $\Theta_{4,k}$. Let us consider the function $G(x,\theta)=\left<\nabla_{\theta}\bar{g}(\theta),\nabla_{\theta}g(x,\theta)-\nabla_{\theta}\bar{g}(\theta)\right>$. The function $G(x,\theta)$ satisfies the centering condition (A.1) of Theorem A.1. Therefore, the Poisson equation (A.2) with right hand side $G(x,\theta)$ will have a unique smooth solution, say $v(x,\theta)$, that grows at most polynomially in $x$. Let us apply Itô formula to the function $u(t,x,\theta)=\alpha_{t}v(x,\theta)$ that is solution to this Poisson equation.

Rearranging the previous Itô formula yields

Following the exact same steps as in the proof of Lemma 3.1 gives us that $\lim_{k\rightarrow\infty}\left\lVert\Theta_{4,k}\right\rVert\rightarrow 0$ almost surely.

We now return to $\bar{g}(\theta_{\sigma_{k}})-\bar{g}(\theta_{\tau_{k}})$ and provide an upper bound which is negative. For sufficiently large $k$, we have that:

Choose $\epsilon=\min\{\frac{\lambda\kappa^{2}}{32},\frac{\lambda}{32}\}$. On the one hand, if $\|\nabla\bar{g}(\theta_{\tau_{k}})\|\geq 1$:

On the other hand, if $\|\nabla\bar{g}(\theta_{\tau_{k}})\|\leq 1$, then

Finally, let $\gamma=\frac{\kappa^{2}}{32}\lambda$ and the proof of the lemma is complete.

Assume Conditions 2.1, 2.2 and 2.3. Suppose that there are an infinite number of intervals $I_{k}=[\tau_{k},\sigma_{k})$. There is a fixed constant $\gamma_{1}<\gamma$ such that for $k$ large enough,

First, recall that $\left\lVert\nabla\bar{g}(\theta_{t})\right\rVert\leq\kappa$ for $t\in J_{k}=[\sigma_{k-1},\tau_{k}]$. Similar to before, we have that:

The right hand side (RHS) of equation (3.8) converges almost surely to as $k\rightarrow\infty$ as a consequence of similar arguments as given in Lemma 3.4. Indeed, the treatment of the second and third terms on the RHS of (3.8) are exactly the same as in Lemma 3.4. It remains to show that the first term on the RHS of (3.8) converges almost surely to as $k\rightarrow\infty$.

As shown in Lemma 3.4, $\int_{\sigma_{k-1}}^{\tau_{k}}\alpha_{s}\left<\frac{\nabla\bar{g}(\theta_{s})}{R_{s}},\nabla_{\theta}f(X_{s},\theta_{s})\sigma^{-1}dW_{s}\right>\rightarrow 0$ as $k\rightarrow\infty$ almost surely. Finally, note that $\left\lVert\nabla\bar{g}(\theta_{\sigma_{k-1}})\right\rVert\leq\kappa$ (except when $\sigma_{k-1}=\tau_{k}$, in which case the interval $J_{k}$ is length 0 and hence the integral (3.9) over $J_{k}$ is ). Then, $\int_{\sigma_{k-1}}^{\tau_{k}}\alpha_{s}\left<\nabla\bar{g}(\theta_{s}),\nabla_{\theta}f(X_{s},\theta_{s})\sigma^{-1}dW_{s}\right>\rightarrow 0$ as $k\rightarrow\infty$ almost surely.

Therefore, with probability one, $\bar{g}(\theta_{\tau_{k}})-\bar{g}(\theta_{\sigma_{k-1}})\leq\gamma_{1}<\gamma$ for sufficiently large $k$.

Choose a $\kappa>0$. First, consider the case where there are a finite number of times $\tau_{k}$. Then, there is a finite $T$ such that $\left\lVert\nabla\bar{g}(\theta_{t})\right\rVert<\kappa$ for $t\geq T$. Now, consider the other case where there are an infinite number of times $\tau_{k}$ and use Lemmas 3.4 and 3.5. With probability one,

for sufficiently large $k$. Choose a $K$ such that (3.10) holds for $k\geq K$. This leads to:

Let $n\rightarrow\infty$ and then $\bar{g}(\theta_{\tau_{n+1}})\rightarrow-\infty$. However, we also have that by definition $\bar{g}(\theta)\geq 0$. This is a contradiction, and therefore almost surely there are a finite number of times $\tau_{k}$.

Consequently, there exists a finite time $T$ (possibly random) such that almost surely $\left\lVert\nabla\bar{g}(\theta_{t})\right\rVert<\kappa$ for $t\geq T$. Since the original $\kappa>0$ was arbitrarily chosen, this shows that $\left\lVert\nabla\bar{g}(\theta_{t})\right\rVert\rightarrow 0$ as $t\rightarrow\infty$ almost surely. ∎

Estimating the Coefficient Function of the Diffusion Term and Generalizations

The stochastic gradient descent update in continuous time follows the stochastic differential equations:

We assume that $\sigma^{\ast}(x)$ is such that the process $X_{t}$ is ergodic with a unique invariant measure (for example one may assume that it is non-degenerate, i.e., bounded away from zero and bounded by above). In addition, we assume that $w(x,\nu)$ satisfies the same assumptions as $g(x,\theta)$ does in Condition 2.3.

From the previous results in Section 3, $\lim_{t\rightarrow\infty}\|\nabla\bar{g}(\theta_{t})\|=0$ as $t\rightarrow\infty$ with probability 1. Let’s study the convergence of the stochastic gradient descent algorithm (4) for $\nu_{t}$. By Itô’s formula,

Applying exactly the same procedure as in Section 3, $\lim_{t\rightarrow\infty}\|\nabla\bar{w}(\nu_{t})\|=0$ as $t\rightarrow\infty$ with probability 1. We omit the details as the proof is exactly the same as in Section 3.

Notice also that $\sigma^{\ast}(x)$ is not identifiable; for example, $X_{t}$ has the same distribution under the diffusion coefficient $-\sigma^{\ast}(x)$. Only $\sigma^{\ast}(x)\sigma^{\ast,\top}(x)$ is identifiable. We are therefore essentially estimating a model $\sigma(x,\nu)\sigma^{\top}(x,\nu)$ for $\sigma^{\ast}(x)\sigma^{\ast,\top}(x)$.

We close this section with the following remark.

Model Estimation: Numerical Analysis

We implement SGDCT for several applications and numerically analyze the convergence. Section 5.1 studies continuous-time stochastic gradient descent for the Ornstein-Uhlenbeck process, which is widely used in finance, physics, and biology. Section 5.2 studies the multidimensional Ornstein-Uhlenbeck process. Section 5.3 estimates the diffusion coefficient in Burger’s equation with continuous-time stochastic gradient descent. Burger’s equation is a widely-used nonlinear partial differential equation which is important to fluid mechanics, acoustics, and aerodynamics. Burger’s equation is extensively used in engineering. In Section 5.4, we show how SGDCT can be used for reinforcement learning. In the final example, the drift and volatility functions for the multidimensional CIR process are estimated. The CIR process is widely used in financial modeling.

For the numerical experiments, we use an Euler scheme with a time step of $10^{-2}$. The learning rate is $\alpha_{t}=\min(\alpha,\alpha/t)$ with $\alpha=10^{-2}$. We simulate data from (5.1) for a particular $\theta^{\ast}$ and the stochastic gradient descent attempts to learn a parameter $\theta_{t}$ which fits the data well. $\theta_{t}$ is the statistical estimate for $\theta^{\ast}$ at time $t$. If the estimation is accurate, $\theta_{t}$ should of course be close to $\theta^{\ast}$. This example can be placed in the form of the original class of equations (1.1) by setting $f(x,\theta)=c(m-x)$ and $f^{\ast}(x)=f(x,\theta^{\ast})$.

We study $10,500$ cases. For each case, a different $\theta^{\ast}$ is generated uniformly at random in the range $\times$. For each case, we solve for the parameter $\theta_{t}$ over the time period $[0,T]$ for $T=10^{6}$. To summarize:

Generate a random $\theta^{\ast}$ in $\times$

Simulate a single path of $X_{t}$ given $\theta^{\ast}$ and simultaneously solve for the path of $\theta_{t}$ on $[0,T]$

The accuracy of $\theta_{t}$ at times $t=10^{2},10^{3},10^{4},10^{5}$, and $10^{6}$ is reported in Table 1. Figures 1 and 2 plot the mean error in percent and mean squared error (MSE) against time. In the table and figures, the “error” is $|\theta_{t}^{n}-\theta^{\ast,n}|$ where $n$ represents the $n$-th case. The “error in percent” is $100\times\frac{|\theta_{t}^{n}-\theta^{\ast,n}|}{|\theta^{\ast,n}|}$. The “mean error in percent” is the average of these errors, i.e. $\frac{100}{N}\sum_{n=1}^{N}\frac{|\theta_{t}^{n}-\theta^{\ast,n}|}{|\theta^{\ast,n}|}$.

Finally, we also track the objective function $\bar{g}(\theta_{t})$ over time. Figure 3 plots the error $\bar{g}(\theta_{t})$ against time. Since the limiting distribution $\pi(x)$ of (5.2) is Gaussian with mean $m^{\ast}$ and variance $\frac{1}{2c^{\ast}}$, we have that:

2 Multidimensional Ornstein-Uhlenbeck process

For the numerical experiments, we use an Euler scheme with a time step of $10^{-2}$. The learning rate is $\alpha_{t}=\min(\alpha,\alpha/t)$ with $\alpha=10^{-1}$. We simulate data from (5.2) for a particular $\theta^{\ast}=(M^{\ast},A^{\ast})$ and the stochastic gradient descent attempts to learn a parameter $\theta_{t}$ which fits the data well. $\theta_{t}$ is the statistical estimate for $\theta^{\ast}$ at time $t$. If the estimation is accurate, $\theta_{t}$ should of course be close to $\theta^{\ast}$. This example can be placed in the form of the original class of equations (1.1) by setting $f(x,\theta)=M-Ax$ and $f^{\ast}(x)=f(x,\theta^{\ast})$.

The matrix $A^{\ast}$ must be generated carefully to ensure that $X_{t}$ is ergodic and has a stable equilibrium point. If some of $A^{\ast}$’s eigenvalues have negative real parts, then $X_{t}$ can become unstable and grow arbitrarily large. Therefore, we randomly generate matrices $A^{\ast}$ which are strictly diagonally dominant. $A^{\ast}$’s eigenvalues are therefore guaranteed to have positive real parts and $X_{t}$ will be ergodic. To generate random strictly diagonally dominant matrices $A^{\ast}$, we first generate $A_{i,j}^{\ast}$ uniformly at random in the range $$for$i\neq j$. Then, we set$A_{i,i}^{\ast}=\sum_{j\neq i}A_{i,j}^{\ast}+U_{i,i}$where$U_{i,i}$is generated randomly in$$.$M_{i}^{\ast}$for$i=1,\ldots,d$is also generated randomly in$$.

We study $525$ cases and analyze the error in Table 2. Figures 4 and 5 plot the error over time.

3 Burger’s Equation

The stochastic Burger’s equation that we consider is given by:

where $x\in$ and $W(t,x)$ is a Brownian sheet. The finite-difference discretization of (5.3) satisfies a system of nonlinear stochastic differential equations (for instance, see or ). We use continuous-time stochastic gradient descent to learn the diffusion parameter $\theta$.

We use the following finite difference scheme for Burger’s equation:

For our numerical experiment, the boundary conditions $u(t,x=0)=0$ and $u(t,x=1)=1$ are used and $\sigma=0.1$. (5.4) is simulated with the Euler scheme (i.e., we solve Burger’s equation with explicit finite difference). A spatial discretization of $\Delta x=.01$ and a time step of $10^{-5}$ are used. The learning rate is $\alpha_{t}=\min(\alpha,\alpha/t)$ with $\alpha=10^{-3}$. The small time step is needed to avoid instability in the explicit finite difference scheme. We simulate data from (5.3) for a particular diffusion coefficient $\theta^{\ast}$ and the stochastic gradient descent attempts to learn a diffusion parameter $\theta_{t}$ which fits the data well. $\theta_{t}$ is the statistical estimate for $\theta^{\ast}$ at time $t$. If the estimation is accurate, $\theta_{t}$ should of course be close to $\theta^{\ast}$.

This example can be placed in the form of the original class of equations (1.1). Let $f_{i}$ be the $i$-th element of the function $f$. Then, $f_{i}(u,\theta)=\theta\frac{u(t,x_{i+1})-2u(t,x_{i})+u(t,x_{i-1})}{\Delta x^{2}}-u(t,x_{i})\frac{u(t,x_{i+1})-u(t,x_{i-1})}{2\Delta x}$. Similarly, let $f_{i}^{\ast}$ be the $i$-th element of the function $f^{\ast}$. Then, $f_{i}^{\ast}(u)=f_{i}(u,\theta^{\ast})$.

We study $525$ cases. For each case, a different $\theta^{\ast}$ is generated uniformly at random in the range $[.1,10]$. This represents a wide range of physical cases of interest, with $\theta^{\ast}$ ranging over two orders of magnitude. For each case, we solve for the parameter $\theta_{t}$ over the time period $[0,T]$ for $T=100$.

The accuracy of $\theta_{t}$ at times $t=10^{-1},10^{0},10^{1},$ and $10^{2}$ is reported in Table 3. Figures 6 and 7 plot the mean error in percent and mean squared error against time. The convergence of $\theta_{t}$ to $\theta^{\ast}$ is fairly rapid in time.

4 Reinforcement Learning

We consider the classic reinforcement learning problem of balancing a pole on a moving cart (see ). The goal is to balance a pole on a cart and to keep the cart from moving outside the boundaries via applying a force of $\pm 10$ Newtons.

The position $x$ of the cart, the velocity $\dot{x}$ of the cart, angle of the pole $\beta$, and angular velocity $\dot{\beta}$ of the pole are observed. The dynamics of $s=(x,\dot{x},\beta,\dot{\beta})$ satisfy a set of ODEs (see ):

where $g$ is the acceleration due to gravity, $m_{c}$ is the mass of the cart, $m$ is the mass of the pole, $2l$ is the length of the pole, $\mu_{c}$ is the coefficient of friction of the cart on the ground, $\mu_{p}$ is the coefficient of friction of the pole on the cart, and $F_{t}\in\{-10,10\}$ is the force applied to the cart.

For this example, $f^{\ast}(s)=(\dot{x},\ddot{x},\dot{\beta},\ddot{\beta})$. The model $f(s,\theta)=(f_{1}(s,\theta),f_{2}(s,\theta),f_{3}(s,\theta),f_{4}(s,\theta))$ where $f_{i}(s,\theta)$ is a single-layer neural network with rectified linear units.

The boundary is $x=\pm 2.4$ meters and the pole must not be allowed to fall below $\beta=\frac{24}{360\pi}$ radians (the frame of reference is chosen such that the perfectly upright is radians). A reward of $+1$ is received every $0.02$ seconds if $\left\lVert x\right\rVert\leq 2.4$ and $\left\lVert\theta\right\rVert\leq\frac{24}{360\pi}$. A reward of $-100$ is received (and the episode ends) if the cart moves beyond $x=\pm 2.4$ or the pole falls below $\beta=\frac{24}{360\pi}$ radians. The sum of these rewards across the entire episode is the reward for that episode. The initial state $(x,\dot{x},\beta,\dot{\beta})$ at the start of an episode is generated uniformly at random in $[-.05,.05]^{4}$. For our numerical experiment, we assume that the rule for receiving the rewards and the distribution of the initial state are both known. An action of $\pm 10$ Newtons may be chosen every $0.02$ seconds. This force is then applied for the duration of the next $0.02$ seconds. The system (5.5) is simulated using an Euler scheme with a time step size of $10^{-3}$ seconds.

The goal, of course, is to statistically learn the optimal actions in order to achieve the highest possible reward. This requires both: 1) statistically learning the physical dynamics of $(x,\dot{x},\beta,\dot{\beta})$ and 2) finding the optimal actions given these dynamics in order to achieve the highest possible reward. The dynamics $(x,\dot{x},\beta,\dot{\beta})$ satisfy the set of ODEs (5.5); these dynamics can be learned using continuous-time stochastic gradient descent. We use a neural network for $f$. Given the estimated dynamics $f$, we use a policy gradient method to estimate the optimal actions. The approach is summarized below.

For time $[0,T_{\textrm{end of episode}}]$:

Update the model $f(s,\theta)$ for the dynamics using continuous-time stochastic gradient descent.

Periodically update the optimal policy $\mu(s,a,\theta^{\mu})$ using policy gradient method. The optimal policy is learned using data simulated from the model $f(s,\theta)$. Actions are randomly selected via the policy $\mu$.

The policy $\mu$ is a neural network with parameters $\theta^{\mu}$. We use a single hidden layer with rectified linear units followed by a softmax layer for $\mu(s,a,\theta^{\mu})$ and train it using policy gradients.Let $r_{e,t}$ be the reward for episode $e$ at time $t$. Let $R_{t,e}=\sum_{t^{\prime}=t+1}^{T_{\textrm{end of episode}}}\gamma^{t^{\prime}-t}r_{e,t^{\prime}}$ be the cumulative discounted reward from episode $e$ after time $t$ where $\gamma\in$ is the discount factor. Stochastic gradient descent is used to learn the parameter $\theta^{\mu}$: $\theta^{\mu}\leftarrow\theta^{\mu}+\eta_{e}R_{t,e}\frac{\partial}{\partial\theta^{\mu}}\log\mu(s_{t},a_{t},\theta^{\mu})$ where $\eta_{e}$ is the learning rate. In practice, the cumulative discounted rewards are often normalized across an episode. The policy $\mu(s,a,\theta^{\mu})$ gives the probability of taking action $a$ conditional on being in the state $s$.

$525$ cases are run, each for $25$ hours. The optimal policy is learned using the estimated dynamics $f(s,\theta)$ and is updated every $5$ episodes. Table 4 reports the results at fixed episodes using continuous-time stochastic gradient descent. Table 5 reports statistics on the number of episodes required until a target episodic reward ($100$, $500$, $1000$) is first achieved.

Alternatively, one could directly apply policy gradient to learn the optimal action using the observed data. This approach does not use continuous-time stochastic gradient descent to learn the model dynamics, but instead directly learns the optimal policy from the data. Again using 525 cases, we report the results in Table 6 for directly learning the optimal policy without using continuous-time stochastic gradient descent to learn the model dynamics. Comparing Tables 4 and 6, it is clear that using continuous-time stochastic gradient descent to learn the model dynamics allows for the optimal policy to be learned significantly more quickly. The rewards are much higher when using continuous-time stochastic gradient descent (see Table 4) than when not using it (see Table 6).

5 Estimating both the drift and volatility functions for the multidimensional CIR process

American Options

High-dimensional American options are extremely computationally challenging to solve with traditional numerical methods such as finite difference. Here we propose a new approach using statistical learning to solve high-dimensional American options. SGDCT achieves a high accuracy on two benchmark problems with 100 dimensions.

Before describing the SGDCT algorithm for American options, it is important to note that traditional stochastic gradient descent faces certain difficulties in this class of problems. Some brief remarks are provided below regarding this fact; the authors plan to elaborate on these issues in more detail in a future work. The well-known Q-learning algorithm uses stochastic gradient descent to minimize an approximation to the discrete-time Hamilton-Jacobi-Bellman equation. To demonstrate the challenges and the issues that arise, consider using Q-learning to estimate the value function:

where $\gamma>0$ is a discount factor and $r(x)$ is a reward function. The function $Q(x,\theta)$ is an approximation for the value function $V(x)$. The parameter $\theta$ must be estimated. The traditional approach would discretize the dynamics (6.1) and then apply a stochastic gradient descent update to the objective function:

This results in the stochastic gradient descent algorithm:

Although now computationally efficient, the Q-learning algorithm (6.4) is now biased (due to ignoring the inner expectations). Furthermore, when $\Delta\rightarrow 0$, the Q-learning algorithm (6.4) blows up. A quick investigation shows that the term $\frac{1}{\Delta}(W_{t+\Delta}-W_{t})^{2}=O(1)$ arises while all other terms are $O(\Delta)$ or $O(\sqrt{\Delta})$.

The SGDCT algorithm is unbiased and computationally efficient. It can be directly derived by letting $\Delta\rightarrow 0$ and using Itô’s formula in (6.3):

Note that computationally challenging terms in (6.3) become differential operators in (6.7), which are usually easier to evaluate. This is one of the advantages of developing the theory in continuous time for continuous-time models. Once the continuous-time algorithm is derived, it can be appropriately discretized for numerical solution.

2 SGDCT for American Options

$\mathcal{L}_{x}$ is the infinitesimal generator for the $X$ process. The continuous-time algorithm (6.7) is run for many iterations $n=0,1,2,\ldots$ until convergence. See the authors’ paper for implementation details on pricing American options with deep learning.

We implement the SGDCT algorithm (6.7) using a deep neural network for the function $Q(t,x;\theta)$. Two benchmark problems are considered where semi-analytic solutions are available. The SGDCT algorithm’s accuracy is evaluated for American options in $d=100$ dimensions, and the results are presented in Table 8.

Appendix A On a related Poisson equation

We recall the following regularity result from on the Poisson equations in the whole space, appropriately stated to cover our case of interest.

and that for some positive constants $K$ and $q$,

Let $\mathcal{L}_{x}$ be the infinitesimal generator for the $X$ process. Then the Poisson equation

References