Code:sde util

From FinMath

Utilities for simulation of stochastic differential equations. Lots to add here, for example "Exact Simulation of Stochastic Volatility and other Affine Jump Diffusion Processes", Broadie and Kaya (2004).

1 Demonstration
2 cirIndependent.m
3 cirSurvival.m
4 cirTimeVarying.m
5 vasIndependent.m
6 exactHeston.m

[edit] Demonstration

Cut and paste the following at the matlab prompt:

  fmuse('/sde_util','cirIndependent.m','vasIndependent.m');
  cirIndependent('demo');figure;
  vasIndependent('demo');

This will only work if you have already visited Quick start. As with all code on this site you can also manually download from www.financialmathematics.com/code/sde_util.

[edit] cirIndependent.m

Straightfoward simulation of paths for a simple square root process with jumps, drawing from non-central chi-square distribution and (optionally) an independent Poisson jump process with intensity $l$ and mean $μ$ .

$dx^i_t = \kappa \left( \theta^i - x^i_t \right) dt + \sigma^i \sqrt{x^i_t} dW^i_t + dJ^i_t$

Index $i$ might represent a path or a choice of factor, for example. Requires Matlab statistics toolbox. Returns as many cases as there are initial values supplied in $x 0$ , one in each row of the output matrix $x$ . The output does not contain a copy of the initial condition $x 0$ but rather starts with $x t = δ t$ . Thus the $k$ 'th column of output $x$ is equal to $x t = k δ t$ . If a second output argument $y$ is supplied, an estimate of the integrated process is computed by trapezoid rule with the same time convention. The $k$ 'th column of output $y$ is equal to $y t = k δ t$ .

Algorithm takes a subdivision of the interval based on the highest intensity process and may not be efficient. This function is merely a special case of cirTimeVarying.m but is easier to read and devoid of argument size checking.

function [x,y] = cirIndependent(x0,dt,nPeriods,kapp,thet,sigm,l,mu)
%
% Edit, discuss at http://www.financialmathematics.com/wiki/Code:sde_util
%
% kapp,thet,sigm,l,mu    nCases x 1 or 1 x 1 
% x0                     nCases x 1 
% dt                     1 x 1
% nPeriods               1 x 1 
% x                      nCases x nPeriods
 
if nargin==1 && isstr(x0) 
  switch x0,
      case 'demo',
         x = demo_cirIndependent;
      case 'test'
         x = test_cirIndependent;
      otherwise
         error(['Unrecognized command passed to cirIndependent']);
  end
else
  nCases = length(x0);
  x0_x = zeros(nCases,nPeriods+1);
  x0_x(:,1) = x0; 
  v = 4*kapp.*thet./(sigm.^2);
  dJ = zeros(nCases,nPeriods);
  jumping = nargin>=7;
  if jumping,
      dT = 0.1/max(l);
      nSubIntervals = ceil(dt/dT);
      dT = dt/nSubIntervals; 
      for j=1:nPeriods,
        for k=1:nSubIntervals,
          didJump = rand(nCases,1)<l*dT; 
          howHigh = -mu*log(rand(sum(didJump),1)); 
          dJ(didJump>0,j) = dJ(didJump>0,j)+ howHigh;
        end  
      end      
  end
  for j=1:nPeriods,
     ncp = 4*kapp.*exp(-kapp.*dt).*x0_x(:,j)./(sigm.^2.*(1-exp(-kapp.*dt)));
     x0_x(:,j+1) = (sigm.^2.*(1-exp(-kapp.*dt)))./(4*kapp).*ncx2rnd(v,ncp) + dJ(:,j);
  end
  x = x0_x(:,2:nPeriods+1);
  if nargout>1,
     middle = dt*0.5*(x+x0_x(:,1:nPeriods));
     y = cumsum(middle,2);
  end
end
 
function result = demo_cirIndependent
%
% Call with cirIndependent('demo')
 
[x0,dt,nPeriods,kapp,thet,sigm,l,mu] = eg_cirIndependent;
[x,y] = cirIndependent(x0,dt,nPeriods,kapp,thet,sigm,l,mu);
plot([dt:dt:dt*nPeriods],x(1:3,:)); title('3 Realizations of the CIR curve.  Demonstration cirIndependent.m is finished'); 
xlabel('Years'); 
result = NaN;
 
function [x0,dt,nPeriods,kapp,thet,sigm,l,mu] = eg_cirIndependent
nSamples = 200;dt = 1/100;nPeriods = 500; 
kapp = 1;thet = 0.1; sigm = 0.05;        
x0 = thet*ones(nSamples,1);             % Initial value
l = 0.5;mu = 0.05;                      % Jump intensity and mean
 
