Monte Carlo simulation (Matlab)

From LiteratePrograms
Jump to: navigation, search

This program is under development.
Please help to debug it. When debugging
is complete, remove the {{develop}} tag.


Contents


[edit] Default Matlab functions

In the Matlab Finance Toolbox, a function called portsim.m allows to perform Monte Carlo simulations. It uses the mvnrnd.m function (from the Statistics Toolbox) which returns a random vector given a chosen variance and expectation.

Those two functions can be used to perform Monte Carlo simulations of simple ITo diffusion like:

(D.1)\;\;\; dS = S\, ( r \, dt + \sigma\, (dW + \lambda\, dt))

where r, σ and λ are constants of time and S.

[edit] The simplest Monte Carlo Simulation with Matlab

Here is the simplest Monte Carlo possible using the basic functions of Matlab (to understand the mechanisms of this simulation).

[edit] The parameters of a Monte Carlo simulation in finance

The goal is the simulation of trajectories of several stochastic process defined in (D.1). S is the price of the stocks, S0 their prices at the begining of the simulation, \mu = r + \sigma\, \lambda, σ their volatilities (annual), T the length of the simulations (step the time step), and bn_traj the number of trajectories to simulate.

<<test_simplest_montecarlo0.m>>=
S0      = 100
mu      = .05
sigma   = .12
T       =   5
nb_traj = 100
step    = 1/255

[S, t]  = simplest_montecarlo0( sigma, T, nb_traj, S0, mu, step);

figure('Color',[0.9412 0.9412 0.9412 ]);
plot(t,S, 'linewidth',2);
axis([t(1) t(end) min(S(:)) max(S(:))]);
xlabel('time');ylabel('value');
title(sprintf('%d trajectories', size(S,2)));

[edit] Simplest Monte Carlo from scratch

Obtained simulation
We will use the fact that a solution of (D.1) is:
S_t = S_0\exp\left[ \left( \mu - {\sigma^2\over 2}\right) t + \sigma W_t\right]

Of course we will have to take into account the size of the time steps.

<<simplest_montecarlo0.m>>=
function [S, t] = simplest_montecarlo0( sigma, T, nb_traj, S0, mu, step)
% SIMPLEST_MONTECARLO - the simplest Monte Carlo simulation with Matlab
% use:
%  S = simplest_montecarlo( sigma, T, nb_traj, S0, mu, step)
% example:
%  [S, t] = simplest_montecarlo( .12, 5, 50, 100, .05, 1/255);

