function CMC05ATOCV
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
% Using log-uniform jump-diffusion model: DS(t)/S(t-) = mu dt + sigma dw + J(Q)dP;  
% The Option pricing method is used by risk-neutral pricing method; 
% The Numerical Method is Conditional Monte Carlo Method  with Antithetic and 
% Optimal Control Variates Techniques;  
% For details, please see Zongwu Zhu and Floyd B. Hanson's working paper: 
% <<Option Pricing for Log-Uniform Jump-Diffusion Model>>.
% The purpose of this program is for numerical test part of the paper, not for practical use.  
% For practical use, the codes need some modifications.
% Programmed by Zongwu Zhu 
% Date 12/28/2004; 
% Updated: Jan 16, 2005
% Updated: April 22, 2005 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc;
clf % clear figures
fprintf('\nfunction CMC05ATOCV InPut/OutPut: (%s):\n',datestr(now))
% Model parameters 
S0=100; T = 0.2; lam = 100; qa = -0.3; qb = 0.2; sigma = 0.20; K = 105;
r=0.05; %ZZ: changed 8/15/2005
%n = 4e5; % In order to get almost exact value, otherwise n = 1e4 suggested 
n=1e4;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input test parameters: 
R = K/S0; % input('the ration of K/S0 :\n');
% sigma = input('the \sigma :\n'); 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
K = S0*R; 
tic;    % time start;
t1=toc;
% Get the Black-Scholes' Call Price
BSC = BS(S0, K, T, r, sigma); 
%
lamt = lam*T;
meanj = (exp(qb)-exp(qa))/(qb - qa) - 1; 
mean2j = (exp(2*qb)-exp(2*qa))/(2*qb - 2*qa) - 1; 
rand('state',1); % Ignore Depreciated Warning
X = zeros(1,n);
Y = zeros(1,n);
for i =1:n 
       % Simulate NJ: the jumps during the period T
       p = exp(-lamt); 
      Phi = p; 
      NJ = 0; 
      uj = rand(1); 
      while uj > Phi
        NJ = NJ + 1; 
        p = p*lamt/NJ;
        Phi = Phi + p; 
      end 
     % Generate the jump part;
     ZJ = rand(NJ,1); 
     AZJ = 1-ZJ; 
     jump1 = sum(qa+ZJ*(qb-qa));
     jump2 = sum(qa + AZJ*(qb-qa)); 
     % The simulated Black-Scholes' price is: 
     BSMt1 = BS(S0*exp(-lamt*meanj+jump1),K,T,r,sigma);
     BSMt2 = BS(S0*exp(-lamt*meanj+jump2),K,T,r,sigma);
     % Generate antithetic variate X_i and control variate Y_i
     X(i) = 0.5*(BSMt1 + BSMt2); 
     Y(i) = 0.5*(exp(jump1) + exp(jump2)); 
end % end for loop  
mx = mean(X);
my = mean(Y); 
vc= 0.5*(exp(lamt*mean2j) - 2*exp(lamt*2*meanj)+exp(lamt*(exp(qa+qb)-1))); % VC = Var(Y)
alpha =n/(n-1)*(mean(X.*Y)-mx*my)/vc;
BSM = X - alpha*(Y-exp(lamt*meanj));
Cn =mean(BSM);  
% Make a correction 
Z = (2*my -Y).*Y; 
mz = mean(Z);
b = 1/(n-1)*(mean(X.*Z)-mx*mz)/vc; 
Cn = Cn - b; 
% Get the standard error estimation
variance = var(BSM);
stderror =std(BSM)/sqrt(n); 
% Get the put value by put-call parity
Pn = Cn + K*exp(-r*T) - S0; 
t2=toc; 
%
fprintf('The European call option price by Monte Carlo method with AOCV is:%8.4f;\n', Cn); 
fprintf('The European put option price is:%8.4f;\n', Pn); 
fprintf('The Black-Scholes Price is %8.4e\n', BSC);
fprintf('variance of Conditional Black Scholes formula is %8.3f.\n', variance); 
fprintf('The standard error  is %8.4f\n', stderror);
fprintf('The computing time is: %8.3e seconds, the sample size is %i.\n', t2-t1, n);
fprintf('The sqrt of computing time multiply error is: %8.3e.\n', sqrt((t2-t1))*stderror); 
%%%%%
function Phinorm = Phinorm(d)  % the standard accumulated distribution function
Phinorm = 1 - 0.5*erfc(d/sqrt(2)); 
%%%%%
function BS = BS(S0,K,T,r,sigma) % Black-Scholes option pricing formula
d1 = (log(S0/K) + (r + 0.5*sigma^2)*T)/(sigma*sqrt(T)); 
d2 = d1 - sigma*sqrt(T); 
BS = S0*Phinorm(d1)-K*exp(-r*T)*Phinorm(d2);
%%%%%%%%%