function result = test_cirIndependent
% Compare monte carlo and analytic mean survival
 
fmuse('/sde_util','cirSurvival.m');
[x0,dt,nPeriods,kapp,thet,sigm,l,mu] = eg_cirIndependent;
x0 = x0(1).*ones(5000,1); % Moderately large sample. 
T = nPeriods*dt;
[x,y] = cirIndependent(x0,dt,nPeriods,kapp,thet,sigm,l,mu); 
result.monteCarlo = mean(exp(-y(:,end)));
result.analytic = cirSurvival(x0(1),dt*nPeriods,kapp,thet,sigm,l,mu);

[edit] cirSurvival.m

Survival function and affine coefficient functions when the stochastic process

$dx_t = \kappa \left( \theta- x_t \right) dt + \sigma \sqrt{x_t} dW_t + dJ_t$

represents a hazard rate. Here $d J t$ is an independent Poisson jump process with $exp(l)$ distributed jumps. [Reference]. The probability of survival from time zero to time $T$ is $P(\tau > T) = E \left[ \exp( -\int_0^T x_s ds) \right]$ $= \exp\left( A + B x_0 \right)$ where [Reference]

$A = \frac{ \kappa \theta \gamma } { b c_1 d_1 } \log \left( \frac{c_1 + d_1 e^{bT} }{ -\gamma } \right) + \frac{\kappa \theta}{c_1} T + \frac{ l(a c_2 - d_2)}{ bc_2 d_2} \log \left( \frac{c_2 + d_2 e^{bT} } {c_2+d_2} \right) + \frac{l-c_2 l}{c_2} T$

and

$B = \frac{1-e^{bT}} {c_1+d_1 e^{bT}}$

Here

$\gamma = \sqrt{\kappa^2+2\sigma^2}, \ c_1 = -(\gamma+\kappa)/2$

and

$d_1 = c_1 + \kappa,\ c_2 = 1-\mu/c_1, \ b = d_1 + (\kappa c_1 - \sigma^2)/\gamma, \ a = d_1/c_1$

function [S,A,B] = cirSurvival(x0,T,kapp,thet,sigm,l,mu);
 
% Edit, discuss at http://www.financialmathematics.com/wiki/Code:sde_util
 
gamm = sqrt(kapp.^2+2.*sigm.^2);
c1 = -(gamm+kapp)/2;
d1 = c1 + kapp;
c2 = 1-mu./c1;
d2 = (d1+mu)./c1;
b = d1+(kapp*c1-sigm.^2)./gamm;
a = d1./c1;
B = (1-exp(b.*T))./(c1+d1*exp(b.*T));
term1 = (kapp.*thet.*gamm)./(b.*c1.*d1).*log( (c1+d1.*exp(b.*T) )./( -gamm )  );
term2 = T.*(kapp.*thet)./c1;
term3 = l.*(a.*c2-d2)./(b.*c2.*d2).*log(  (c2+d2.*exp(b.*T))./(c2+d2  )   );
term4 = T.*(l-c2*l)./c2;
A = term11+term2+term3+term4;
S = exp(A+B.*x0);

[edit] cirTimeVarying.m

Very similar to cirIndependent.m but allowing time and case varying parameters.

function [x,y] = cirTimeVarying(x0,dt,nPeriods,kapp,thet,sigm,l,mu)
%
% Edit, discuss at http://www.financialmathematics.com/wiki/Code:sde_util
%
% kapp,thet,sigm,l,mu    nCases x nPeriods, nCases x 1, 1 x nPeriods, or 1 x 1 
% x0                     nCases x 1 
% dt                     nCases x nPeriods, nCases x 1, 1 x nPeriods, or 1 x 1
% nPeriods               1 x 1 
% x                      nCases x nPeriods
% y                      nCases x nPeriods
 
fmuse('/argument_passing_util','makeSameSize.m');
 
if nargin==1 && isstr(x0) 
  switch x0,
      case 'demo',
         x = demo_cirTimeVarying;
      case 'test'
         x = test_cirTimeVarying;
      otherwise
         error(['Unrecognized command passed to cirTimeVarying']);
  end
