% ------------------------------------------------------------------------
% UNIVERSIDAD NACIONAL AUTÓNOMA DE MÉXICO
% FACULTAD DE CIENCIAS
% ------------------------------------------------------------------------
% This code is part from Actuarial Science Bachelor's dissertation 
% ------------------------------------------------------------------------
% Objective:
% Stochastic simulation code using Metropolis-Hastings algorithm with 
% independence chain for ARCH(2) model. The volatility posterior samples 
% are implemented in the Black & Scholes model to come by samples of 
% posterior distribution for call or put prices.
% ------------------------------------------------------------------------
% Reasoning:
% First, the parameters (alpha0, alpha1 and alpha2) in ARCH(2) model are 
% estimated for the financial time series used. The corresponding standard 
% errors are also estimated. Second, once volatility is estimated we have
% the first state in our Markov chain. So, a MCMC technique for the bayesian
% outcomes is developed. Finally, all the results are presented. 
% ------------------------------------------------------------------------
% Outcomes:
% A) Samples from posterior distribution for the volatility. 
% B) Samples from posterior distribution for Black & Scholes model. 
% C) Figures
% ------------------------------------------------------------------------
% INTRUCTIONS: 
% The function is called bsmcmc(c,iter) where c represents
% the calibration constant for the variance-covariance matrix and iter
% is the number of iterations for the stochastic chain.
% ***In this code you should read your own data and calibrating your own
% variance-covariance matrix ***
% ------------------------------------------------------------------------
% This code was written by Carlos Vladimir Rodríguez Caballero with the
% theoretical help of Alejandro Villagrán-Hernández in 2005.
% Questions or comments? vladimir.rodriguez@ciencias.unam.mx
% ------------------------------------------------------------------------


function Results=bsmcmc(c,iter);

x=textread('ipc200.txt');       % Mexican Stock Market index
                                % You should change your data
ipc=x(2500:2750);		% Date range 

y=price2ret(ipc);               % Return series
n=length(y);                    % Length of return series

spec=garchset('P',0,'Q',2,'Display','off');   % Specification Model 
[Coeff,Errors,LLF,Innovations,Sigma,Summary]=garchfit(spec,y); 
                                              
alpha0= garchget(Coeff,'k');                  % alpha 0 estimation
alphas= garchget(Coeff,'ARCH');
alpha1= alphas(1);                            % alpha 1 estimation
alpha2= alphas(2);                            % alpha 2 estimation

s0=garchget(Errors,'k');	% Statistical errors
ss=garchget(Errors,'ARCH');
s1=ss(1);                          
s2=ss(2);                         

cov=[s0,0,0;0,s1,0;0,0,s2]^2;       % Variance Covariance Matrix 
ccov=c.*cov;               % Calibration Var-Cov Matrix 
			   % You should calibrate your own matrix

mu=[alpha0,alpha1,alpha2]';             
fuera=0;                   % Initial counter for the reject rate
cont=0;                    % Initial counter for the acceptance rate
hp=Sigma.^2;

% Volatility simulation considering an ARCH(2) process

h(1)=alpha0;
h(2)=alpha0+alpha1*y(1)^2;
for t=3:n        
      h(t)=alpha0+alpha1*y(t-1)^2+alpha2*y(t-2)^2;     
end

 ha=h';
 thetap=mu;
 T(:,1)=mu;                % First step in the stochastic chain

% ---------------- 
% From here the MCMC technique is developed
% The algoritm used is the Independence Metropolis-Hastings 
% ----------------

 for i=2:iter        
    u = unifrnd(0,1);      	% Uniform random variable generator                                
    v = log(u);            
    Theta = [mvnrnd(mu,ccov,1)]';     % Instrumental distribution
    
% In this work the proposal distribution is defined as
% trivariate normal distribution. The mean vector and covariance
% matrix have already defined before in the code. 
    
% Here, the parameters are defined for the actual and past state 
% in the Markov Chain. The subindex "a" represents the actual state
% and the subindex "p" represents the last state.

    sum1=sum(ha);      
    sum2=sum(hp);
    sum3=sum((y.^2)./ha);
    sum4=sum((y.^2)./hp);
    hp=ha;
    thetap=T(:,i-1);                 
    Va=(Theta-mu)'*inv(ccov)*(Theta-mu);
    Vp=(thetap-mu)'*inv(ccov)*(thetap-mu);
        
% Here, the ARCH(2) model restrictions are displayed and the logratio
% is defined.

    if (Theta(1)>0)&(Theta(2)>0)&(Theta(3)>0)&(Theta(2)<1)&(Theta(3)<1)
        logratio = -0.25*(sum1-sum2+sum3-sum4+Vp-Va);    
    else 
     logratio=0;
     fuera=fuera+1;
end
   
% Here, the Metropolis-Hastings condition is implemented  

  if v<=min(logratio,0)
    if (Theta(1)>0)&(Theta(2)>0)&(Theta(3)>0)&(Theta(2)<1)&(Theta(3)<1)
      cont=cont+1;
      T(:,i)=Theta;     
  else 
      T(:,i)=T(:,i-1);
  end
      
  elseif v>min(logratio,0)
      T(:,i)=T(:,i-1);
  end
    ha(1)=T(1,i);
    ha(2)=T(1,i)+T(2,i)*y(1)^2;
 
    for t=3:n
        ha(t)=T(1,i)+T(2,i)*y(t-1)^2+T(3,i)*y(t-2)^2;
    end   
