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?

Answers


horchler February 2016

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

Post Status

Asked in February 2016
Viewed 1,982 times
Voted 8
Answered 1 times

Search




Leave an answer