else  
  nCases = length(x0);
  dummy = zeros(nCases,nPeriods);
  [dt,kapp,thet,sigm,l,mu,dummy] = makeSameSize(dt,kapp,thet,sigm,l,mu,dummy);
  x0_x = zeros(nCases,nPeriods+1);
  x0_x(:,1) = x0; 
  dJ = zeros(nCases,nPeriods);
  jumping = nargin>=7;
  if jumping, 
      for j=1:nPeriods,
        dT = 0.1/max(l(:,j));
        nSubIntervals = max(ceil(dt(:,j)./dT));
        dT = dt(:,j)./nSubIntervals;   
        for k=1:nSubIntervals,
          didJump = rand(nCases,1)<l(:,j).*dT; 
          howHigh = -mu(didJump>0,j).*log(rand(sum(didJump),1)); 
          dJ(didJump>0,j) = dJ(didJump>0,j)+ howHigh;
        end  
      end      
  end
  for j=1:nPeriods,
     v = 4*kapp(:,j).*thet(:,j)./(sigm(:,j).^2); 
     ncp = 4*kapp(:,j).*exp(-kapp(:,j).*dt(:,j)).*x0_x(:,j)./(sigm(:,j).^2.*(1-exp(-kapp(:,j).*dt(:,j))));
     x0_x(:,j+1) = (sigm(:,j).^2.*(1-exp(-kapp(:,j).*dt(:,j))))./(4*kapp(:,j)).*ncx2rnd(v,ncp) + dJ(:,j);
  end
  x = x0_x(:,2:nPeriods+1);
  if nargout>1,
     middle = dt.*0.5.*(x+x0_x(:,1:nPeriods));
     y = cumsum(middle,2);
  end
end
 
function result = demo_cirTimeVarying
%
% Call with cirTimeVarying('demo')
 
