Your answer is one click away!

Mari February 2016
### Stochastic Differential Equations (SDE) in 2 dimensions

I am working on stochastic differential equations for the first time. I am looking to simulate and solve a stochastic differential equations in two dimensions.

The model is as follows:

dp=F(t,p)dt+G(t,p)dW(t)

where:

- p is a 2-by-1 vector: p=(theta(t); phi(t))
- F is a column vector: F=(sin(theta)+Psi* cos(phi); Psi* cot(theta)*sin(phi))
- G is a 2-by-2 matrix: G=(D 0;0 D/sin(theta))
- Psi is a parameter and D is the diffusion constant

I wrote code as follows:

```
function MDL=gyro_2dim(Psi,D)
% want to solve for 2-by-1 vector:
%p=[theta;phi];
%drift function
F=@(t,theta,phi) [sinth(theta)+Psi.*cos(phi)-D.*cot(theta);Psi.*cot(theta).*sin(phi)];
%diffusion function
G=@(t,theta,phi) [D 0;0 D./sin(theta)];
MDL=sde(F,G)
end
```

Then I call the function with the following script:

```
params.t0 = 0; % start time of simulation
params.tend = 20; % end time
params.dt =0.1; % time increment
D=0.1;
nPeriods=10; % # of simulated observations
Psi=1;
MDL=gyro_2dim(Psi,D);
[S,T,Z]=simulate(MDL, nPeriods,'DeltaTime',params.dt);
plot(T,S)
```

When I run the code, I receive this error message:

Drift rate invalid at initial conditions or inconsistent model dimensions.

Any idea how to fix this error?

From the documentation for `sde`

:

User-defined drift-rate function, denoted by

`F`

.`DriftRate`

is a function that returns an`NVARS`

-by-1 drift-rate vector when called with two inputs:

- A real-valued scalar observation time`t`

.

- An`NVARS`

-by-1 state vector`Xt`

.

A similar specification is provided for the `Diffusion`

function. However, you're passing in the elements of your state vector as scalars and thus have three, rather than two, inputs. You can try changing your model creation function to:

```
function MDL=gyro_2dim(Psi,D)
% State vector: p = [theta;phi];
F = @(t,p)[sin(p(1))+Psi.*cos(p(2))-D.*cot(p(1));
Psi.*cot(p(1)).*sin(p(2))]; % Drift
G = @(t,p)[D 0;
0 D./sin(p(1))]; % Diffusion
MDL = sde(F,G);
MDL.StartTime = 0; % Set initial time
MDL.StartState = ... % Set initial conditions
```

I also changed `sinth(theta)`

to `sin(p(1))`

as there is no `sinth`

function. I can't test this as I don't have the Financial toolbox (few do).

Asked in February 2016

Viewed 1,982 times

Voted 8

Answered 1 times

Viewed 1,982 times

Voted 8

Answered 1 times