Ornstein-Uhlenbeck Process simulation bug

Question

I think I found a bug in a programm somebody posted but I can't fix it. It is about the simulation of an Ornstein-Uhlenbeck Process. The problem is from this [article][1] & and from wikipedia from my understanding. I'll quote the important part anyway:

2. Solution in terms of integral

The equation in your question is in terms of a stochastic integral

$$x_t = x_0 e^{-\theta t} + \mu (1-e^{-\theta t}) + \sigma e^{-\theta t}\int_0^t e^{\theta s} \mathrm{d}W_s$$

To obtain a numerical solution in Matlab with this, you'll need need to numerically approximate (discretize) the integral term using an SDE integration scheme like Euler–Maruyama described above:

th = 1;
mu = 1.2;
sig = 0.3;
dt = 1e-2;
t = 0:dt:2;             % Time vector
x0 = 0;                 % Set initial condition
rng(1);                 % Set random seed
W = zeros(1,length(t)); % Allocate integrated W vector
for i = 1:length(t)-1
    W(i+1) = W(i)+sqrt(dt)*exp(th*t(i+1))*randn;
end
ex = exp(-th*t);
x = x0*ex+mu*(1-ex)+sig*ex.*W;
figure;
plot(t,x);

Found the bug through using other constants:

Set dt = 0:dt:500; (not necassary but makes it easier to see) and furthermore set th = 2;

When so when we run the programm as follows (you can copy paste):

th = 2;
mu = 1.2;
sig = 0.3;
dt = 1e-2;
t = 0:dt:500;             % Time vector
x0 = 0;                 % Set initial condition
rng(1);                 % Set random seed
W = zeros(1,length(t)); % Allocate integrated W vector
for i = 1:length(t)-1
    W(i+1) = W(i)+sqrt(dt)*exp(th*t(i+1))*randn;
end
ex = exp(-th*t);
x = x0*ex+mu*(1-ex)+sig*ex.*W;
figure;
plot(t,x);
axis([0 500 0 2]);

The plot is not complete anymore. Its only drawn till the half of the points. If you use th = 3; its approximatly a third of the original plot.

Why is that? -> Problem 1. solved.

I understand the problem we get and the implementation james gave me from below works for bigger $th$ just fine.

Problem 2: I tried to modify James programm. What i want to do to is that $\mu$ is not constant anymore it is a real valued function given by: $\mu(t):= a + b * cos(c*t)$. with strictly positiv constants $a,b,c$. If you just use $\mu(t)$ in james programm it doesn't work like in my own program where we have to problem with $0 \cdot \infty$. So I'll just post my own work that you can see what I want to accomplish (this is copy pastable):

theta = 2;
sigma =500;
dt = 1e-2; 
T  = 700;    
grid = 0:dt:T;
D0 = 6000;
W = zeros(1,length(grid));     
for i = 1:length(grid)-1
     W(i+1) = W(i)+sqrt(dt)*exp(theta*grid(i+1))*randn;
end
mu= 6000 + 900* (cos(0.012*grid)); 
ex = exp(-theta*grid); 
D = D0*ex+mu.*(1-ex)+sigma*ex.*W;
figure;
plot(grid,[mu;D]);
axis([-40 T+40 0 10000]);

Answer 1

If the time vector is t = 0:dt:500; there is an overflow in exp(th*t(i+1)) in the for loop, as t gets large.

Perhaps something like this:

th = 2;
mu = 1.2;
sig = 0.3;
dt = 1e-2;
t = 0:dt:500;           % Time vector
x0 = 0;                 % Set initial condition
rng(1);                 % Set random seed
W = zeros(1,length(t)); % Allocate integrated W vector

W(1) = x0;
for i = 1 : length(t) - 1
    W(i+1) = W(i)*exp(-th*dt) + mu*(1-exp(-th*dt)) + ...
        sqrt(sig^2/(2*th)*(1-exp(-2*th*dt)))*randn;
end

figure;
plot(t,W);
axis([0 5 0 2]);

See if this works for the second part. It's borrowed from horchler's post.

theta = 2;
sigma =500;
dt = 1e-2; 
T  = 700;    
grid = 0:dt:T;
mu= 6000 + 900* (cos(0.012*grid)); 
x = zeros(1,length(grid)); 
x(1) = 6000;
for i = 1:length(grid)-1
     x(i+1) = x(i)+theta*(mu(i)-x(i))*dt+sigma*sqrt(dt)*randn;
end

figure;
plot(grid,[mu;x]);
axis([-40 T+40 0 10000]);

Answer 2

Your actual stochastic integral term (in your analytical solution) is huge. (Exercise: use the Ito isometry to compute its variance.) Its contribution to the overall quantity is not huge, because it is multiplied by a tiny factor $e^{-\theta t}$ on the outside. But that doesn't matter in computer arithmetic: to a computer, your stochastic integral is so huge and your exponential is so small that it is as if you have tried to compute $0 \cdot \infty$. There is a straightforward fix: simply rewrite the stochastic integral term as $\sigma \int_0^t e^{-\theta(t-s)} dW_s$, and estimate this convolution integral using an Euler-Maruyama type scheme.

Also, note that in this particular case where the stochastic integral is of a smooth deterministic function, one can use the Paley-Wiener integral to avoid this difficulty entirely. This technique amounts to integration by parts, and in this case it says that

$$\int_0^t e^{\theta s} dW_s = e^{\theta t} W_t - \int_0^t \theta e^{\theta s} W_s ds.$$

Using this trick, the integral to be computed after multiplying through by $e^{-\theta t}$ is an ordinary convolution of continuous functions (rather than a "formal" convolution of a function with white noise).