Ornstein–Uhlenbeck process

A temporally homogeneous Ornstein–Uhlenbeck process can be represented as a scaled, time-transformed Wiener process:

${\displaystyle x_{t}={\frac {\sigma }{\sqrt {2\theta }}}e^{-\theta t}W_{e^{2\theta t}}}$

where ${\displaystyle W_{t}}$ is the standard Wiener process. This is roughly Theorem 1.2 in Doob 1942. Equivalently, with the change of variable ${\displaystyle s=e^{2\theta t}}$ this becomes

${\displaystyle W_{s}={\frac {\sqrt {2\theta }}{\sigma }}s^{1/2}x_{(\ln s)/(2\theta )},\qquad s>0}$

Using this mapping, one can translate known properties of ${\displaystyle W_{t}}$ into corresponding statements for ${\displaystyle x_{t}}$. For instance, the law of the iterated logarithm for ${\displaystyle W_{t}}$ becomes[1]

${\displaystyle \limsup _{t\to \infty }{\frac {x_{t}}{\sqrt {(\sigma ^{2}/\theta )\ln t}}}=1,\quad {\text{with probability 1.}}}$

The stochastic differential equation for ${\displaystyle x_{t}}$ can be formally solved by variation of parameters.[7] Writing

${\displaystyle f(x_{t},t)=x_{t}e^{\theta t}\,}$

we get

${\displaystyle {\begin{aligned}df(x_{t},t)&=\theta \,x_{t}\,e^{\theta t}\,dt+e^{\theta t}\,dx_{t}\\[6pt]&=e^{\theta t}\theta \,\mu \,dt+\sigma \,e^{\theta t}\,dW_{t}.\end{aligned}}}$

Integrating from ${\displaystyle 0}$ to ${\displaystyle t}$ we get

${\displaystyle x_{t}e^{\theta t}=x_{0}+\int _{0}^{t}e^{\theta s}\theta \,\mu \,ds+\int _{0}^{t}\sigma \,e^{\theta s}\,dW_{s}\,}$

whereupon we see

${\displaystyle x_{t}=x_{0}\,e^{-\theta t}+\mu \,(1-e^{-\theta t})+\sigma \int _{0}^{t}e^{-\theta (t-s)}\,dW_{s}.\,}$

From this representation, the first moment (i.e. the mean) is shown to be

${\displaystyle \operatorname {E} (x_{t})=x_{0}e^{-\theta t}+\mu (1-e^{-\theta t})\!\ }$

assuming ${\displaystyle x_{0}}$ is constant. Moreover, the Itō isometry can be used to calculate the covariance function by

${\displaystyle {\begin{aligned}\operatorname {cov} (x_{s},x_{t})&=\operatorname {E} [(x_{s}-\operatorname {E} [x_{s}])(x_{t}-\operatorname {E} [x_{t}])]\\[5pt]&=\operatorname {E} \left[\int _{0}^{s}\sigma e^{\theta (u-s)}\,dW_{u}\int _{0}^{t}\sigma e^{\theta (v-t)}\,dW_{v}\right]\\[5pt]&=\sigma ^{2}e^{-\theta (s+t)}\operatorname {E} \left[\int _{0}^{s}e^{\theta u}\,dW_{u}\int _{0}^{t}e^{\theta v}\,dW_{v}\right]\\[5pt]&={\frac {\sigma ^{2}}{2\theta }}\,e^{-\theta (s+t)}(e^{2\theta \min(s,t)}-1)\\[5pt]&={\frac {\sigma ^{2}}{2\theta }}\left(e^{-\theta |t-s|}-e^{-\theta (t+s)}\right).\end{aligned}}}$

Since the Itô integral of deterministic integrand is normally distributed, it follows that

${\displaystyle x_{t}=x_{0}e^{-\theta t}+\mu (1-e^{-\theta t})+{\tfrac {\sigma }{\sqrt {2\theta }}}W_{1-e^{-2\theta t}}}$

The infinitesimal generator of the process is[8]

