Search code examples
matlabstochastic-process

Extract the values for some parameters from stochastic differential equation solution


I am solving stochastic differential equation in matlab. For example: consider the stochastic differential equation

dx=k A(x,t)dt+  B(x,t)dW(t)

where k is constants, A and B are functions, and dW(t) is Wiener process.

I plot the solution for all t in [0,20]. We know that dW(t) is randomly generated. My question is: I want to know the value of A(x,t), B(x,t), dW(t) for a particular value of t and for particular sub-interval, say [3,6]. What command in Matlab I can use?

Here is the code I used based on a paper by D.Higham:

clear all
close all

t0   = 0;     % start time of simulation
tend = 20;   % end time

m=2^9; %number of steps in each Brownian path

deltat= tend/m; % time increment for each Brownian path
D=0.1; %diffsuion 
R=4;
dt =  R*deltat; 
dW=sqrt( deltat)*randn(2,m);
theta0=pi*rand(1);
phi0=2*pi*rand(1);
P_initial=[ theta0; phi0];

L =  m/ R;
pem=zeros(2,L);
EM_rescale=zeros(2,L);
ptemp=P_initial;

for j=1:L
   Winc = sum(dW(:,[ R*(j-1)+1: R*j]),2);
   theta=ptemp(1);% updating theta
   phi=ptemp(2); % updating phi
   %psi=ptemp(3); % updating psi

   A=[  D.*cot(theta);...
       0];% updating the drift

   B=[sqrt(D)  0 ;...
      0 sqrt(D)./sin(theta) ]; %% updating the diffusion function

   ptemp=ptemp+ dt*A+B*Winc;
   pem(1,j)=ptemp(1);%store theta
   pem(2,j)=ptemp(2);%store phi

   EM_rescale(1,j)=mod(pem(1,j),pi); % re-scale theta
   EM_rescale(2,j)=mod(pem(2,j),2*pi);   % re-scale phi
end

plot([0:dt:tend],[P_initial,EM_rescale],'--*')  

Suppose I want to know all parameters (including random: Brownian) at each specific time point or for any time interval. How to do that?


Solution

  • I'm doing my best to understand your question here, but it's still a bit unclear to me.

    Change the loop to:

    for ii=1:L
       Winc = sum(dW(:,[ R*(ii-1)+1: R*ii]),2);
    
       theta=ptemp(1);% updating theta
       phi=ptemp(2); % updating phi
    
       A{ii}=[  D.*cot(theta);...
           0];% updating the drift
    
       B{ii}=[sqrt(D)  0 ;...
          0 sqrt(D)./sin(theta) ]; %% updating the diffusion function
    
       ptemp = ptemp + dt*A{ii}+B{ii}*Winc;
       pem(:,ii) = ptemp;
    
       EM_rescale(1,ii) = mod(pem(1,ii),pi); % re-scale theta
       EM_rescale(2,ii) = mod(pem(2,ii),2*pi);   % re-scale phi
    
    end
    

    Now, you can get the values of A and B this way:

    t = 3;
    t_num = round(m/tend*t);
    
    A{t_num}
    B{t_num}
    ans =
            0.0690031455719538
                             0
    ans =
             0.316227766016838                         0
                             0          0.38420611784333