Diagnosing Heston Calibration Errors Against Option Bid–Ask Prices

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


# 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;
[data] = xlsread('data.xlsx');
% 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!