%Estimation of the value at risk
%Portfolio constitute of only one risky asset
%Use of Black et Scoles model
%Use of the Gaussian approximation in the computation of the confidence intervals

%COVERAGE RATE

clear all; clc;

%*************************************
%Parameters of the Black Scholes model
%Geometrical brownian motion for the asset price
%*************************************
mu = 0.1;	%return
nu = 0.2;	%intantaneous volatility
s0 = 100;	%price of the risky asset at time t = 0

%*************************************
%Time horizon of the VaR
%*************************************
t = 1;		%maturity date

%*************************************
%Parameter of the quantile x_alpha:
%Theta = F(x_alpha)
%In this particular case, theta = alpha
%*************************************
theta = 0.0100;
z_theta = norminv(theta,0,1);
%Theoritical value of the quantile for comparaison purpose
Quantile_Th = s0*exp( (mu-nu^2/2)*t + nu*sqrt(t)*z_theta );

%*************************************
%Parameters of the confidence interval:
%Confidence level (1 - beta)
%*************************************
beta = [0.01 0.05 0.10 0.25 0.50 0.75]; %Confidence level
z_beta = norminv(1-beta/2,0,1);




%*************************************
%Parameters of the simulation
%*************************************
%Sample sizes for each simulation (row vector)
nb_traj = [100, 500, 1000, 5000, 10000, 50000];

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

%matrix that will contain the no. of the observations taken for the estimations
no_obs_MonteCarlo1 = zeros(nb_simulations,length(beta));
no_obs_MonteCarlo2 = zeros(nb_simulations,1);
no_obs_MonteCarlo3 = zeros(nb_simulations,length(beta));
%matrix that will contain the Monte Carlo estimations of the quantile
x_alpha_MonteCarlo1 = zeros(nb_simulations,length(beta));
x_alpha_MonteCarlo2 = zeros(nb_simulations,1);
x_alpha_MonteCarlo3 = zeros(nb_simulations,length(beta));


%*************************************
%Parameters of the REsimulation
%*************************************
nb_resim = 1000; %Number of iteration
CoverageRates = zeros(nb_simulations,length(beta));


for j=1:nb_resim
   j
   
   %*************************************
   %Simulation
   %*************************************
   for i=1:nb_simulations
      n = nb_traj(1,i);
   
      %Simulation of the sample of size n 
      S = s0*exp( (mu-nu^2/2)*t + nu*sqrt(t)*normrnd(0,1,n,1) );
      S_rank = sort(S);
   
      %Ponctual estimation
      no_obs_MonteCarlo2(i,1) = min( max( ceil( n*theta ), 0), n);
      if no_obs_MonteCarlo2(i,1)>0
         x_alpha_MonteCarlo2(i,1) = S_rank(no_obs_MonteCarlo2(i,1),1);
      end;
      
      for k=1:length(beta)
         %Confidence interval
         no_obs_MonteCarlo1(i,k) = min( max( floor( n*theta - z_beta(1,k)*sqrt(n*theta*(1-theta)) ), 0), n);
         no_obs_MonteCarlo3(i,k) = min( max( ceil( n*theta + z_beta(1,k)*sqrt(n*theta*(1-theta)) ), 0), n);  
      
         if no_obs_MonteCarlo1(i,k)>0
            x_alpha_MonteCarlo1(i,k) = S_rank(no_obs_MonteCarlo1(i,k),1);
         end
         if no_obs_MonteCarlo3(i,k)>0
            x_alpha_MonteCarlo3(i,k) = S_rank(no_obs_MonteCarlo3(i,k)); 
         end
         
         %Coverage rates
         if no_obs_MonteCarlo1(i,k)>0
            if abs( x_alpha_MonteCarlo2(i,1) - Quantile_Th) < (x_alpha_MonteCarlo3(i,k)-x_alpha_MonteCarlo1(i,k))/2;
               CoverageRates(i,k) = CoverageRates(i,k) + 1;
            end
         end
      end;
   end

end
CoverageRates = CoverageRates/nb_resim;

clc;

%*************************************
%Printing
%*************************************
fprintf('\n\n');
fprintf('Parameters of the geometrical brownian motion \n');
fprintf('---------------------------------------------------------------\n\n');
fprintf('     S0                mu          sigma      time horizon \n');
fprintf('\n');
fprintf('---------------------------------------------------------------\n\n');
fprintf(' %12.5f ',s0, mu, nu, t);
fprintf('\n');
fprintf('---------------------------------------------------------------\n\n');
fprintf('\n\n');


fprintf('\n\n');
fprintf('Theoritical quantile of order theta = ');
fprintf(' %12.5f ',theta);
fprintf('\n');
fprintf('---------------------------------------------------------------\n\n');
fprintf(' %12.10f ',Quantile_Th);
fprintf('\n');
fprintf('---------------------------------------------------------------\n\n');
fprintf('\n\n');

fprintf('Confidence interval of level 1 - beta  ');
fprintf('\n');
fprintf('--------------------------------------------------------------------------\n');
fprintf(' %8.3f ',1 - beta);
fprintf('\n');
fprintf('--------------------------------------------------------------------------\n\n');
fprintf('CoverageRates                                                - Sample size\n');
fprintf('Number of resimulation = ');
fprintf(' %5.0f ',nb_resim);
fprintf('\n');
fprintf('--------------------------------------------------------------------------\n\n');
for i=1:nb_simulations
   fprintf(' %8.3f ',CoverageRates(i,:) );
   fprintf(' %8.0f ', nb_traj(1,i) );
   fprintf('\n');
end
fprintf('\n\n');