# Heston Model Calibration with MatLab: model prices do not fall in the bid-ask range

I am currently implementing the MatLab code reported below for the calibration of Heston Model. The code seems fine and, by reading the paper where I took the code, I was able to calibrate and price obtaining exactly the same results as the paper.

Then I have tried to use the same code for call options on SP500 but I could not obtain the same results as it was from the data of the paper. In particular, when I try to price the model prices always fall outside the bid-ask range.

My question then is why this happens? I guess something must be wrong with the data (attached find both paper data and my SP500 data)

function cost = costf(x)
global data; global finalcost;

% Compute individual differences
% Sum of squares performed by Matlab's lsqnonlin
for i=1:length(data)
cost(i)= data(i,5) - call_heston_cf(data(i,1),x(1), x(2), (x(5)+x(3)^2)/(2*x(2)), x(3), data(i,4), x(4), data(i, 2), data(i,3));
end
% Show final cost
finalcost=sum(cost)^2;
end

function y = chfun_heston(s0, v0, vbar, a, vvol, r, rho, t, w)
% Heston characteristic function.
% Inputs:
% s0: stock price
% v0: initial volatility (v0^2 initial variance)
% vbar: long-term variance mean
% a: variance mean-reversion speed
% vvol: volatility of the variance process
% r : risk-free rate
% rho: correlation between the Weiner processes for the stock price and its variance
% w: points at which to evaluate the function
% Output:
% Characteristic function of log (St) in the Heston model
% Interim calculations
alpha = -w.*w/2 - i*w/2;
beta = a - rho*vvol*i*w;
gamma = vvol*vvol/2;
h = sqrt(beta.*beta - 4*alpha*gamma); rplus = (beta + h)/vvol/vvol;
rminus = (beta - h)/vvol/vvol; g=rminus./rplus;
% Required inputs for the characteristic function
C = a * (rminus * t - (2 / vvol^2) .* log((1 - g .* exp(-h*t))./(1-g)));
D = rminus .* (1 - exp(-h * t))./(1 - g .* exp(-h*t));

% Characteristic function evaluated at points w
y = exp(C*vbar + D*v0 + i*w*log(s0*exp(r*t)));

function y = call_heston_cf(s0, v0, vbar, a, vvol, r, rho, t, k)
% Heston call value using characteristic functions.
% y = call_heston_cf(s0, v0, vbar, a, vvol, r, rho, t, k)
% Inputs:
% s0: stock price
% v0: initial volatility (v0^2 initial variance)
% vbar: long-term variance mean
% a: variance mean-reversion speed
% vvol: volatility of the variance process
% r: risk-free rate
% rho: correlation between the Weiner processes of the stock price and its variance
% t: time to maturity
% k: option strike
% chfun_heston: Heston characteristic function
% 1st step: calculate pi1 and pi2

% Inner integral 1
int1 = @(w, s0, v0, vbar, a, vvol, r, rho, t, k) real(exp(-i.*w*log(k)).*chfun_heston(s0, v0, vbar, a, vvol, r, rho, t, w-i)./(i*w.*chfun_heston(s0, v0, vbar, a, vvol, r, rho, t, -i)));% inner integral1
int1 = integral(@(w)int1(w,s0, v0, vbar, a, vvol, r, rho, t, k),0,100); % numerical integration
pi1 = int1/pi+0.5; % final pi1

% Inner integral 2:
int2 = @(w, s0, v0, vbar, a, vvol, r, rho, t, k) real(exp(-i.*w*log(k)).*chfun_heston(s0, v0, vbar, a, vvol, r, rho, t, w)./(i*w));
int2 = integral(@(w)int2(w,s0, v0, vbar, a, vvol, r, rho, t, k),0,100);int2 = real(int2);
pi2 = int2/pi+0.5; % final pi2
% 2rd step: calculate call value
y = s0*pi1-exp(-r*t)*k*pi2;
end

% Heston calibration, local optimization (Matlab's lsqnonlin)
% Input on data.txt
% Data = [So, t, k, r, mid price, bid, ask]
clear all;clc;
global data; global cost; global finalcost;
% Initial parameters and parameter bounds
% Bounds [v0, Vbar, vvol, rho, 2*a*vbar - vvol^2]
% Last bound include non-negativity constraint and bounds for mean-reversion
x0 = [.5,.5,1,-0.5,1];
lb = [0, 0, 0, -1, 0];
ub = [2, 2, 5, 1, 20];
% Optimization: calls function costf.m:
tic;
x = lsqnonlin(@costf,x0,lb,ub);
toc;
% Solution:
Heston_sol = [x(1), x(2), x(3), x(4), (x(5)+x(3)^2)/(2*x(2))]
x
min = finalcost


UPDATE: I have tried on new data taken from another thesis, here attached, and the code just got 7 out of 24 prices correctly. It seems that somehow the code does not work properly (the code is actually on a publication). Does anyone have any idea on where the problem comes from? Thanks again! • Hi, could you please add the result columns to the tables as well? Or maybe a plot of the data and your results? Thanks! Feb 7 at 5:30
• Another question: shouldn’t the goal read sum(c^2) instead of sum(c)^2 ? Feb 7 at 10:53
• Yes, you right, it should be c^2. However I obtain very small improvement on the parameters and overall. Feb 7 at 15:56
• The parameters are now v0=0.0395, vbar=0.0538, vvol=0.8733, rho=-0.9539 and a=7.0832. Let's say now I want to check the calibration, so I try to price the call with maturity= 0.1972603, strike= 3800 and rate=0.0433%. Then I input in the Matlab command window the following :call_heston_cf(3830.17, 0.0395, 0.0538, 7.0832, 0.8731, 0.000433, -0.9539, 0.1972603, 3800). But Ans= 153.8467 and as you can see from the table is outside the bid-ask range. Feb 7 at 16:04
• In the table I uploaded green is when model prices are in the bid - ask range, yellow when they are higher and blue when they are lower than the range. Hope this helps. Feb 7 at 16:06