How can I write an SDE in Matlab?

Question

My professor would like me to solve a system similar to the following: $$ dx_i=[f_i(x_1,x_2,...x_n)]dt + g_ix_idW_i$$

Where $g_i$ are positive constants that measure the amplitude of the random perturbations, and $W_i$ are random variables normally distributed.

Im not sure how I can implement this in Matlab for ode45 to solve. What is throwing me off is the $dt$ tacked on to $f_i(x_1,...)$ and $dW_i$.

Is it as simple as coding

    function dxdt=money(t,x,a,b,c)

x1=x(1);
x2=x(2);
x3=x(3);

dx1=x1.*(x2-a)+x3 + 10*rand();
dx2=1-b*x2-x1.^2 + 10*rand;
dx3=-x1-c*x3 +10*rand();

dxdt=[dx1;dx2;dx3];

end

Answer 1

You absolutely cannot use ode45 to simulate this. ODEs and SDEs are different beasts. Also note that rand in Matlab is uniformly, not normally, distributed, so you're also not even generating proper Wiener increments. Before proceeding, you should take some time to read up on SDEs and how to solve them numerically. I recommend this paper, which includes many Matlab examples:

Desmond J. Higham, 2001, An Algorithmic Introduction to Numerical Simulation of Stochastic Differential Equations, SIAM Rev. (Educ. Sect.), 43 525–46. http://dx.doi.org/10.1137/S0036144500378302

The URL to the Matlab files in the paper won't work; use this one.

If you are familiar with ode45 you might look at my SDETools Matlab toolbox on GitHub. It was designed to be fast and has an interface that works very similarly to Matlab's ODE suite. Here is how you might code up your example using the Euler-Maruyma solver and anonymous functions:

a = ...
b = ...
c = ...
f = @(t,x)[ x(1)*(x(2)-a)+x(3);
            1-b*x(2)-x(2)^2;
           -x(1)-c*x(3)]; % drift function
g = 10; % or @(t,x)10;    % diffusion function, constant in this case
dt = 1e-3;                % time step
t = 0:dt:1;               % time vector
x0 = [...];               % 3-by-1 initial condition
opts = sdeset('RandSeed',1);  % Set random seed
y = sde_euler(f,g,t,x0,opts); % Integrate
figure;
plot(t,y);

In your example, the drift function, $g(t,x)$, doesn't depend on the state variable, $x$. This is so-called additive noise. In this case, more complex SDE integration schemes like Milstein reduce to Euler-Maruyma. More importantly, if you do use more complex diffusion functions (e.g., multiplicative noise where $g(t,x) = g_1(t)x$ like in the equation at the beginning of your question), you need to specify what SDE formulation you're using, Itô (most common in finance) or Stratonovich. My SDETools currently defaults to Stratonovich, but you can change it via opts = sdeset('SDEType','Ito');.

Answer 2

Unfortunately, I don't think you can use ode45 to solve it, whereas, since $dW$ distributes like $N(0,\sqrt{dt})$ you must use the following discretisation: $$ x_{n+1,i}=x_{n,i}+f_i(x_{n,i})\Delta t+g_i(x_{n,i})x_{n,i}\sqrt{\Delta t}\cdot\epsilon_i, $$ $$ \epsilon_i=\mbox{randn()} $$ $$ x_{n,i}\approx x_i(t) $$

Answer 3

Does an exact solution exist for your SDE? If not, you must use a numerical method.

What 7Raiden7 suggests is called the Euler-Maruyama method and is based on the truncated Ito-Taylor expansion.

Another simple numerical method would be the Milstein Scheme, which contains additional terms from the Ito-Taylor expansion. Its implementation is easy to program in MATLAB and exhibits a higher order of convergence than the Euler-Maruyama method.

These will not solve your system, but they might get you started.