%
% Monte Carlo Euler method for the approximation
% of option prices depending on one asset.
%  Numerical SDE course at FSU, Fall, 2006.
% Data definition
clear all
clc
close all
rho = -0.3;   % correlation coefficient
r   = 0.04;  % short rate
alpha = 2e2; % mean reversion coefficient
% Initial values
S0    = 100;
Y0    = -1;
V0    = -1;
% Payoff definition
g = inline('max(x-100,0)','x'); % payoff of final contract.
T = 0.75; % expiry time of the option
% Run Parameters
M = 1e4; % Number of realizations.
Co = 1.65; % Confidence interval coeficient.
N = 1e4; % number of time steps;
NN = floor(N/2);
%--------------------------
timev = linspace(0,T,N+1);
gvals = zeros(M,1);
gvals2 = zeros(M,1);
gvals1 = zeros(M,1);
gvals3 = zeros(M,1);
gvals3d = zeros(M,1);
gvals3V= zeros(M,1);
dt = T/N; %step size
dtt= 2*dt;
S = ones(N+1,1)*S0;
S2 = ones(N+1,1)*S0;
S1 = ones(N+1,1)*S0;
S3 = ones(NN+1,1)*S0;
Y = ones(N+1,1)*Y0;
Yp= ones(NN+1,1)*Y0;
V = ones(N+1,1)*V0;
tmp = 0.4*sqrt(alpha);
for m = 1:M
  dW = sqrt(dt) * randn(N,1);
  dZ = sqrt(dt) * randn(N,1);
  dhZ = rho*dW+sqrt(1-rho^2)*dZ; 
  for n = 1:N
    S2(n+1)= (1+ dt*r + dW(n)*exp(-1)  )*S2(n);
    S(n+1) = (1+ dt*r + dW(n)*exp(Y(n)))*S(n);
    S1(n+1) = (1+ dt*r + dW(n)*exp(V(n)))*S1(n);
    Y(n+1) = (Y(n) + (-alpha+tmp*sqrt(1-rho^2))*dt + tmp * dhZ(n))/(1+alpha*dt);
    V(n+1) = V(n)+(-alpha*V(n)-alpha+tmp*sqrt(1-rho^2))*dt + tmp * dhZ(n);
  end  
  gvals(m) = exp(-r*T)*g(S(n+1));  
  gvals2(m) = exp(-r*T)*g(S2(n+1));
  gvals1(m) = exp(-r*T)*g(S1(n+1));
  
  for n=1:NN
      nn=2*n;
      dWW=(dW(nn)+dW(nn-1))*sqrt(1);
      dZZ=(dZ(nn)+dZ(nn-1))*sqrt(1);
      dhZZ=rho*dWW+sqrt(1-rho^2)*dZZ;
      S3(n+1)= (1+ dtt*r + dWW*exp(Yp(n)))*S3(n);
      Yp(n+1)=(Yp(n) + (-alpha+tmp*sqrt(1-rho^2))*dtt + tmp * dhZZ)/(1+alpha*dtt);
  end
  gvals3(m) = exp(-r*T)*g(S3(n+1));
  gvals3d(m) = abs (gvals3(m) - gvals(m));
  gvals3V(m) = gvals3d(m)^2;
  if mod(m,1000)==0
     m 
  end

end
figure, subplot(511),
plot(timev,S),hold on, plot(timev,S2,'r'), hold on, plot(timev,S1,'g'),
legend('Stochastic Volatility Implicit','Black-Scholes','Stochastic Volatility Explicit',-1);
title 'S pathwise approximations'
subplot(512)
plot(timev,exp(Y),'r',timev,exp(V),'g')
ylabel '\sigma(Y) = e^Y'
legend('Stochastic Volatility Implicit','Stochastic Volatility Explicit',-1);
subplot(513)
plot(timev,Y,'r',timev,V,'g')
xlabel 'time',ylabel 'Y'
legend('Approx. Implicit Y Process','Approx. Explicit Y Process',-1);
subplot(514)
plot(timev,Y,'r')
xlabel 'time',ylabel 'Y'
legend('Approx. Implicit Y Process',-1);
subplot(515)
plot(timev,V,'g')
xlabel 'time',ylabel 'Y'
legend('Approx. Explicit Y Process',-1);

figure
subplot(311)
hist(gvals)
title 'Sample histogram, Implicit Stochastic volatility model'
subplot(312)
hist(gvals2)
title 'Sample histogram, Black-Scholes model'
subplot(313)
hist(gvals1)
title 'Sample histogram, Explicit Stochastic volatility model'

implicitbs=mean(gvals) 
bs=mean(gvals2)
explicitbs=mean(gvals1)
implicitbs_big_delta= mean(gvals3)
mean_error=mean(gvals3d)
desv_dt=sqrt(mean(gvals3V) - mean_error^2);
error_dt=Co* desv_dt / sqrt(M) 
[Call, Put] = blsprice(S0, 100, r, T, exp(-1));
Call

