%
%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 = 1;		%maturity date

%Strike prices (row vector)
K = [0.9, 1, 1.1]*s0; 
K = [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; 
%number of simulations we perform
nb_sim = length(nb_traj);
%number of resimulation for the estimation of the variance of the moment
%matching methods
nb_resim = 100;

%%
%*************************************
%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('Margin of error')
    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')


%%
close all;
%*************************************
%Monte Carlo simulation + antithetic variable
%*************************************
price_Anti = zeros(nb_sim,length_K);     %matrix that will contain the Monte Carlo estimations of the options' prices
std_price_Anti = zeros(nb_sim,length_K); %matrix that will contain the sample standard deviation of Monte Carlo estimations
precision_Anti = zeros(nb_sim,length_K); %matrices that will contain 1/2 of the length of the 95% confidence intervals
time_Anti = 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
   Z=randn(n,1);            %noise
                            %(n x length_K) options discounted cashflows
   C1 = DF*max( (s0*exp( (r-0.5*sigma*sigma)*T + vol*Z ))*ones(1,length_K) - ones(n,1)*K , zeros(n,length_K) );
   C2 = DF*max( (s0*exp( (r-0.5*sigma*sigma)*T - vol*Z ))*ones(1,length_K) - ones(n,1)*K , zeros(n,length_K) );
   price_Anti(i,:) = mean(0.5*(C1+C2));     %Monte carlo option prices
   std_price_Anti(i,:) = std(0.5*(C1+C2));  %Sample standard deviation
   precision_Anti(i,:) = 1.96*std_price_Anti(i,:)/sqrt(n);
   time_Anti(i,:) = toc;                    %Computation time
   clear Z C1 C2
end

%Plots - precision
%*************************************
figure
for j = 1:length_K
    subplot(length_K,1,j);
    plot(precision_MC(:,j),time_MC,'--b*')
    hold on
    plot(precision_Anti(:,j),time_Anti,'--ro')
    hold off
    xlabel('Margin of error')
    ylabel('Computation time')
    title(['Precision of the Monte Carlo estimators of a Call price with K =', num2str(K(1,j))])
end
    legend('Naive','Antithetic')

%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*')
    plot(log(nb_traj),price_Anti(:,j),'--ro')
    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','Antithetic')


%%
close all;
%*************************************
%Monte Carlo simulation + control variable
%*************************************
price_Cont = zeros(nb_sim,length_K);     %matrix that will contain the Monte Carlo estimations of the options' prices
std_price_Cont = zeros(nb_sim,length_K); %matrix that will contain the sample standard deviation of Monte Carlo estimations
precision_Cont = zeros(nb_sim,length_K); %matrices that will contain 1/2 of the length of the 95% confidence intervals
time_Cont = zeros(nb_sim,1);             %column vector that will contain the computation times
for i=1:nb_sim
   tic;
   n = nb_traj(1,i);
   Exp_S = s0 / DF;     %Expected value of the stock price (for the control variable)
   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 (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
   betas=ones(length_K,1); %Estimation of the Betas
   for j=1:length_K
        A = cov( S,C(:,j) ) / cov(S);  betas(j,1) = A(1,2);
   end
   price_Cont(i,:) = mean(C - (S-Exp_S)*(betas'));     %Monte carlo option prices
   std_price_Cont(i,:) = std(C - (S-Exp_S)*(betas'));  %Sample standard deviation
   precision_Cont(i,:) = 1.96*std_price_Cont(i,:)/sqrt(n);
   time_Cont(i,:) = toc;                                  %Computation time
   clear A S C betas;
end

%Plots - precision
%*************************************
figure
for j = 1:length_K
    subplot(length_K,1,j,'align');
    plot(precision_MC(:,j),time_MC,'--b*')
    hold on
    plot(precision_Anti(:,j),time_Anti,'--ro')
    plot(precision_Cont(:,j),time_Cont,'--gs')
    hold off
    xlabel('Margin of error')
    ylabel('Computation time')
    title(['Precision of the Monte Carlo estimators of a Call price with K =', num2str(K(1,j))])
end
    legend('Naive','Antithetic','Control')

%Plots - prices
%*************************************
figure
for j = 1:length_K
    subplot(length_K,1,j,'align');
    plot(log(nb_traj),Th_Prices(1,j)*ones(nb_sim,1),'-k')
    hold on
    plot(log(nb_traj),price_MC(:,j),'--b*')
    plot(log(nb_traj),price_Anti(:,j),'--ro')
    plot(log(nb_traj),price_Cont(:,j),'--gs')
    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','Antithetic','Control')
    
%%
close all;
%*************************************
%Monte Carlo simulation + antithetic variable + control variable
%*************************************
price_Cont_Anti = zeros(nb_sim,length_K);     %matrix that will contain the Monte Carlo estimations of the options' prices
std_price_Cont_Anti = zeros(nb_sim,length_K); %matrix that will contain the sample standard deviation of Monte Carlo estimations
precision_Cont_Anti = zeros(nb_sim,length_K); %matrices that will contain 1/2 of the length of the 95% confidence intervals
time_Cont_Anti = zeros(nb_sim,1);             %column vector that will contain the computation times
for i=1:nb_sim
   tic;
   n = nb_traj(1,i);
   Exp_S = s0 / DF;     %Expected value of the stock price (for the control variable)
   randn('state',0);    %initialize the normal random number generator
   Z = randn(n,1);      %noise terms
   S1 = s0*exp( (r-0.5*sigma*sigma)*T + vol*Z );  %Simulated sample of stock prices (of size n)
   S2 = s0*exp( (r-0.5*sigma*sigma)*T - vol*Z );  %Simulated sample of stock prices (of size n) 
   C1 = DF*max( S1*ones(1,length_K) - ones(n,1)*K , zeros(n,length_K) );  %(n x length_K) options discounted cashflows
   C2 = DF*max( S2*ones(1,length_K) - ones(n,1)*K , zeros(n,length_K) );  %(n x length_K) options discounted cashflows
   betas1=ones(length_K,1); %Estimation of the Betas
   betas2=ones(length_K,1); %Estimation of the Betas
   for j=1:length_K
        A = cov( S1,C1(:,j) ) / cov(S1);  betas1(j,1) = A(1,2);
        A = cov( S2,C2(:,j) ) / cov(S2);  betas2(j,1) = A(1,2);
   end
   price_Cont_Anti(i,:) = mean(0.5*(C1 - (S1-Exp_S)*(betas1')+ C2 - (S2-Exp_S)*(betas2')));    %Monte carlo option prices
   std_price_Cont_Anti(i,:) = std(0.5*(C1 - (S1-Exp_S)*(betas1')+ C2 - (S2-Exp_S)*(betas2')));  %Sample standard deviation
   precision_Cont_Anti(i,:) = 1.96*std_price_Cont_Anti(i,:)/sqrt(n);
   time_Cont_Anti(i,:) = toc;  %Computation time
   clear Z S1 S2 C1 C2 A betas1 betas2;
end

%Plots - precision
%*************************************
figure
for j = 1:length_K
    subplot(length_K,1,j,'align');
    plot(precision_MC(:,j),time_MC,'--b*')
    hold on
    plot(precision_Anti(:,j),time_Anti,'--ro')
    plot(precision_Cont(:,j),time_Cont,'--gs')
    plot(precision_Cont_Anti(:,j),time_Cont_Anti,'--cs')
    hold off
    xlabel('Margin of error')
    ylabel('Computation time')
    title(['Precision of the Monte Carlo estimators of a Call price with K =', num2str(K(1,j))])
end
    legend('Naive','Antithetic','Control','Anti + control')

%Plots - prices
%*************************************
figure
for j = 1:length_K
    subplot(length_K,1,j,'align');
    plot(log(nb_traj),Th_Prices(1,j)*ones(nb_sim,1),'-k')
    hold on
    plot(log(nb_traj),price_MC(:,j),'--b*')
    plot(log(nb_traj),price_Anti(:,j),'--ro')
    plot(log(nb_traj),price_Cont(:,j),'--gs')
    plot(log(nb_traj),price_Cont_Anti(:,j),'--cs')
    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','Antithetic','Control','Anti + control')

    
%%
close all;
%*************************************
%Monte Carlo simulation + EMS
%*************************************
price_EMS = zeros(nb_sim,length_K);     %matrix that will contain the Monte Carlo estimations of the options' prices
std_price_EMS = zeros(nb_sim,length_K); %matrix that will contain the sample standard deviation of Monte Carlo estimations
precision_EMS = zeros(nb_sim,length_K); %matrices that will contain 1/2 of the length of the 95% confidence intervals
time_EMS = zeros(nb_sim,1);             %column vector that will contain the computation times
for i=1:nb_sim
   tic;
   n = nb_traj(1,i);
   Exp_S = s0 / DF;     %Expected value of the stock price (for the control variable)
   randn('state',0);    %initialize the normal random number generator
   Z = randn(n,1);      %noise terms
   S = s0*exp( (r-0.5*sigma*sigma)*T + vol*Z );  %Simulated sample of stock prices (of size n)
   C = max( S*ones(1,length_K) - ones(n,1)*K , zeros(n,length_K) ); 
   Sstar = Exp_S*S/mean(S);
   Cstar = max( Sstar*ones(1,length_K) - ones(n,1)*K , zeros(n,length_K) );  %(n x length_K) options discounted cashflows
   Phi = mean( (S*ones(1,length_K)>ones(n,1)*K) .* (S*ones(1,length_K))) / Exp_S;
   term1 = (std(C)).^2;
   term2 = (Phi.*std(S)).^2;
   term3 = Phi .* ( mean( C .* (S*ones(1,length_K)) ) - mean(C).*mean(S) );
   price_EMS(i,:) = DF*mean(Cstar);    %Monte carlo option prices
   std_price_EMS(i,:) = DF*sqrt(term1 + term2 -2*term3);  %Sample standard deviation
   precision_EMS(i,:) = 1.96*std_price_EMS(i,:)/sqrt(n);
   time_EMS(i,:) = toc;  %Computation time
   clear Z S C Sstar Cstar term1 term2 term3 Phi Exp_S;
end

%Plots - precision
%*************************************
figure
for j = 1:length_K
    subplot(length_K,1,j,'align');
    plot(precision_MC(:,j),time_MC,'--b*')
    hold on
    plot(precision_Anti(:,j),time_Anti,'--ro')
    plot(precision_Cont(:,j),time_Cont,'--gs')
    plot(precision_Cont_Anti(:,j),time_Cont_Anti,'--cs')
    plot(precision_EMS(:,j),time_EMS,':mx')
    hold off
    xlabel('Margin of error')
    ylabel('Computation time')
    title(['Precision of the Monte Carlo estimators of a Call price with K =', num2str(K(1,j))])
end
    legend('Naive','Antithetic','Control','Anti + control','EMS')

%Plots - prices
%*************************************
figure
for j = 1:length_K
    subplot(length_K,1,j,'align');
    plot(log(nb_traj),Th_Prices(1,j)*ones(nb_sim,1),'-k')
    hold on
    plot(log(nb_traj),price_MC(:,j),'--b*')
    plot(log(nb_traj),price_Anti(:,j),'--ro')
    plot(log(nb_traj),price_Cont(:,j),'--gs')
    plot(log(nb_traj),price_Cont_Anti(:,j),'--cs')
    plot(log(nb_traj),price_EMS(:,j),':mx')
    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','Antithetic','Control','Anti + control','EMS')   
    
    
 %%
close all;
%*************************************
%Monte Carlo simulation + importance sampling - ad hoc choice of the
%importance density
%*************************************
clc

price_IM = zeros(nb_sim,length_K);     %matrix that will contain the Monte Carlo estimations of the options' prices
std_price_IM = zeros(nb_sim,length_K); %matrix that will contain the sample standard deviation of Monte Carlo estimations
precision_IM = zeros(nb_sim,length_K); %matrices that will contain 1/2 of the length of the 95% confidence intervals
time_IM = 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
   Z = randn(n,1);                                                                      %Standard normal vector
   th = 1 + (s0./K).*(s0 < K) - (K/s0).*(s0 > K);                                       %ad hoc choice of the importance density   
   L = exp(-sqrt(T)*Z*th-0.5*ones(n,1)*(th.*th)*T);                                     %Likelihood ratio - Randon Nikodym derivative    
   S = s0*exp( (r-0.5*sigma*sigma)*T + ones(n,1)*sigma*th*T+ vol*Z*ones(1,length_K) );  %Simulated sample of stock prices
   C = DF*max( S - ones(n,1)*K , zeros(n,length_K) );                                   %(n x length_K) options discounted cashflows

   price_IM(i,:) = mean(C.*L);                                                          %Monte carlo option prices
   std_price_IM(i,:) = std(C.*L);                                                       %Sample standard deviation
   precision_IM(i,:) = 1.96*std_price_IM(i,:)/sqrt(n);
   time_IM(i,:) = toc;                                                                  %Computation time
   clear Z S C th L;
end

%Plots - precision
%*************************************
figure
for j = 1:length_K
    subplot(length_K,1,j,'align');
    plot(precision_MC(:,j),time_MC,'--b*')
    hold on
    plot(precision_Anti(:,j),time_Anti,'--ro')
    plot(precision_Cont(:,j),time_Cont,'--gs')
    plot(precision_Cont_Anti(:,j),time_Cont_Anti,'--cs')
    plot(precision_EMS(:,j),time_EMS,':mx')
    plot(precision_IM(:,j),time_IM,':b<')
    hold off
    xlabel('Margin of error')
    ylabel('Computation time')
    title(['Precision of the Monte Carlo estimators of a Call price with K =', num2str(K(1,j))])
end
    legend('Naive','Antithetic','Control','Anti + control','EMS','IS ad hoc')

%Plots - prices
%*************************************
figure
for j = 1:length_K
    subplot(length_K,1,j,'align');
    plot(log(nb_traj),Th_Prices(1,j)*ones(nb_sim,1),'-k')
    hold on
    plot(log(nb_traj),price_MC(:,j),'--b*')
    plot(log(nb_traj),price_Anti(:,j),'--ro')
    plot(log(nb_traj),price_Cont(:,j),'--gs')
    plot(log(nb_traj),price_Cont_Anti(:,j),'--cs')
    plot(log(nb_traj),price_EMS(:,j),':mx')
    plot(log(nb_traj),price_IM(:,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','Antithetic','Control','Anti + control','EMS','IS ad hoc')  
    
    
 
 %%
close all;
%*************************************
%Monte Carlo simulation + Fu & Su
%*************************************
clc
N1 = 20; N2 = 500; epsilon = 0.001;

price_FS = zeros(nb_sim,length_K);     %matrix that will contain the Monte Carlo estimations of the options' prices
std_price_FS = zeros(nb_sim,length_K); %matrix that will contain the sample standard deviation of Monte Carlo estimations
precision_FS = zeros(nb_sim,length_K); %matrices that will contain 1/2 of the length of the 95% confidence intervals
time_FS = zeros(nb_sim,1);             %column vector that will contain the computation times

for i=1:nb_sim
   tic; 
   
   %Choice of the optimal theta
      lambda = (log(K/s0))/T;  %initialization of lambda
      cond = 0; n1 = 0;
      while cond<0.5;       %Verification of the stopping criteria
         n1 = n1 + 1;
         W = sqrt(T)*randn(N2,1)*ones(1,length_K);
         S = s0*exp( ones(N2,1)*(lambda-0.5*sigma*sigma)*T - sigma*W);
         C = DF*max( S - ones(N2,1)*K, zeros(N2,length_K)); %option cashflow
         th = ones(N2,1)*(lambda - r)/sigma;                %theta
         gn = -mean( C.*C.*W.*exp(-2*th.*W - th.*th*T) );   %gradient estimation
         if n1==1
           a0 = 1./abs(gn);
         end
       %update
       an = n1^(-0.75).*a0;
       lambda_new = lambda + sigma*(an.*gn);  %th_new = th + an*gn. Hence, lamb_new - r = lamb - r + sig*an*gn
       if abs(lambda-lambda_new)<0.25
           lambda=lambda_new;
           n1
       end
       if n1>N1
           cond = 1;
       end
       if max(abs(an.*gn))<epsilon
           cond = 1;
       end
   end
   n1
   lambda
   
   %tic
   n = nb_traj(1,i);
   randn('state',0); %initialize the normal random number generator
   th = ones(n,1)*(lambda - r)/sigma;
   W = sqrt(T)*randn(n,1)*ones(1,length_K);                               %Standard normal vector                                   %ad hoc choice of the importance density   
   L = exp(-W.*th - 0.5*(th.*th)*T);                        %Likelihood ratio - Randon Nikodym derivative    
   S = s0*exp( ones(n,1)*(lambda-0.5*sigma*sigma)*T + sigma*W);     %Simulated sample of stock prices
   C = DF*max( S - ones(n,1)*K , zeros(n,length_K) );             %(n x length_K) options discounted cashflows

   price_FS(i,:) = mean(C.*L);                                    %Monte carlo option prices
   std_price_FS(i,:) = std(C.*L);                                 %Sample standard deviation
   precision_FS(i,:) = 1.96*std_price_FS(i,:)/sqrt(n);
   time_FS(i,:) = toc;                                            %Computation time
   %clear Z S C L N1 N2 an gn;
end
   
%Plots - precision
%*************************************
figure
for j = 1:length_K
    subplot(length_K,1,j,'align');
    plot(precision_MC(:,j),time_MC,'--b*')
    hold on
    plot(precision_Anti(:,j),time_Anti,'--ro')
    plot(precision_Cont(:,j),time_Cont,'--gs')
    plot(precision_Cont_Anti(:,j),time_Cont_Anti,'--cs')
    plot(precision_EMS(:,j),time_EMS,':mx')
    plot(precision_IM(:,j),time_IM,':b<')
    plot(precision_FS(:,j),time_FS,':r<')
    hold off
    xlabel('Margin of error')
    ylabel('Computation time')
    title(['Precision of the Monte Carlo estimators of a Call price with K =', num2str(K(1,j))])
end
    legend('Naive','Antithetic','Control','Anti + control','EMS','IS ad hoc','Fu Su')

%Plots - prices
%*************************************
figure
for j = 1:length_K
    subplot(length_K,1,j,'align');
    plot(log(nb_traj),Th_Prices(1,j)*ones(nb_sim,1),'-k')
    hold on
    plot(log(nb_traj),price_MC(:,j),'--b*')
    plot(log(nb_traj),price_Anti(:,j),'--ro')
    plot(log(nb_traj),price_Cont(:,j),'--gs')
    plot(log(nb_traj),price_Cont_Anti(:,j),'--cs')
    plot(log(nb_traj),price_EMS(:,j),':mx')
    plot(log(nb_traj),price_IM(:,j),':b<')
    plot(log(nb_traj),price_FS(:,j),':r<')
    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','Antithetic','Control','Anti + control','EMS','IS ad hoc', 'Fu et Su')  
    
%price = [Th_Prices; price_MC; price_Anti; price_Cont; price_Cont_Anti; price_EMS; price_IM; price_FS]