% Monte Carlo simulation to compare the skews for call options, from the following three models:
% (1) Black-Scholes
% (2) Jump
% (3) Heston
function skew_comp(M)

T = 1;           % expiration of the call
r = 0.01;        % risk-free short rate
strikes = [45 50 55];
sigma_imp = zeros(3,3);

global N dt

N = 125;         % number of steps
dt = T/N;        % time step

sigma = 0.25;    % volatility for the Black-Scholes model

lambda = 5;      % jump intensity

theta = 0.1;     % the following is a list of parameters used in Heston's model
kappa = 1.3;
xi = 0.25;
vt0 = 0.25;
rho = -0.75;

S0 = 50;         % initial stock price

for k=1:3
sigma_imp(1,k) = impvol(bs_mc(M,S0,strikes(k),T,sigma,r),S0,strikes(k),r,T);
sigma_imp(2,k) = impvol(jump_mc(M,S0,strikes(k),T,lambda,r),S0,strikes(k),r,T);
sigma_imp(3,k) = impvol(heston_mc(M,S0,strikes(k),T,theta,kappa,xi,vt0,rho,r),S0,strikes(k),r,T);
end

plot(strikes/S0,sigma_imp(1,:),'*',strikes/S0,sigma_imp(2,:),'+',strikes/S0,sigma_imp(3,:),'.');
title('Volatility Skews');
axis([0.8 1.2 0 1]);
xlabel('Moneyness'); ylabel('Implied Vol');
legend('Black-Scholes','Jump','Heston');

% Monte Carlo simulation to price a call in Black-Scholes model:
function price = bs_mc(M,S0,strike,T,sigma,r)
global N dt

stock = zeros(N+1,1);
payoff = zeros(M,1);

for k=1:M

stock(1) = S0;
logretd = log(1 + r*dt + sqrt(dt)*sigma*randn(N,1));     % generate normally distributed random numbers
stock(2:N+1) = S0*exp(cumsum(logretd));

payoff(k) = max(stock(N+1)-strike,0);

end

price = exp(-r*T)*mean(payoff);

% Monte Carlo simulation to price a call in the jump model:
function price = jump_mc(M,S0,strike,T,lambda,r)
global N dt

Nj_max = 50;

payoff = zeros(M,1);

for k=1:M

tau = -log(1-rand(Nj_max,1))/lambda;             % simulate the time periods in between jumps
jumptime = cumsum(tau);
njump = length(find(jumptime < T));              % jump times
tjump = jumptime(1:njump);
jjump = 1 + floor(tjump/dt+0.5);                 % jump time indices j
newvol = rand(njump,1);

stock = S0*ones(N+1,1);
vol = 0.5*ones(N+1,1);

for jp=1:njump
vol(jjump(jp):N+1) = newvol(jp);
end

logretd = log(1 + r*dt + sqrt(dt)*vol(1:N).*randn(N,1));
stock(2:N+1) = S0*exp(cumsum(logretd));

payoff(k) = max(stock(N+1)-strike,0);

end

price = exp(-r*T)*mean(payoff);


% Monte Carlo simulation to price a call in Heston model:
function price = heston_mc(M,S0,strike,T,theta,kappa,xi,vt0,rho,r)
global N dt

stock = zeros(N+1,1);
payoff = zeros(M,1);
vt = zeros(N+1,1);

for k=1:M

w1 = zeros(N,1);
vt(1) = vt0;

for j=1:N
vtry = -1;
while vtry < 0
wz = randn(2,1);
wz(2) = rho*wz(1) + sqrt(1-rho^2)*wz(2);
vtry = vt(j) + kappa*(theta-vt(j))*dt + xi*sqrt(vt(j)*dt)*wz(2);
end
w1(j) = wz(1);
vt(j+1) = vtry;

end

logretd = log(1 + r*dt + sqrt(dt*vt(1:N)).*w1);
stock(2:N+1) = S0*exp(cumsum(logretd));

payoff(k) = max(stock(N+1)-strike,0);

end

price = exp(-r*T)*mean(payoff);
