I am trying to calibrate Heston model on some data and I have the following code. Code is supposed, after it reads the data, to give back 5 parameters. However, I get an empty answer from MatLab. Does somebody know what's wrong with this code? Thanks in advance.
function call = HestonCallQuadl(kappa,theta,sigma,rho,v0,r,T,s0,K)
warning off;
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);
function ret = HestonP(kappa,theta,sigma,rho,v0,r,T,s0,K,type)
ret = 0.5 + 1/pi*quadl(@HestonPIntegrand,0,100,[],[],kappa,theta,sigma,rho,v0,r,T,s0,K,type);
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);
function ret = HestonCallDifferences(input)
global NoOfIterations; global NoOfOptions; global PriceDifference;
NoOfIterations = NoOfIterations + 1;
%counts the no of iterations run to calibrate model
for i = 1:NoOfOptions
PriceDifference(i) = (OptionData(i,5)-HestonCallQuadl((input(1)+input(3)^2)/(2*input(2)),input(2),input(3),input(4),input(5),OptionData(i,1),OptionData(i,2),OptionData(i,3),OptionData(i,4)))/sqrt((abs(OptionData(i,6)-OptionData(i,7))));
%input matrix = [kappa theta sigma rho v0]
end
ret = PriceDifference';
end
clear;
global data; global OptionData;
NoOfIterations = 0;
[data,text]=xlsread('OptionData.xlsx');
%OptionData = format [r-q,T,S0,K,Option Value,bid,offer]
Size = size(data);
NoOfOptions = Size(1);
OptionData = data;
%input sequence in initial vectors [2*kappa*theta - sigma^2,theta,sigma,rho,v0]
x0 = [1 0.5 1 -0.5 0.5];
lb = [0 0 0 -1 0];
ub = [20 1 5 0 1];
options = optimset('MaxFunEvals',20000);
%sets the max no. of iteration to 20000 so that termination doesn't take place early.
tic;
Calibration = lsqnonlin(@HestonCallDifferences,x0,lb,ub);
toc;
Solution = [(Calibration(1)+Calibration(3)^2)/(2*Calibration(2)), Calibration(2:5)];