Monte Carlo simulation (Matlab)

From LiteratePrograms

Contents 1 Default MATLAB functions 1.1 The simplest Monte Carlo Simulation with Matlab 1.1.1 The parameters of a Monte Carlo simulation in finance 1.1.2 Simplest Monte Carlo from scratch 1.2 Mutlivariate Monte Carlo from scratch 1.2.1 Simple Multivariate Monte Carlo 1.2.2 Tests simplest simulations 2 A more generic Monte Carlo simulation in MATLAB

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:

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

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

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, S₀ is their prices at the begining of the simulation, , σ their volatilities (annual), T is the length of the simulations (step is the time step), and nb_traj is the number of trajectories to simulate.

<<test_simplest_montecarlo0.m>>=
S0      = 1000
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)));

Simplest Monte Carlo from scratch

We will use the fact that a solution of (D.1) is:

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

Mutlivariate Monte Carlo from scratch

Obtained mutli assets simulation

Now we will simulate this kind of diffusion (for N correlated assets): 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: In this expression, δt is the time step, 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

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.

• Then the Cholesky decomposition (chol.m Matlab function) is used to be able to build a correlated multivariate gaussian noise.

• The randn function can generate multivariate independant gaussian variables, here we need nT sample of 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 S_d matrix ; the second one give them as arguments to the blkdiag function. those two lines can be read as:

• Thanks to that trick the product of our generated gaussians by S_db is:

id est:

so we went from nb_traj n-dimensionnal i.i.d. gaussian noises dX to nb_traj 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

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

A more generic Monte Carlo simulation in MATLAB

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

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 (W₁,W₂);
2. to simulate the diffusion of σ;
3. to finally obtain a simulation of S.

Download code