Hello,
I'm trying to do the optimization for Heston parameters using the optmum package (I already tested it on Black-Scholes and it works). The problem on Heston is that, after giving the initial values, the optimization procedure runs for ~10 minutes but the values are not optimized. The code is as follows:
new;
library pgraph,optmum;
pqgwin many;
declare matrix _Heston_sigma;
declare matrix _Heston_b;
declare matrix _Heston_a;
declare matrix _Heston_sigmav;
declare matrix _Heston_rho;
load data[]=D:\Disertatie\Data.txt;
data=reshape(data,rows(data)./6,6);
data=delif(data,data[.,6].<=6);
data=delif(data,data[.,2].<=3/8);
data=delif(data,data[.,3]./data[.,4].>=5);
r=data[.,5]/100;
tau=data[.,6]/365;
S=data[.,3];
K=data[.,4];
pret_real=data[.,2];
/* Formula de pricing a unei optiuni europene folosind modelul Heston - Volatilitate Stochastica */
proc HESTON_CALL(&f,&g,param);
local c,P1,P2,lmax,lmin;
local f:proc;
local g:proc;
_Heston_sigma = param[1];
_Heston_b = param[2];
_Heston_a = param[3];
_Heston_sigmav = param[4];
_Heston_rho = param[5];
lmax=1000;
lmin=0;
P1 = 1/2+(1/pi).*intquad1(&f,lmax|lmin);
P2 = 1/2+(1/pi).*intquad1(&g,lmax|lmin);
c=S.*P1-K.*exp(-r.*tau).*P2;
retp(c);
endp;
proc f1(phi);
local sigma,a,b,rho,sigmav;
local i,fct,g,lambda1,lambda2,term1,term2;
sigma = _Heston_sigma;
b = _Heston_b; /*long term mean */
a = _Heston_a; /* mean reversion speed */
sigmav = _Heston_sigmav;
rho = _Heston_rho;
i=sqrt(-1);
lambda1=-(1+i.*phi).*(-0.5+rho.*a./sigmav+0.5.*(1+i.*phi).*(1-rho.^2));
lambda2=-(1+i.*phi).*rho./sigmav;
g=sqrt(a^2+2.*(sigmav.^2).*lambda1);
term1=((2.*g.*exp((a+g).*tau./2))./(sigmav.^2.*lambda2.*(exp(g.*tau)-1)+(g-a)+exp(g.*tau).*(a+g)))^(2.*a.*b./(sigmav.^2));
term2=exp(-(sigma.^2).*(lambda2.*(a+g+exp(g.*tau).*(g-a))+2.*lambda1.*(exp(g.*tau)-1))./(sigmav.^2.*lambda2.*(exp(g.*tau)-1)+(g-a)+exp(g.*tau).*(a+g)));
fct=exp(i.*phi.*(ln(s)+r.*tau)-(rho./sigmav).*(1+i.*phi).*(sigma.^2+a.*b.*tau)).*term1.*term2;
retp(real(fct.*exp(-i.*phi.*ln(k))./(i.*phi)));
endp;
proc f2(phi);
local sigma,a,b,rho,sigmav;
local i,fct,g,lambda1,lambda2,term1,term2;
sigma = _Heston_sigma;
b = _Heston_b;
a = _Heston_a;
sigmav = _Heston_sigmav;
rho = _Heston_rho;
i=sqrt(-1);
lambda1=-(i.*phi).*(-0.5+rho.*a./sigmav+0.5.*(i.*phi).*(1-rho^2));
lambda2=-(i.*phi).*rho./sigmav;
g=sqrt(a^2+2.*(sigmav.^2).*lambda1);
term1=((2.*g.*exp((a+g).*tau./2))./(sigmav.^2.*lambda2.*(exp(g.*tau)-1)+(g-a)+exp(g.*tau).*(a+g)))^(2.*a.*b./(sigmav.^2));
term2=exp(-(sigma.^2).*(lambda2.*(a+g+exp(g.*tau).*(g-a))+2.*lambda1.*(exp(g.*tau)-1))./(sigmav.^2.*lambda2.*(exp(g.*tau)-1)+(g-a)+exp(g.*tau).*(a+g)));
fct=exp(i.*phi.*(ln(s)+r.*tau)-(rho./sigmav).*(i.*phi).*(sigma.^2+a.*b.*tau)).*term1.*term2;
retp(real(fct.*exp(-i.*phi.*ln(k))./(i.*phi)));
endp;
prm =0.1|0.2|0.3|0.4|0.5;
x=HESTON_CALL(&f1,&f2,prm);
proc eroare(&f,&g,&h,param);
local f:proc;
local g:proc;
local h:proc;
param[1]=abs(param[1]);
param[2]=abs(param[2]);
param[3]=abs(param[3]);
param[4]=abs(param[4]);
param[5]=sin(param[5]);
retp(f(&g,&h,param)-pret_real);
endp;
proc eroare_patrat_heston(param);
retp(sumc(eroare(&HESTON_CALL,&f1,&f2,param).^2));
endp;
proc implicit_heston(param);
local x,y,gradient,retcode;
{x,y,gradient,retcode}=optmum(&eroare_patrat_heston,param);
x[1]=abs(x[1]);
x[2]=abs(x[2]);
x[3]=abs(x[3]);
x[4]=abs(x[4]);
x[5]=sin(x[5]);
retp(x);
endp;
xheston=0.18|0.04|1|0.1|-0.64;
para=0;
heston=0;
i=data[1,1]; /* prima zi din esantion */
datas=data;
do while i<=data[rows(data),1];
datam=selif(data,data[.,1].==i);
datas=delif(datas,datas[.,1].==i);
i=datas[1,1];
s=datam[.,3];
k=datam[.,4];
r=datam[.,5]./100;
tau=datam[.,6]./365;
pret_real=datam[.,2];
parame=implicit_heston(xheston);
para=para|parame;
heston=heston|HESTON_CALL(&f1,&f2,parame);
endo;
heston=trim(heston,1,0);
vol=trim(para,1,0);
vol=reshape(vol,rows(vol)./5,5);
volatilitate_implicita=vol[.,1];
LT_mean=vol[.,2];
mean_reversion=vol[.,3];
sigma_vol=vol[.,4];
rho=vol[.,5];
err=heston-data[.,2];
dataf=data~err;
s=dataf[.,3];
k=dataf[.,4];
pret_real=dataf[.,2];
graphset;
title("Volatilitatea implicita Heston");
_pltype=6;_pdate = ""; _pnum = 1; xlabel("Timp");ylabel("Volatilitate implicita");
xy(seqa(1,1,rows(volatilitate_implicita)),volatilitate_implicita);
graphset;
title("Eroare de evaluare Heston");
_pltype=6;_plctrl={-1,0};_pstype=14;_pdate = ""; _pcolor={10,12};_pnum = 1; xlabel("Moneyness (S/K)");ylabel("Eroarea");
xy(s./k,err~zeros(rows(err),1));
You have the advanced thank you from a masters student who is on the breach of going crazy with his dissertation paper 🙂 !