Code covered by the BSD License

# Monte Carlo simulations using MATLAB

11 Dec 2007 (Updated )

Demonstrations of Monte Carlo simulations in MATLAB

WebinarScript.m
```%% Classic Monte Carlo simualtion
% Vincent Leclercq, The MathWorks, 2007, vincent.leclercq@mathworks.fr
%

clear all;
close all;

NbTrials = 10000;

%RunMode = 'LogNomrality';

RunMode = 'OptionPricing';
%% Load Data (retrieved originally from Thomson Datastream)

PastDate = today()-1 * 365;

%% Retrieve the dates in a numeric format

Dates   = cellfun(@(x)(datenum(x,'yyyy-mm-ddTHH:MM:SS')),Equities.DATE);
AssetPrices = cellfun(@(x) (str2double(x)),Equities.P );

%% Plot the Series

plot(Dates,AssetPrices);set(gcf,'WindowStyle','Docked');
legend(Equities.DISPNAME);
datetick('x','mmmyy');
xlim([PastDate, today()]);
grid on;

%% Compute the returns
% We can compute the returns for one or many stocks at the same time using
% matrix computation and matlab easy syntax

Suez_Returns = tick2ret(AssetPrices);
DailyVol = std(Suez_Returns);
AnnualVol = DailyVol * sqrt(252);

SpotPrice = AssetPrices(end);

InterestRate = 0.0375;

%% Portfolio simulation (Monte Carlo) using the Financial toolbox function
% For help, one can use the doc portsim function
% We call the portfolio simulation using Financial toolbox to
% simulate 10000 scenarios. Of course, correlation are preserved
% We assume an horizon of 6 * 22 trading days, ie 6 month maturity

%% Using Annual statistics

SimulatedRetsAnnual = portsim(InterestRate,AnnualVol^2, 12, 1/12, NbTrials,'Expected'); % NbStep * TimeStep = 1 (in years !!!)

%% Using Daily statistics
NumberOfSimulationSteps = 12;
SimulatedRetsDaily = portsim(InterestRate./252, DailyVol^2, 12, 252./12, NbTrials,'Expected');% NbStep * TimeStep = 252 (in days!!!)

%% Generate the Prices and plot them

SimulatedPricesAnnual = ret2tick(squeeze(SimulatedRetsAnnual) ,SpotPrice);
SimulatedPricesDaily  = ret2tick(squeeze(SimulatedRetsDaily)  ,SpotPrice );

figure;hist(SimulatedPricesAnnual(end,:),40);title('Prices, annual timestep used');set(gcf,'WindowStyle','Docked');
figure;hist(SimulatedPricesDaily(end,:),40);title('Prices, daily timestep used');set(gcf,'WindowStyle','Docked');

%% Check For LogNormality of the Price series

ExpectedVariance = (SpotPrice^2)  * (exp(AnnualVol^2) - 1)* exp(2*InterestRate);

disp(['Mean Price (Annual Parameters) -> ', num2str(mean(SimulatedPricesAnnual(end,:))) ' , Theoric value (Hull) :' num2str(SpotPrice * exp(InterestRate))]);
disp(['Mean Price (Daily Parameters) -> ', num2str(mean(SimulatedPricesDaily(end,:))) ' , Theoric value (Hull) :' num2str(SpotPrice * exp(InterestRate))]);

disp(['Expected Variance (Annual Parameters) -> ', num2str(var(SimulatedPricesAnnual(end,:))) ' , Theoric value (Hull) :' num2str(ExpectedVariance)]);
disp(['Expected Variance (Daily Parameters) -> ', num2str(var(SimulatedPricesDaily(end,:))) ' , Theoric value (Hull) :' num2str(ExpectedVariance)]);

%% Parameter sweep
% Now that we have done this, we can ccompute the same thing for different
% Exercise prices
if strcmp(RunMode,'OptionPricing')
k = 1;
NumberOfSteps = 400;
ExercisePrices= linspace(0.8 * SpotPrice,1.2 * SpotPrice,NumberOfSteps);

VanillaPriceAnnual    = zeros(NumberOfSteps,1);
VanillaPriceDaily    = zeros(NumberOfSteps,1);

ProbabilityITMAnnual  = zeros(NumberOfSteps,1);
ProbabilityITMDaily  = zeros(NumberOfSteps,1);
CIAnnual    =     zeros(NumberOfSteps,2);
CIDaily     =     zeros(NumberOfSteps,2);
BLSPrices       = zeros(NumberOfSteps,1);

%%
TimeInYear = 1;
for i = 1 : NumberOfSteps
[BLSPrices(i),dummy]   = blsprice(SpotPrice, ExercisePrices(i), InterestRate, TimeInYear, AnnualVol, 0);
[VanillaPriceAnnual(i), ProbabilityITMAnnual(i) ,CIAnnual(i,:)] = GetOptionPrice(SimulatedPricesAnnual,ExercisePrices(i),TimeInYear,InterestRate,'Vanilla');
[VanillaPriceDaily(i), ProbabilityITMDaily(i) ,  CIDaily(i,:)] = GetOptionPrice(SimulatedPricesDaily,ExercisePrices(i),TimeInYear,InterestRate,'Vanilla');

end;

%%

h =    figure;
[AX,H1,H2] = plotyy(ExercisePrices,[VanillaPriceAnnual BLSPrices CIAnnual] , ExercisePrices,ProbabilityITMAnnual);

xlabel('Exercise Price');
title('Option prices for a Vanilla option using a 1 year - Annual volatility');
Axes_YLabels = get(AX,'Ylabel');
set(Axes_YLabels{1},'String','Option Price') ;

set(H1(1),'LineStyle','-');
set(H1(1),'Color','r');
set(H1(1),'LineWidth',2);

set(H1(2),'LineStyle','-');
set(H1(2),'Color','b');
set(H1(2),'LineWidth',2);

set(H1(3),'LineStyle','--');
set(H1(3),'Color','r');
set(H1(3),'LineWidth',1);

set(H1(4),'LineStyle',':');
set(H1(4),'Color','r');
set(H1(4),'LineWidth',1);

set(H2,'LineStyle','-');
set(H2,'Color','g');
set(H2,'LineWidth',2);

set(get(AX(2),'Ylabel'),'String','Probability of being In the Money');
legend(H1,{['Option Price (Monte Carlo)'], ['Option Price (Black Scholes)'], ['99% Confidence interval (Lower)'],['99% Confidence interval (Upper)']},'Location','NorthEast');
legend(H1(2),{'Option Price (Black Shcoles)'},'Location','NorthEast');
legend(H2,{'Probability'},'Location','SouthWest');
grid on;
set(h,'WindowStyle','Docked');

%%
h=     figure;
[AX,H1,H2] = plotyy(ExercisePrices,[VanillaPriceDaily BLSPrices CIDaily] , ExercisePrices,ProbabilityITMDaily);

xlabel('Exercise Price');
title('Option prices for a Vanilla option using a 1 year - Daily volatility');
Axes_YLabels = get(AX,'Ylabel');
set(Axes_YLabels{1},'String','Option Price') ;

set(H1(1),'LineStyle','-');
set(H1(1),'Color','r');
set(H1(1),'LineWidth',2);

set(H1(2),'LineStyle','-');
set(H1(2),'Color','b');
set(H1(2),'LineWidth',2);

set(H1(3),'LineStyle','--');
set(H1(3),'Color','r');
set(H1(3),'LineWidth',1);

set(H1(4),'LineStyle',':');
set(H1(4),'Color','r');
set(H1(4),'LineWidth',1);

set(H2,'LineStyle','-');
set(H2,'Color','g');
set(H2,'LineWidth',2);

set(get(AX(2),'Ylabel'),'String','Probability of being In the Money');
legend(H1,{['Option Price (Monte Carlo)'], ['Option Price (Black Scholes)'], ['99% Confidence interval (Lower)'],['99% Confidence interval (Upper)']},'Location','NorthEast');
legend(H2,{'Probability'},'Location','SouthWest');
grid on;
set(h,'WindowStyle','Docked');
end;```