%%
%Monte Carlo estimators for the prices of call options

%clear memory and screen, close the figures' windows
clear all; clc; close all;

%*************************************
%Parameters of the Black Scholes model
%*************************************
r = 0.05;       %risk free rate
sigma = 0.2;	%intantaneous volatility
s0 = 100;       %price of the risky asset at time t = 0

%*************************************
%Parameters of the call option
%*************************************
T = 0.25;		%maturity date

%Strike prices (row vector)
K = [0.9, 1, 1.1]*s0; %K = [0.5, 0.8, 0.9, 1, 1.1, 1.2]*s0; 

%Usefull constants that we would like to compute only once
DF = exp(-r*T);                 %Risk free discount factor
vol = sigma*sqrt(T);

%Theoritical prices of the call options (row vector)
length_K = length(K);
d = (log(s0)*ones(1,length_K) - log(K) + (r+sigma^2/2)*T*ones(1,length_K)) / vol;
Th_Prices = s0*normcdf(d) - K*DF.* normcdf(d - vol*ones(1,length_K));

%Theritical standard deviations of the prices of the call options
term1 = s0^2 * exp(sigma^2*T) * normcdf(d+vol);
term2 = 2*s0*DF * K .* normcdf(d);
term3 = K.^2 * exp(-2*r*T) .* normcdf(d-vol);
Th_std = sqrt(term1 - term2 + term3 - Th_Prices.^2); 
clear term1 term2 term3 d;

%*************************************
%Parameters of the simulation
%*************************************
%Sample sizes for each simulation (row vector)
nb_traj = [5000, 7500, 10000, 12500, 15000, 17500, 20000, 25000]*10; %, 500000, 75000, 100000, 250000, 500000];

%number of simulations we perform
nb_sim = length(nb_traj);

%%
%*************************************
%Naive Monte Carlo simulation
%*************************************
price_MC = zeros(nb_sim,length_K);      %matrix that will contain the Monte Carlo estimations of the options' prices
std_price_MC = zeros(nb_sim,length_K);  %matrix that will contain the sample standard deviation of Monte Carlo estimations
precision_MC = zeros(nb_sim,length_K);  %matrices that will contain 1/2 of the length of the 95% confidence intervals
time_MC = zeros(nb_sim,1);              %column vector that will contain the computation times
for i=1:nb_sim
   tic;
   n = nb_traj(1,i);
   randn('state',0);                                        %initialize the normal random number generator 
   S = s0*exp( (r-0.5*sigma*sigma)*T + vol*randn(n,1) );    %Simulated sample of stock prices (vector of size n)
   C = DF*max( S*ones(1,length_K) - ones(n,1)*K , zeros(n,length_K) ); %(n x length_K) options discounted cashflows
   price_MC(i,:) = mean(C);                                 %Monte carlo option prices
   std_price_MC(i,:) = std(C);                              %Sample standard deviation
   precision_MC(i,:) = 1.96*std_price_MC(i,:)/sqrt(n);
   time_MC(i,:) = toc;                                      %Computation time
   clear S C;
end

%Plots - precision
%*************************************
figure
for j = 1:length_K
    subplot(length_K,1,j);
    plot(precision_MC(:,j),time_MC,'--b*')
    xlabel('Precision')
    ylabel('Computation time')
    title(['Precision of the Monte Carlo estimators of a Call price with K =', num2str(K(1,j))])
end

    legend('Naive')

%Plots - prices
%*************************************
figure
for j = 1:length_K
    subplot(length_K,1,j);
    plot(log(nb_traj),Th_Prices(1,j)*ones(nb_sim,1),'-k')
    hold on
    plot(log(nb_traj),price_MC(:,j),'--b*')
    hold off
    xlabel('logarithm of the number of trajectories')
    ylabel('Prices')
    title(['Monte Carlo estimators of a Call price with K =', num2str(K(1,j))])
end
    legend('True','Naive')