Feeds:
Posts

## MATLAB CODE FOR HESTON MODEL

This is from the 2012 thesis-fastcalibration-heston-model, and extremely helpful for us.

MATLAB Code
B.1
Calibration
function [x] = run ()
%Initial Parameter Set
%x0 = [kappa, theta, sigma, rho, v0]
x0 = [0.2542 0.9998 0.8096 -0.0001 0.0119]
%Lower and upper bound for optimisation
lb = [0 0 0 -1 0];
ub = [20 1 5 0 1];
datensatz = xlsread( ‘ Satensatz.xlsx ‘ , ‘ all ‘ , ‘ A2:I251 ‘ );
% Maturity
T = datensatz(:,2);
% Strike price
strike = datensatz(:,4);
% Underlying price
Under = datensatz(:,3);
%Implied Volatility
implvol = datensatz(:,9);
tic; % start timer
% Optimisation
options = optimset( ‘ TolFun ‘ , 1e-8);
101102
APPENDIX B.
MATLAB CODE
[x,resnorm,unused,exitflag]
= lsqnonlin(@(x) costf3(x,T,strike,Under,implvol),x0,lb,ub,options);
%value of objective function
fprintf( ‘ resnorm = %.20f\n ‘ , resnorm);
%3: Change in the residual was less than the specified tolerance
fprintf( ‘ exitflag = %d\n\n ‘ , exitflag);
fprintf( ‘ kappa=%.20f\n theta=%.20f\n sigma=%.20f\n ‘ ,x(1),x(2),x(3));
fprintf( ‘ rho=%.20f\n v0= %.20f\n ‘ ,x(4),x(5));
fprintf( ‘ \nRechenzeit: %d s\n ‘ , int16(toc)); % stop timer and output
end
function [L] = costf3(x, T, strike, Under, implvol)
for i = 1 : length(T)
%L(i) depends on the regime one wants to use:
L(i)=((implvol(i)^2-SmallTime(x(1),x(2),x(3),x(4),x(5),money(i),T(i)))/implvol(i)^2);
((implvol(i)^2-SmallTimeLDT(x(1),x(2),x(3),x(4),x(5),money(i),T(i)))/implvol(i)^2);
%Large-time, Large-strike:
((implvol(i)^2-LargeTime(x(1),x(2),x(3),x(4),x(5),money(i),T(i)))/implvol(i)^2);
%Small-strike:
(implvol(i)-SmallStrike(x(1),x(2),x(3),x(4),x(5),T(i),strike(i),Under(i)))/implvol(i);
%Large-strike:
(implvol(i)-LargeStrike(x(1),x(2),x(3),x(4),x(5),T(i),strike(i),Under(i)))/implvol(i);
end
endB.2.
B.2
SMALL-MATURITY
Small-maturity
function H = SmallTime(kappa,theta,sigma,rho,v0,money,T)
x = log(money);
alpha = kappa*theta;
term1 = 1 + (rho/2)*(sigma*x/v0) + (1-(7/4)*rho^2)* (sigma^2*x^2)/(12*v0^2);
term2a = (rho*sigma*v0)/2 – (sigma^2/6)*(1-rho^2/4) + alpha – kappa*v0;
term2b = sigma^2*(1-rho^2) + rho*sigma*v0 – 2*alpha – 2*kappa*v0;
term3a = (176 – 712*rho^2 + 521*rho^4)*sigma^2 + 40*sigma*rho^3*v0;
term3b = 80*(13*rho^2 – 6)*alpha – 80*kappa*rho^2*v0;
H1 = v0 * term1 + (term2a + (rho*sigma)/(12*v0) * term2b * x) * (T/2);
H2 = (sigma^2/(7680*v0^2)) * (term3a + term3b) * x^2*T;
H = H1 + H2;
end
B.3
function H = SmallTimeLDT(kappa,theta,sigma,rho,v0,money,T)
x = log(money);
alpha = kappa*theta;
term1 = 1 + (rho*sigma*x)/(4*v0);
term2 = (1/24)*(1 – 5*rho^2/2)*(sigma^2*x^2/v0^2);
H = v0*(term1 + term2)^2;
end
B.4
Large-maturity
function H = LargeTime(kappa,theta,sigma,rho,v0,money,T)
103104
APPENDIX B.
MATLAB CODE
x = (1/T)*log(money);
kappa2 = kappa – rho*sigma;
theta2 = (kappa*theta)/kappa2;
%p corresponds to p*
pterm1 = (sigma-2*kappa*rho)/(2*(1-rho^2)*sigma);
pterm2a = (sigma^2+4*kappa^2 – 4*kappa*rho*sigma);
pterm2b = (x^2*sigma^2+2*x*kappa*theta*rho*sigma+kappa^2*theta^2);
pterm2 = sqrt(pterm2a/pterm2b);
p = pterm1 + (kappa*theta*rho + x*sigma)/(2*(1 – rho^2)*sigma) * pterm2;
Vterm = (kappa – rho*sigma*p – sqrt((kappa-rho*sigma*p)^2 + sigma^2*(p-p^2)));
V = (kappa*theta/sigma^2) * Vterm;
Vstern = p*x – V;
if ((x > theta2/2) || (x < -theta/2) )
H = 2*(2*Vstern – x – 2*sqrt(Vstern^2 – Vstern*x));
else
H = 2*(2*Vstern – x + 2*sqrt(Vstern^2 – Vstern*x));
end
end
B.5
Small-strike
function [H] = SmallStrike(kappa, theta, sigma, rho, v0, T, strike, Under)
a
b
c
k
=
=
=
=
kappa*theta;
-kappa;
sigma;
log(strike/Under);
%root finding procedure
NullstelleGefunden = false;
Start = 100;
options = optimset( ‘ Display ‘ , ‘ off ‘ );
while ~NullstelleGefundenB.5.
SMALL-STRIKE
105
fs = ExplosionTime(Startwert, rho, c, b, T);
if isfinite(fs) && isreal(fs)
[solution, fval, exitflag] = fzero(@(s) ExplosionTime(s, rho, c, b, T), Start, options);
if (exitflag == 1)
if isreal(solution) & solution >= -1
NullstelleGefunden = true;
s = solution;
end;
end;
end;
Start = Start + 100;
if (Start >= 5000)
fprintf( ‘ Startwert schon bei %d\n ‘ , Start);
disp(rho);
disp(c);
disp(b);
disp(T);
end;
end
B3 = -(s+1);
B2 = 2* (2*v0)^(1/2)/(c * slope^(1/2));
rho1= 2^(1/2) * ((B3+2)^(1/2) – (B3+1)^(1/2));
rho2= B2*2^(-1/2) * ( 1/(B3+1)^(1/2) – 1/(B3+2)^(1/2));
rho3= 2^(-1/2) * (1/4 – a/c^2) * ( 1/(B3+2)^(1/2) – 1/(B3+1)^(1/2));
H = (1/T^(1/2)) * (rho1*(-k)^(1/2) + rho2 + (rho3 * log(-k))/((-k)^(1/2)));
end
function Tstern = ExplosionTime(s, rho, c, b, T)
term = (atan(sqrt(-(s*rho*c+b)^2 + c^2*(s^2-s))/(s*rho*c+b)));
Tstern = ((2/sqrt(-(s*rho*c+b)^2 + c^2*(s^2-s))) * term)-T;
end106
APPENDIX B.
B.6
MATLAB CODE
Large-strike
function [H] = LargeStrike(kappa, theta, sigma, rho, v0, T, strike, Under)
a
b
c
k
=
=
=
=
kappa*theta;
-kappa;
sigma;
log(strike/Under);
%root finding procedure
NullstelleGefunden = false;
Start = 100;
options = optimset( ‘ Display ‘ , ‘ off ‘ );
while ~NullstelleGefunden
fs = ExplosionTime(Start, rho, c, b, T);
if isfinite(fs) && isreal(fs)
[solution, fval, exitflag] = fzero(@(s) ExplosionTime(s, rho, c, b, T), Start, options);
if (exitflag == 1)
if isreal(solution) & solution >= -1
NullstelleGefunden = true;
s = solution;
end;
end;
end;
Start = Start + 100;
if (Startwert >= 5000)
fprintf( ‘ Startwert schon bei %d\n ‘ , Start);
disp(rho);
disp(c);
disp(b);
disp(T);
end;
end
term1a = T*c^2*s*(s-1) * (c^2*(2*s-1) – 2*rho*c*(s*rho*c + b));
term1b = 2*(s*rho*c +b) * (c^2*(2*s-1) – 2*rho*c*(s*rho*c + b));
term1c = 4*rho*c * (c^2*s*(s-1) – (s*rho*c + b)^2);
R1 = term1a – term1b + term1c;B.7.
HESTON MODEL
107
R2 = 2*c^2*s * (s-1) * (c^2*s*(s-1) – (s*rho*c + b)^2);
%slope corresponds to sigma in formula
slope = R1 / R2;
A2 = 2*(sqrt(2*v0))/(c*sqrt(slope));
A3 = s + 1;
beta1 = sqrt(2)*(sqrt(A3-1) – sqrt(A3-2));
beta2 = (A2/sqrt(2)) * (1/sqrt(A3-2) – 1/sqrt(A3-1));
beta3 = (1/sqrt(2)) * (1/4 – a/c^2) * (1/sqrt(A3-1) – 1/sqrt(A3-2));
H = (1/sqrt(T)) * (beta1*k^(1/2) + beta2 + (beta3 * log(k))/k^(1/2));
end
function Tstern = ExplosionTime(s, rho, c, b, T)
term = (atan(sqrt(-(s*rho*c+b)^2 + c^2*(s^2-s))/(s*rho*c+b)) + pi);
Tstern = ((2/sqrt(-(s*rho*c+b)^2 + c^2*(s^2-s))) * term)-T;
end
B.7
Heston Model
The MATLAB code for the semi-closed form solution of the Heston model is based on [51].
warning off;
if(opt==1)
call = s0*HestonP(kappa,theta,sigma,rho,v0,r,T,s0,K,1)
– K*exp(-r*T)*HestonP(kappa,theta,sigma,rho,v0,r,T,s0,K,2);
else
call= s0*HestonP(kappa,theta,sigma,rho,v0,r,T,s0,K,1)
– K*exp(-r*T)*HestonP(kappa,theta,sigma,rho,v0,r,T,s0,K,2)-s0+K*exp(-r*T);
end
function ret = HestonP(kappa,theta,sigma,rho,v0,r,T,s0,K,type)
APPENDIX B.
MATLAB CODE
function ret = HestonPIntegrand(phi,kappa,theta,sigma,rho,v0,r,T,s0,K,type)
ret = real(exp(-i*phi*log(K)).*Hestf(phi,kappa,theta,sigma,rho,v0,r,T,s0,type)./(i*phi));
function f = Hestf(phi,kappa,theta,sigma,rho,v0,r,T,s0,type);
if type == 1
u = 0.5;
b = kappa – rho*sigma;
else
u = -0.5;
b = kappa;
end
a = kappa*theta; x = log(s0);
d = sqrt((rho*sigma*phi.*i-b).^2-sigma^2*(2*u*phi.*i-phi.^2));
g = (b-rho*sigma*phi*i + d)./(b-rho*sigma*phi*i – d);
C = r*phi.*i*T + a/sigma^2.*((b- rho*sigma*phi*i + d)*T -2*log((1-g.*exp(d*T))./(1-g)));
D = (b-rho*sigma*phi*i + d)./sigma^2.*((1-exp(d*T))./(1-g.*exp(d*T)));
f = exp(C + D*v0 + i*phi*x);