end

tasa = cont/iter;          % Acceptance rate

% The burn-in-period is not used for the bayesian inference. 
% The amount of burn-in-period (7000 in this case) might be
% changed
  TF1=T(1,10000:iter);      
  TF2=T(2,10000:iter);
  TF3=T(3,10000:iter);
  z=length(TF1);
  
% Results
  mean0=mean(TF1);                          % Alpha 0 posterior mean
  mean1=mean(TF2);                          % Alpha 1 posterior mean
  mean2=mean(TF3);                          % Alpha 2 posterior mean
  desv0 = std(TF1,0);                       % Alpha 0 posterior STD
  desv1 = std(TF2,0);                       % Alpha 1 posterior STD
  desv2 = std(TF3,0);                       % Alpha 2 posterior STD
  alpha0real=(alpha0);
  alpha1real=(alpha1);
  alpha2real=(alpha2);
  
  par_lik_bay=[alpha0real,mean0,desv0;alpha1real,mean1,desv1;alpha2real,mean2,desv2];
   
% Volatility estimation using likelihood estimation
      hver(1)=alpha0;
      hver(2)=alpha0+alpha1*y(1)^2;
      for t=3:n
          hver(t)=alpha0+alpha1*y(t-1)^2+alpha2*y(t-2)^2;
      end
% Volatility estimation using mean posterior  
     hpost(1)=mean0;
     hpost(2)=mean0+mean1*y(1)^2;
     for t=3:n
         hpost(t)=mean0+mean1*y(t-1)^2+mean2*y(t-2)^2;
     end
    
% Credibility intervals at 95 %
 percup0=prctile(TF1,97.5);
 percup1=prctile(TF2,97.5);
 percup2=prctile(TF3,97.5);
 percdown0=prctile(TF1,2.5);
 percdown1=prctile(TF2,2.5);
 percdown2=prctile(TF3,2.5);
     perchup(1)= percup0;
     perchdown(1)=percdown0;
     perchup(2)=percup0+percup1*y(1)^2;
     perchdown(2)=percdown0+percdown1*y(1)^2;
     for q=3:n
         perchup(q)=percup0+percup1*y(q-1)^2+percup2*y(q-2)^2;
         perchdown(q)=percdown0+percdown1*y(q-1)^2+percdown2*y(q-2)^2;
     end
 
% Full posterior distribution sampling for the volatility
      for j=1:z
          vol(j)=TF1(j)+TF2(j)*y(n-1)^2+TF3(j)*y(n-2)^2;
      end

% The parameters for the Black & Scholes model are defined
% Of course you should change this parameter considering your own data

 S0=5531;      		%   Actual Mexican Stock Market Index
 X=5550;           	%   Strike Price
 r=.10;           	%   Free risk rate 
 T=.04166666;           %   Duration period
  
% The posterior volatility is implemented in B&S model.
 
 for a=1:z
      d1(a)=(log(S0/X)+(r+(vol(a)^2/2)*T))/(vol(a)*sqrt(T));
      d2(a)=d1(a)-vol(a)*sqrt(T);
      call(a)=S0.*normcdf(d1(a))-X*exp(-r*T)*normcdf(d2(a));
      put(a)=X*exp(-r*T)*normcdf(-d2(a))-S0.*normcdf(-d1(a));
 end
  
% Nonparametric density estimation 
 
 [fa,xa] = ksdensity(TF1);          % Kernel density for alpha 0
 [ga,ya] = ksdensity(TF2);          % Kernel density for alpha 1
 [ha,za] = ksdensity(TF3);          % Kernel density for alpha 2
 [ff,xx] = ksdensity(vol);   % Kernel density for the volatility
 [hh,zz] = ksdensity(call);  %Kernel density for the call option 
 [ii,ss] = ksdensity(put);   %Kernel density for the put option 

% Output 
 Results.rate=tasa;
 Results.est=par_lik_bay;
 
% Figures
% Main figure
  figure(1)
  subplot(3,2,1)
  plot(TF1)
  title('Iteration for alpha0')
  subplot(3,2,2)
  plot(xa,fa)
  title('alpha0 density')
  subplot(3,2,3)
  plot(TF2)
  title('Iteration for alpha1')
  subplot(3,2,4)
  plot(ya,ga)
  title('alpha1 density')
  subplot(3,2,5)
  plot(TF3)
  title('Iteration for alpha1')
  subplot(3,2,6)
  plot(za,ha)
  title('alpha2 density')

% Comparison between volatility estimation using likelihood posterior mean
% and credibility intervals
  figure(2)
  plot(hpost,'b--')
  hold on
  plot(hver,'k-.')
  hold on
  plot(perchup,'r:')
  hold on
  plot(perchdown,'r:')

% Time series, returns and its volatility
  figure(3)
  subplot(3,1,1)
  plot(ipc)
  title('Financial time series')
  subplot(3,1,2)
  plot(y)
  title('Returns')
  subplot(3,1,3)
  plot(hver)
  title('Volatility')
    
% Volatility posterior distribution 
  figure(4)
  plot(xx,ff)
  title('Posterior volatility density (Bayesian ARCH(2)model)'); 
  
