# A discrete time Ornstein-Uhlenbeck type process

Let's consider a stochastic process $$X_{t+1}=X_{t}+\Delta X_t$$ where $$\Delta X_t$$ is a Gaussian with mean $$-\lambda X_t \Delta t$$ and variance $$2\Delta t$$ where $$\lambda$$ and $$\Delta t$$ are parameters of the model.

We can then write the transition probability density as $$$$P(X_{t+1}=x|X_{t}=y)=\frac{1}{\sqrt{4\pi\Delta t}}\exp\left[\frac{[x-y(1-\lambda \Delta t)]^2}{4\Delta t}\right]$$$$

We can find the stationary distribution of this process which will take the form $$$$P_{\infty}(x)=\sqrt{\frac{a}{2\pi}}e^{-ax^2/2}$$$$ where $$a$$ depends on the parameters $$\lambda$$ and $$\Delta t$$. It is easy to see that if we scale time by $$\Delta t$$ and in the limit $$\Delta t$$ goes to $$0$$, we get something like the Ornstein-Uhlenbeck process.

If we initialise the process in state $$X_0$$ distributed as the stationary state, and then want the $$t$$-step propagator of this process, which is, $$P(X_t=x|X_0=0)$$, it is easy to see that this will also be a Gaussian (as it is a sum of Gaussians). But is there a way to compute the distribution exactly?

So you could start with the propagator, chain them together and integrate over the intermediate values and you'll end with a discretized version of the path integral. This is, however, rather tedious, so let's just instead write the stochastic difference equation: $$\Delta X_n = -\lambda \Delta t + \sqrt{2} \Delta W_n$$ Here $$\Delta W_n$$ are independent Gaussian increments with variance $$\Delta t$$ and mean $$0$$; $$W_0 = 0$$. Note then that $$\sum_{n=0}^{N-1} \Delta W_n = W_N - W_0 = W_N \sim \mathcal{N}(0, \Delta t N)$$ as variance is additive. So I guess you can see where we're going with this: $$X_N - X_0 = \sum_{n=0}^{N-1} \Delta X_n = \sum_{n=0}^{N-1} (-\lambda \Delta t + \sqrt{2} \Delta W_n) = -\lambda \Delta t N + \sqrt{2} W_N$$ So $$X_N$$ is normally distributed with mean $$X_0 -\lambda \Delta t N$$ and variance $$2\Delta t N$$.

Now none of this is related to Orstein--Uhlenbeck. Maybe instead you meant to say: $$\Delta X_n = -\lambda X_n \Delta t + \sqrt{2} \Delta W_n$$ This, like the above, can be solved much like the corresponding SDE. Start by the simpler: $$\Delta Z_n = -\lambda Z_n \Delta t \Rightarrow Z_{n+1} = (1-\lambda \Delta t) Z_n \Rightarrow Z_n = (1-\lambda \Delta t)^n Z_0$$, and so we make the substitution $$X_n = (1-\lambda \Delta t)^n Y_n$$ to the original equation and are left with: $$(1-\lambda \Delta t)^{n+1}Y_{n+1} - (1-\lambda \Delta t)^nY_n = -\lambda(1-\lambda \Delta t)^n Y_n \Delta t + \sqrt{2}\Delta W_n$$ or rearranging $$\Delta Y_n = \frac{\sqrt{2}}{(1-\lambda\Delta t)^{n+1}}\Delta W_n$$ Thus $$Y_N$$ is again just a sum of normals, so mean $$Y_0$$ and variance $$2\sum_{n=0}^{N-1}\left(\frac{1}{(1-\lambda\Delta t)^2}\right)^{n+1}\Delta t = 2\Delta t \frac{1-(1-\lambda\Delta t)^{-2N}}{(1-\lambda\Delta t)^2 - 1}$$, and so $$X_N \sim \mathcal{N}\left(X_0(1-\lambda \Delta t)^N,2\Delta t \frac{(1-\lambda\Delta t)^{2N} - 1}{(1-\lambda\Delta t)^{2} - 1}\right)$$

In closing, some python code:

import numpy as np
import matplotlib.pyplot as plt

dt = .01
N  = 100
l  = .05
X0 = .75

M = 1000000
res = np.zeros(M) + X0
for _ in range(N):
res += -l*res*dt + np.sqrt(2*dt) * np.random.randn(M)

yy, xx = np.histogram(res, 128, normed=True)
plt.plot(.5*(xx[1:]+xx[:-1]), yy, 'b')
mu    = X0 * (1-l*dt)**N
sigma = np.sqrt(2*dt*(((1-l*dt)**(2*N)-1) / ((1-l*dt)**2 - 1)))
plt.plot(xx, 1./np.sqrt(2*np.pi*sigma**2)*np.exp(-(xx-mu)**2/(2*sigma**2)), 'r-.')


And the graphs go on top of each other:

• You're right, I don't know how I missed that crucial part of the mean of $\Delta X_t$ depending on $X_t$. Commented Nov 10, 2019 at 7:53
• Is there any advantage in considering $\sqrt 2 \Delta W_n$ where $\Delta W_n$ has variance $\Delta t$ over considering $\Delta W_n$ with variance $2\Delta t$? Commented Nov 10, 2019 at 8:38
• @StatisticalMechanic No advantage, just the convention (Brownian motion is defined with $dt$ scaling) and as such, more explicit and easier to follow. Commented Nov 10, 2019 at 8:43