simplest_monovariate_montecarlo
<<simplest_monovariate_montecarlo>>=
nT = ceil(T/step);
W  = sigma * sqrt(step) * cumsum(randn(nT, nb_traj));
c  = repmat((mu - sigma^2/2) *step * (1:nT)',1,nb_traj);
S  = [repmat(S0,1,nb_traj); S0 * exp( c + W)];
if nargout > 1
   t = [0;step * (1:nT)'];
end

[edit] Mutlivariate Monte Carlo from scratch

Obtained mutli assets simulation

Now we will simulate this kind of diffusion (for N correlated assets):

(D.1')\;\;\; \forall \,1\leq n\leq N,\;  dS_n = S\,( r\, dt + \sigma_n\, ( dW_n + \lambda_n\, dt) ),\;\; d\langle W_k,W_l\rangle = \rho_{k,l} \, dt


Here we won't solve the diffusion process, but rather simulate a discretized version of it.

That is we will write each on the N equation of (D.1') as:

\begin{matrix}
           & S_n-S_{n-1} &=& S_{n-1}\,( \mu_n\, \delta t + \sigma_n\,  \delta W_n  ) &,\, \mu_n=r+\sigma_n\lambda_n\\
\Rightarrow& S_n         &=& S_{n-1}\,(1+\mu_n\, \delta t + \sigma_n\,  \delta W_n  )\\
\Rightarrow& S_n         &=& S_0\,\prod_{k=1}^n (1+\mu_k\, \delta t + \sigma_k\,  \delta W_k  )
\end{matrix}


In this expression, δt is the time step, (\delta W_k)_{1\leq k\leq T} is a correlated i.i.d. gaussian noise.

<<simplest_montecarlo.m>>=
function [S, t] = simplest_montecarlo( mode, sigma, T, nb_traj, S0, mu, step)
% SIMPLEST_MONTECARLO - the simplest Monte Carlo simulation with Matlab
% use:
%  S = simplest_montecarlo( mode, sigma, T, nb_traj, S0, mu, step)
% example:
%  [S, t] = simplest_montecarlo( 'simplest', .12, 5, 50, 100, .05, 1/255);

switch lower(mode)
    
    case 'test-simple'

        test_simplest_montecarlo

    case 'test-multi'

        test_simplest_multiple_montecarlo
        
    case {'simple', 'simplest'}

        simplest_monovariate_montecarlo

    case {'multi', 'multiple'}
        
        simplest_multivariate_montecarlo
end

[edit] Simple Multivariate Monte Carlo

  • At first we use the corr2cov function, which returns the standard deviation and the correlation matrix associated with a covariance matrix. The only goal (before reverting it back to the covariance matrix), is to take into account the used step.
  • The randn function can generate multivariate independant gaussian variables, here we need nT sample of n_\mbox{assets}\times n_\mbox{traj} realisations.
  • Then we use a typical Matlab trick in two lines:
Sd  = repmat({Sd},1,nb_traj);
Sdb = blkdiag(Sd{:});
the first one build a cellarray (a list) of nb_traj replications of the Sd matrix ; the second one give them as arguments to the blkdiag function. those two lines can be read as:
S_{db} = \mbox{blkdiag}\underbrace{(S_d,S_d,\ldots,S_d)}_{nb_{traj}\; times}
  • Thanks to that trick the product of our generated gaussians by Sdb is:
dW_c := \left[\begin{matrix} 
dX_1(0,1) &\cdots& dX_n(0,1) &\cdots& dX_1(0,nb_{traj}) &\cdots& dX_n(0,nb_{traj})\\
\vdots  &        &   \vdots  &\vdots&     \vdots        &      &\vdots  \\
dX_1(T,1) &\cdots& dX_n(T,1) &\cdots& dX_1(T,nb_{traj}) &\cdots& dX_n(T,nb_{traj})
\end{matrix}\right] \times \left[\begin{matrix} 
S_d & 0   &\cdots & 0\\
0   & S_d &\cdots & 0\\
\vdots & \ddots & &\vdots\\
0 & 0 & \cdots & S_d\end{matrix}\right]
id est:
dW_c=\left[\begin{matrix} 
dW_1(0,1) &\cdots& dW_n(0,1) &\cdots& dW_1(0,nb_{traj}) &\cdots& dW_n(0,nb_{traj})\\
\vdots  &        &   \vdots  &\vdots&     \vdots        &      &\vdots  \\
dW_1(T,1) &\cdots& dW_n(T,1) &\cdots& dW_1(T,nb_{traj}) &\cdots& dW_n(T,nb_{traj})
\end{matrix}\right]
so we went from nbtraj n-dimensionnal i.i.d. gaussian noises dX to nbtraj n-dimensionnal correlated i.i.d. normal noises dW.
<<simplest_multivariate_montecarlo>>=
nT  = ceil(T/step);
[s,c] = cov2corr(sigma);
s   = s * sqrt(step);
sigma = corr2cov(s,c);
Sd  = chol(sigma);
dW  = randn(nT,length(S0)*nb_traj);
Sd  = repmat({Sd},1,nb_traj);
Sdb = blkdiag(Sd{:});
dWc = dW * Sdb;
c   = repmat(1+mu*step, nT, nb_traj);
S   = [repmat(S0, 1, nb_traj); repmat(S0, nT, nb_traj) .* cumprod(c + dWc)];
S   = reshape(S, [nT+1, size(S0,2), nb_traj]);
if nargout > 1
   t = [0;step * (1:nT)'];
end

[edit] Tests simplest simulations

Each of those two blocks only call the main function. They:

  • initialize the parameters
  • call the main function in the correct mode
  • plot the obtained result
<<test_simplest_motecarlo>>=
mode    = 'simplest';
S0      = 100;
mu      = .05;
sigma   = .12;
T       = 5;
nb_traj = 100;
step    = 1/255;
        
[S, t]  = simplest_montecarlo( mode, sigma, T, nb_traj, S0, mu, step);

figure('Color',[0.9412 0.9412 0.9412 ]);
plot(t,S, 'linewidth',2);
axis([t(1) t(end) min(S(:)) max(S(:))]);
xlabel('time');ylabel('value');
title(sprintf('%d trajectories', size(S,2)));
<<test_simplest_multiple_montecarlo>>=
mode    = 'multiple';
T       = 1/255*10;
nb_traj = 2;
step    = 1/255;
S0      = [100, 150];
s       = [.12 .18];
sigma   = corr2cov(s, [1 .3; .3 1]);
mu      = [0.05 0.04];
T       = 5;
nb_traj = 100;
step    = 1/255;
        
[S, t]  = simplest_montecarlo( mode, sigma, T, nb_traj, S0, mu, step);

figure('Color',[0.9412 0.9412 0.9412 ]);
subplot(2,1,1);
plot(t,reshape(S(:,1,:),size(S,1),size(S,3)), 'linewidth',2);
axis([t(1) t(end) min(S(:)) max(S(:))]);
xlabel('time');ylabel('value');
title(sprintf('%d trajectories - stock 1 (\\sigma=%5.2f%%)', size(S,3), s(1)*100));
subplot(2,1,2);
plot(t,reshape(S(:,2,:),size(S,1),size(S,3)), 'linewidth',2);
axis([t(1) t(end) min(S(:)) max(S(:))]);
xlabel('time');ylabel('value');
title(sprintf('%d trajectories - stock 2 (\\sigma=%5.2f%%)', size(S,3), s(2)*100));

[edit] A more generic Monte Carlo simulation in Matlab

We often need to simulate diffusion process more generic than (D.1) like:

(D.2) \;\;\; \left\{ \begin{matrix} 
dS      &=& S\, ( r_t \, dt + \sigma_t\, (dW_1 + \lambda_t\, dt))\\
d\sigma &=& A_t(\sigma)\, dt + B_t(\sigma)\, dW_2
\end{matrix}\right.\,,\, d\langle W_1,W_2\rangle = \rho\, dt


To be able to perform a Monte Carlo simulation of diffusion (D.2), we will need:

  1. to simulate the 2-dimensional brownian diffusion process (W1,W2);
  2. to simulate the diffusion of σ;
  3. to finally obtain a simulation of S.
Download code
hijacker
hijacker
hijacker
hijacker