%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

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.01;
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.05; %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, 100000]*10;

%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
%(lower bond, ponctual, and upper bound)
no_obs_MC = zeros(nb_simulations,3);
%matrix that will contain the Monte Carlo estimations of the quantile
%(lower bond, ponctual, and upper bound)
x_alpha_MC = zeros(nb_simulations,3);


%*************************************
%Simulation
%*************************************
for i=1:nb_simulations
   n = nb_traj(1,i);
   
   %Simulation of the sample of size n 
   randn('state',0);
   S = s0*exp( (mu-nu^2/2)*t + nu*sqrt(t)*normrnd(0,1,n,1) );
   S_rank = sort(S);
   
   %Ponctual estimation
   no_obs_MC(i,2) = ceil( n*theta );
   x_alpha_MC(i,2) = S_rank(no_obs_MC(i,2),1);

   %Confidence interval
   no_obs_MC(i,1) = floor( n*theta - z_beta*sqrt(n*theta*(1-theta)) );
   no_obs_MC(i,3) = ceil( n*theta + z_beta*sqrt(n*theta*(1-theta)) );      
   if no_obs_MC(i,1)>0
      x_alpha_MC(i,1) = S_rank(no_obs_MC(i,1),1);
   end;
   if no_obs_MC(i,3)>0
       if no_obs_MC(i,3) < n+1
           x_alpha_MC(i,3) = S_rank(no_obs_MC(i,3),1);
       end
   end

end




%*************************************
%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.2f ',s0, mu, nu, t);
fprintf('\n');
fprintf('---------------------------------------------------------------\n\n');
fprintf('\n\n');

fprintf('No. of observations                             - Sample size\n');
fprintf('---------------------------------------------------------------\n\n');
for i=1:nb_simulations
   fprintf(' %12.0f ',no_obs_MC(i,:) );
   fprintf(' %12.0f ', nb_traj(1,i) );
   fprintf('\n');
end
fprintf('\n\n');
fprintf('If the no. of observation is smaller than one, than the sample size is too small\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(' %12.3f ',1 - beta);
fprintf('\n');
fprintf('--------------------------------------------------------------------------\n\n');
fprintf('lower bound        quantile      higher bound     - Sample size     length\n');
fprintf('--------------------------------------------------------------------------\n\n');
for i=1:nb_simulations
   fprintf(' %12.10f ',x_alpha_MC(i,:) );
   fprintf(' %12.0f ', nb_traj(1,i) );
   fprintf(' %12.10f ',x_alpha_MC(i,3)-x_alpha_MC(i,1) );
   fprintf('\n');
end
fprintf('\n\n');