[x0,dt,nPeriods,kapp,thet,sigm,l,mu] = eg_cirTimeVarying;
[x,y] = cirTimeVarying(x0,dt,nPeriods,kapp,thet,sigm,l,mu);
plot(x(1:3,:)'); title('Demonstration cirTimeVarying.m is finished'); 
result = NaN;
 
function [x0,dt,nPeriods,kapp,thet,sigm,l,mu] = eg_cirTimeVarying
nSamples = 200;nPeriods = 500;
dt = 0.1+0.02*rand(1,nPeriods);
kapp = 0.1+rand(1,nPeriods);
thet = 0.1; 
sigm = 0.05*ones(nSamples,1);        
x0 = thet*ones(nSamples,1);             % Initial value
l = 0.5;mu = 0.05*rand(nSamples,1);     % Jump intensity and mean
 
function result = test_cirTimeVarying
% Compare monte carlo and analytic mean survival
 
fmuse('/sde_util','cirSurvival.m');
[x0,dt,nPeriods,kapp,thet,sigm,l,mu] = eg_cirTimeVarying;
x0 = x0(1).*ones(5000,1); % Moderately large sample.
[x,y] = cirTimeVarying(x0,dt,nPeriods,kapp,thet,sigm,l,mu); 
result.monteCarlo = mean(exp(-y(:,end)));
result.analytic = 1;
% TODO: Compute recursively

[edit] vasIndependent.m

Straightfoward simulation of paths for an Ornstein-Uhlenbeck process with non-zero mean and jumps.

$dx^i_t = \kappa \left( \theta^i - x^i_t \right) dt + \sigma^i dW^i_t + dJ^i_t$

Here $d J t$ is an independent Poisson jump process with $exp(l)$ distributed jumps. The joint process $(x t, y t)$ is gaussian after conditioning on jumps. Index $i$ might represent a path or a choice of factor, for example. Returns as many cases as there are initial values supplied in $x 0$ , one in each row of the output matrix $x$ . The output does not contain a copy of the initial condition $x 0$ but rather starts with $x t = δ t$ . Thus the $k$ 'th column of output $x$ is equal to $x t = k δ t$ . The $k$ 'th column of output $y$ is equal to $y t = k δ t$ .

function [x,y,yTrapezoid] = vasIndependent(x0,dt,nPeriods,kapp,thet,sigm,l,mu)
%
% Edit, discuss at http://www.financialmathematics.com/wiki/Code:sde_util
%
% kapp,thet,sigm         nCases x 1 or 1 x 1 
% x0                     nCases x 1 
% dt                     1 x 1
% nPeriods               1 x 1 
% x                      nCases x nPeriods
% y                      nCases x nPeriods
 
if nargin==1 && isstr(x0) 
  switch x0,
      case 'demo',
         x = demo_vasIndependent;
      case 'test'
         x = test_vasIndependent;
      otherwise
         error(['Unrecognized command passed to vasIndependent']);
  end
else
  nCases = length(x0);
  x0_x = zeros(nCases,nPeriods+1);
  y0_y = zeros(nCases,nPeriods+1);
  x0_x(:,1) = x0; 
  v = 4*kapp.*thet./(sigm.^2);
  dJ = zeros(nCases,nPeriods);
  intDJ = zeros(nCases,nPeriods);
  jumping = nargin>=7;
  if jumping,
      % Compute incremental effect of jumps taking mean reversion into account
      % A little rough and ready. Assumes jumps occur at mid point of sub-interval
      dT = 0.1/max(l);
      nSubIntervals = ceil(dt/dT);
      dT = dt/nSubIntervals;
      timeToEndOfPeriod = linspace(dt-dT/2,dT/2,nSubIntervals); 
      for j=1:nPeriods,
        for k=1:nSubIntervals,
          didJump = rand(nCases,1)<l*dT;
          if any(didJump),
              timeLeft = timeToEndOfPeriod(k);
              jumpEffect = exp(-kapp.*timeLeft).*didJump;                    % nCases x 1
              jumpIntegralEffect = didJump.*(1-exp(-kapp.*timeLeft))./kapp;  % nCases x 1
              howHigh = -mu.*log(rand(sum(didJump),1));
              dJ(didJump>0,j) = dJ(didJump>0,j)+ howHigh.*(jumpEffect(didJump>0));
              intDJ(didJump>0,j) = intDJ(didJump>0,j)+ howHigh.*(jumpIntegralEffect(didJump>0));
          end
        end  
      end      
  end
  % Simulate gaussian increments for x and y ignoring drift terms
  E = exp(-kapp.*dt);  % nCases x 1 or 1 x 1
  F = (1-E)./kapp;  
  xVar = sigm.^2.*(1-E.^2)./(2.*kapp);
  yVar = sigm.^2.*( (dt-F)./(kapp.^2) - F.^2./(2*kapp) );
  xyCov = (E.*sigm).^2.*(E-1).^2./(2.*kapp.^2);
  C = [xVar,xyCov;xyCov,yVar];
  nSamples = nPeriods.*nCases;
  R = mvnrnd(zeros(nSamples,2),C);
  dx = reshape(R(:,1),nCases,nPeriods);
  dy = reshape(R(:,2),nCases,nPeriods);
 
  for j=1:nPeriods,
    x0_x(:,j+1) = thet + (x0_x(:,j)-thet).*E + dJ(:,j) + dx(:,j);  
    y0_y(:,j+1) = y0_y(:,j) + intDJ(:,j) + dy(:,j) + thet.*dt + F.*(x0_x(:,j)-thet);             
  end
  x = x0_x(:,2:nPeriods+1);
  y = y0_y(:,2:nPeriods+1);
  if nargout>2,
     middle = dt*0.5*(x+x0_x(:,1:nPeriods));
     yTrapezoid = cumsum(middle,2);
  end
end
 
function result = demo_vasIndependent
% Call with vasIndependent('demo')
 
[x0,dt,nPeriods,kapp,thet,sigm,l,mu] = eg_vasIndependent;
[x,y] = vasIndependent(x0,dt,nPeriods,kapp,thet,sigm,l,mu);
plot([dt:dt:dt*nPeriods],x(1:3,:)); title('3 Realizations of the Vasicek curve.  Demonstration vasIndependent.m is finished'); 
xlabel('Years');
result = NaN;
 
function [x0,dt,nPeriods,kapp,thet,sigm,l,mu] = eg_vasIndependent
nSamples = 500;dt = 0.01;nPeriods = 500; 
kapp = 1;thet = 0.1; sigm = 0.05;        
x0 = thet*ones(nSamples,1);             % Initial value
l = 1;mu = 0.02;                      % Jump intensity and mean
 
function result = test_vasIndependent
% Call with vasIndepedent('test')
[x0,dt,nPeriods,kapp,thet,sigm,l,mu] = eg_vasIndependent;
T = nPeriods*dt;
[result.x,result.y,result.yCheck] = vasIndependent(x0,dt,nPeriods,kapp,thet,sigm,l,mu);

[edit] exactHeston.m

Stub: code will be posted soon.

This code simulates the Heston stochastic volatility model exactly. It is an implementation of the algorithm outlined in "Exact Simulation of Stochastic Volatility and other Affine Jump Diffusion Processes", Broadie and Kaya (2004).

Recall that Heston's model stipulates the following dynamics for the stock price and the variance process under the risk-neutral measure.

$dS_t = r S_tdt+\sqrt{V_t}S_t\left[\rho dW_t^{(1)}+\sqrt{1-\rho^2}dW_t^{(2)}\right]$

$dV_t = \kappa(\theta-V_t)dt+\sigma_v\sqrt{V_t}dW_t^{(1)}$

where $θ$ is the long run mean of the variance, $κ$ is the rate of meTehan reversion, $σ v$ determines the volatility of the variance process. $ρ$ is the instanteneous correlation between the return and volatility process and both $W_t^{(1)}$ and $W_t^{(2)}$ are independent Brownian motions.