$${\displaystyle Lf=-\theta (x-\mu )f'+{\frac {1}{2}}\sigma ^{2}f''}$$

If we let ${\displaystyle y=x{\sqrt {\frac {2\mu }{\sigma ^{2}}}}}$, then the eigenvalue equation simplifies to:

$${\displaystyle {\frac {d^{2}}{dy^{2}}}\phi -y{\frac {d}{dy}}\phi -{\frac {\lambda }{\mu }}\phi =0}$$

which is the defining equation for Hermite polynomials. Its solutions are ${\displaystyle \phi (y)=He_{n}(y)}$, with ${\displaystyle \lambda =n\mu }$, which implies that the mean first passage time for a particle to hit a point on the boundary is on the order of ${\displaystyle \mu ^{-1}}$.

The Ornstein–Uhlenbeck process is a prototype of a noisy relaxation process. A canonical example is a Hookean spring (harmonic oscillator) with spring constant ${\displaystyle k}$ whose dynamics is overdamped with friction coefficient ${\displaystyle \gamma }$. In the presence of thermal fluctuations with temperature ${\displaystyle T}$, the length ${\displaystyle x(t)}$ of the spring fluctuates around the spring rest length ${\displaystyle x_{0}}$; its stochastic dynamics is described by an Ornstein–Uhlenbeck process with

${\displaystyle {\begin{aligned}\theta &=k/\gamma ,\\\mu &=x_{0},\\\sigma &={\sqrt {2k_{B}T/\gamma }},\end{aligned}}}$

where ${\displaystyle \sigma }$ is derived from the Stokes–Einstein equation ${\displaystyle D=\sigma ^{2}/2=k_{B}T/\gamma }$ for the effective diffusion constant.[12][13] This model has been used to characterize the motion of a Brownian particle in an optical trap.[13][14]

At equilibrium, the spring stores an average energy ${\displaystyle \langle E\rangle =k\langle (x-x_{0})^{2}\rangle /2=k_{B}T/2}$ in accordance with the equipartition theorem.[15]

The Ornstein–Uhlenbeck process is used in the Vasicek model of the interest rate.[16] The Ornstein–Uhlenbeck process is one of several approaches used to model (with modifications) interest rates, currency exchange rates, and commodity prices stochastically. The parameter ${\displaystyle \mu }$ represents the equilibrium or mean value supported by fundamentals; ${\displaystyle \sigma }$ the degree of volatility around it caused by shocks, and ${\displaystyle \theta }$ the rate by which these shocks dissipate and the variable reverts towards the mean. One application of the process is a trading strategy known as pairs trade.[17][18][19]

A further implementation of the Ornstein–Uhlenbeck process is derived by Marcello Minenna in order to model the stock return under a lognormal distribution dynamics. This modeling aims at the determination of confidence interval to predict market abuse phenomena. [20] [21]

The Ornstein–Uhlenbeck process has been proposed as an improvement over a Brownian motion model for modeling the change in organismal phenotypes over time.[22] A Brownian motion model implies that the phenotype can move without limit, whereas for most phenotypes natural selection imposes a cost for moving too far in either direction. A meta-analysis of 250 fossil phenotype time-series showed that an Ornstein–Uhlenbeck model was the best fit for 115 (46%) of the examined time series, supporting stasis as a common evolutionary pattern.[23] This said, there are certain challenges to its use: model selection mechanisms are often biased towards preferring an OU process without sufficient support, and misinterpretation is easy to the unsuspecting data scientist.[24]

A multi-dimensional version of the Ornstein–Uhlenbeck process, denoted by the N-dimensional vector ${\displaystyle \mathbf {x} _{t}}$, can be defined from

${\displaystyle d\mathbf {x} _{t}=-{\boldsymbol {\beta }}\,\mathbf {x} _{t}\,dt+{\boldsymbol {\sigma }}\,d\mathbf {W} _{t}.}$

where ${\displaystyle \mathbf {W} _{t}}$ is an N-dimensional Wiener process, and ${\displaystyle {\boldsymbol {\beta }}}$ and ${\displaystyle {\boldsymbol {\sigma }}}$ are constant N×N matrices.[28] The solution is

${\displaystyle \mathbf {x} _{t}=e^{-{\boldsymbol {\beta }}t}\mathbf {x} _{0}+\int _{0}^{t}e^{-{\boldsymbol {\beta }}(t-t')}{\boldsymbol {\sigma }}\,d\mathbf {W} _{t'}}$

and the mean is

${\displaystyle \operatorname {E} (\mathbf {x} _{t})=e^{-{\boldsymbol {\beta }}t}\operatorname {E} (\mathbf {x} _{0}).}$

These expressions make use of the matrix exponential.

The process can also be described in terms of the probability density function ${\displaystyle P(\mathbf {x} ,t)}$, which satisfies the Fokker–Planck equation[29]

${\displaystyle {\frac {\partial P}{\partial t}}=\sum _{i,j}\beta _{ij}{\frac {\partial }{\partial x_{i}}}(x_{j}P)+\sum _{i,j}D_{ij}{\frac {\partial ^{2}P}{\partial x_{i}\,\partial x_{j}}}.}$

where the matrix ${\displaystyle {\boldsymbol {D}}}$ with components ${\displaystyle D_{ij}}$ is defined by ${\displaystyle {\boldsymbol {D}}={\boldsymbol {\sigma }}{\boldsymbol {\sigma }}^{T}/2}$. As for the 1d case, the process is a linear transformation of Gaussian random variables, and therefore itself must be Gaussian. Because of this, the transition probability ${\displaystyle P(\mathbf {x} ,t\mid \mathbf {x} ',t')}$ is a Gaussian which can be written down explicitly. If the real parts of the eigenvalues of ${\displaystyle {\boldsymbol {\beta }}}$ are larger than zero, a stationary solution ${\displaystyle P_{\text{st}}(\mathbf {x} )}$ moreover exists, given by

${\displaystyle P_{\text{st}}(\mathbf {x} )=(2\pi )^{-N/2}(\det {\boldsymbol {\omega }})^{-1/2}\exp \left(-{\frac {1}{2}}\mathbf {x} ^{T}{\boldsymbol {\omega }}^{-1}\mathbf {x} \right)}$

where the matrix ${\displaystyle {\boldsymbol {\omega }}}$ is determined from the Lyapunov equation ${\displaystyle {\boldsymbol {\beta }}{\boldsymbol {\omega }}+{\boldsymbol {\omega }}{\boldsymbol {\beta }}^{T}=2{\boldsymbol {D}}}$.[5]