2

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]);
frakChris
  • 79
  • You appear to be using the code from my answer here. As @James has pointed out, there is not bug. Rather, you're running into numerical overflow. You can try to transform the equations, but it would probably be best to just simulate the SDE using Euler-Maruyama unless there is some special reason you need the the analytical solution. – horchler Jul 28 '16 at 14:45
  • @Ian: not sure if that question is addressed to me or not. I don't see where $\theta < 0$. Perhaps you're referring to the rearrangement in the last term which can indeed result in overflow. I'll modify the code in the answer to the other question to avoid this. – horchler Jul 28 '16 at 14:57
  • @horchler No, I made a minor mistake caused by expecting the named variables to mean something different from what they actually meant. Your code is "formally" correct, even though it breaks down in this parameter regime. – Ian Jul 28 '16 at 14:58
  • Thanks for your help guys. I added my problem accordingly @Ian. :) – frakChris Jul 29 '16 at 08:00
  • @frakChris Again, your problem is that your stochastic integral loop is huge and then multiplied by something tiny. I do see one difficulty with my suggestion: $e^{-\theta(t-s)}$ depends on $t$, so with my suggestion, the actual things being summed in the loop depend on $t$. In fact what you are doing is essentially a convolution. In the original form it is a "formal" convolution (convolve $e^{\theta s}$ with white noise), but upon using the Paley-Wiener trick, it is an ordinary convolution between two continuous functions. I'm sure there are standard routines for that. – Ian Jul 29 '16 at 10:59
  • Still, this is now getting sufficiently complicated that it seems easier to just directly simulate the SDE using Euler-Maruyama... – Ian Jul 29 '16 at 11:01
  • Totally agreed with Ian. – James Jul 29 '16 at 16:01
  • Minor typo in the third comment above this one: you are convolving $e^{-\theta t}$ with white noise. The choice of $t$ vs. $s$ doesn't matter of course but that minus sign does matter. – Ian Jul 29 '16 at 16:15
  • Also, I see now that in fact you cannot simply replace $\mu$ with $\mu(t)$ in the formula for the solution of the Ornstein-Uhlenbeck process, because in the derivation of that solution the assumption that $\mu$ is constant is used! Indeed I believe that $dX_t=\theta(\mu(t)-X) dt + \sigma dW_t$ will require direct numerical simulation even for fairly simple $\mu$. – Ian Jul 29 '16 at 16:18
  • Hello @Ian. The $\mu(t)$ i used is from the closed solution of the problem. For better unterstanding you should check my other post where somebody solved the problem for a given $\tilde{\mu}(t)$. See link where $\bar{\mu}(t)$ is our $\mu(t)$ in this post. In fact you helped me on that post too. Anyway its strange that i can't use James programm for my problem? – frakChris Jul 29 '16 at 18:31

2 Answers2

3

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]);
James
  • 1,392
  • I actually need th = 500 for some given approximation. How can i compensated that overflow? Thank you sir for your answer @James – frakChris Jul 28 '16 at 06:36
  • 1
    The original for loop is not quite right.. – James Jul 28 '16 at 14:40
  • thanks for answer. I dont understand though what you really changed. Everything you did was put the terms after in the loop? Sorry :)

    Even though im not so sure why you have the term \begin{align} W(i) \cdot exp(-th \cdot dt) \end{align} in the loop?

     – frakChris Jul 28 '16 at 14:47
  • It was Monte Carlo simulation to the OU process \begin{align} x_{t} & \sim N\left(\mu+e^{-\theta t}\left(x_{0}-\mu\right),\frac{\sigma^{2}}{2\theta}\left(1-e^{-2\theta t}\right)\right). \end{align} The other methods as pointed out are fine too. – James Jul 28 '16 at 15:09
  • I see. My first Problem is solved. Please check above for a short look at my secound problem. edited accordingly! :)! @James – frakChris Jul 29 '16 at 07:57
3

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

Ian
  • 104,572