%  options_Merton.m 
%  Option prices and implied vols for speical case of Merton model (poisson-normal mixture)  
%  Merton (JFE, 1976) is the reference to the general setup 
%  Notation 
%    log re = log r1 + N(mu,sigma^2) + j N(theta,delta^2)
%    j = # poisson jumps, intensity omega 
%  Backus, Chernov, and Martin, "Equity index options"
%  Written:  April 2009 and after 
format compact 
clear all
close all

disp(' ')
disp('Merton option model') 
disp('---------------------------------------------------------------')
disp(' ')
%
%  1. parameter inputs 
disp('Parameter values') 
% data 
logr1 = 0.025;
q1 = exp(-logr1);
ep = 0.0400;
sigmare = 0.1800;

% true distribution (sigma is all we need here) 
omega = 1.5120;
theta = -0.0259;
delta = 0.0407;
mu = ep - omega*theta
sigma = sqrt(sigmare^2-omega*(theta^2+delta^2))

% risk-neutral distribution 
disp('omega* theta* delta*')
omegas = 1.5120; deltas = 0.0981; thetas = log(1-0.0482)-deltas^2/2;    %  Chernov
[omegas thetas deltas] 

% set mu* to satisfy arb condition 
mus = - sigma^2/2 - omegas*(exp(thetas+deltas^2/2)-1)
%mudiff = mu - mus

%  2. Option prices (Merton model) 
disp(' ')
disp('Option prices') 

% strike grid 
ngrid = 10;
bscale = 0.06;
%logbgrid = bscale*[-ngrid:ngrid]'/ngrid;
%bgrid = exp(logbgrid);
bgrid = bscale*[-ngrid:ngrid]'/ngrid + 1;
logbgrid = log(bgrid);

% probs 
nprob = 10;
jgrid = [0:nprob]';
probs = exp(-omegas)*omegas.^jgrid./gamma(jgrid+1);
checksumprobs = sum(probs)  %  should be exactly 1.0 

% option prices (compute for given number of jumps j, then average using Poisson probs) 
% based on put prices (despite what the variable names suggest) 
% vectorize and skip loop?  
callj = zeros(size(bgrid));
qcall = callj;
dgrid = zeros(size(bgrid));
for it = 1:length(probs) 
    j = jgrid(it);
    probj = probs(it);
    kappa1j  = mus + j*thetas;
    kappa2j = sigma^2 + j*deltas^2;
    dgrid = (logbgrid - logr1 - kappa1j)/sqrt(kappa2j);
    callj = q1*bgrid.*normcdf(dgrid) - exp(kappa1j+kappa2j/2)*normcdf(dgrid-sqrt(kappa2j));
    qcall = qcall + probj*callj;
end     

%disp('Strike and put price') 
%[bgrid qcall]

%return 

%  3. Implied vols  
disp(' ')
disp('Implied volatilities') 
vol = zeros(size(qcall)) + sigmare;
kappa1 = zeros(size(qcall));
kappa2 = kappa1;
maxit = 1000;
tol = 1.e-10;    % tolerance set on relative price error  
put = qcall;     % true put prices 

% Newton's method on f(log(vol)) = bs_price - price
% derivative fp is the "vega" divided by vol 
for it=1:maxit
    it;
    kappa2 = vol.^2;
    kappa1 = -kappa2/2;
    dgrid = (logbgrid - logr1 - kappa1)./vol;
    putbsmguess = q1*bgrid.*normcdf(dgrid) - normcdf(dgrid-vol);
    f = putbsmguess - put;
    fp = normpdf(dgrid-vol).*vol;      %  NB:  this is the pdf 
    if max(abs(f./put)) < tol, break, it, end
    logvolnew = log(vol) - f./fp;
    volnew = exp(logvolnew);
    vol = volnew;
end  

disp('Convergence indicators for vols')
it 
maxerror = max(f)
maxrelerror = max(f./put)
volmodel = vol;

disp('Strike, put price, implied vol') 
[bgrid qcall vol]

plot(bgrid,vol)

return