# Monte Carlo simulation (Matlab)

This program is under development.
is complete, remove the {{develop}} tag.

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

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


#### 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


### 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


#### 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


#### 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:

$(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.
hijacker
hijacker
hijacker
hijacker