3

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?

Leos313
  • 5,152
  • 6
  • 40
  • 69
David
  • 77
  • 8

1 Answers1

3

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

horchler
  • 18,384
  • 4
  • 37
  • 73
  • Thank you for your comment. I did modify the function as you said but still I have the same error message: Drift rate invalid at initial conditions or inconsistent model dimensions. – David Feb 09 '16 at 09:22
  • @Mari: How are you setting your initial conditions and initial time? You need to set the `StartState` and `StartTime` properties of the `MDL` object. I have no idea what happens if you don't set these, but your drift function will return `[-Inf;NaN]` if it defaults to `[0;0]` for your initial conditions, which will be bad. – horchler Feb 12 '16 at 21:15