Financial optimization of heston

Hi Matlab,

At the moment I'm pricing options, and I have written the code below to calibrate a pricing-model (the model is shown below here)

function [Call_SV] = Price_SV(Sigma, Kappa, Theta, eta, VIX, vT, vt, K, r, y)

Alpha= (1-exp(-Kappa.*((30./365))))./(Kappa.*(30./365)); 
Beta = (Theta).*(1-Alpha);
%c_2 = (2.*Kappa)./(eta.^2.*(1-exp(-Kappa.*(vT-vt))));
%w = c_2.*(V.*exp(-Kappa.*(vT-vt)));
%q = ((2.*Kappa.*Theta)./(eta.^2))-1; 
%I_q = besseli(q,2.*sqrt(w.*c_2.*y));
%gV = c_2.*exp(-w-c_2.*y).*((c_2./y)./w).^(q/2).*I_q
c_2 = (2.*Kappa)./((eta.^2).*(1-exp(-Kappa.*(vT-vt))));
w = c_2.*(((((VIX.^2)-Beta)./Alpha)).*exp(-Kappa.*(vT-vt)));
q = ((2.*Kappa.*Theta)./(eta.^2))-1; 
X = 2.*sqrt(w.*c_2.*((y.^2 -Beta)./Alpha));
I_q = besseli(q, X);
gV_inv = ((2.*y)./Alpha).*(c_2.*exp(-w-c_2.*((y.^2 -Beta)./Alpha)).*((c_2./((y.^2 -Beta)./Alpha))./w).^(q/2).*I_q);
INT =  integral(@(y) max(y-K,0).*gV_inv, 0, inf, 'RelTol',0,'AbsTol',1e-8);
Call_SV = exp(-r.*(vT-vt)).*INT;
%Feller 2*Kappa*Theta > eta
%con1 = 2.*Kappa.*Theta > eta;
%other condition y > sqrt(Beta) (otherwise zero)
%con2 = y > sqrt(Beta);
end 
%Test the model: Price_SV(0.61, 3.21, 0.19, 0.7, 1, 0.8, 0.7, 2, 0.1, 5)

Now my question is, how do I write code to calibrate my model?

I have to minimize sum of squarred residuals (using lqsnonlin), but how do I write the code?

Attempt:

%Price_SV(Sigma, Kappa, Theta, eta, VIX, vT, vt, K, r, y)
% Estimate model parameters
options = optimoptions('lsqnonlin','Display', 'iter-detailed', 'PlotFcn', 'optimplotresnorm', 'MaxIterations', 10000, 'TolFun', 10^-12 );
lsqnonlin(@Error, [0.5 0.5 0.5 0.7 0 0], [0 0 0 0 0], [10 10 10 10 10], options)

where Error is sum(data-Price_SV).^2

Any help appreciatet, really struling with this one, thanks alot!

Best regards Karnow

[Torsten]

Is y an input to your function

function [Call_SV] = Price_SV(Sigma, Kappa, Theta, eta, VIX, vT, vt, K, r, y)

or an integration variable

INT =  integral(@(y) max(y-K,0).*gV_inv, 0, inf, 'RelTol',0,'AbsTol',1e-8);

?

And a Relative Tolerance of 0 (thus an error of 0 in the integration) will usually not be possible to achieve.

[Simon Christensen]

Hi Torsten,

In this model, y is both: It is the spot price, price of the underlying asset. Once the integration has taken place, y has a value (and also before)

[Simon Christensen]

Perhaps it is not possible to have y placed in both places?

[Simon Christensen]

I just uploaded a image of the model:

[Torsten]

Perhaps it is not possible to have y placed in both places?

Yes, the constant value and the integration variable should be named differently. Otherwise your y in gV_inf is treated as a constant, and I think you want it as a function of y, don't you ? Same for X and Iq, I guess.

gV_inv = @(y) ((2.*y)./Alpha).*(c_2.*exp(-w-c_2.*((y.^2 -Beta)./Alpha)).*((c_2./((y.^2 -Beta)./Alpha))./w).^(q/2).*I_q);

[Simon Christensen]

Yes, that right :-), I would like it as a function of y :-)

[Simon Christensen]

But i think the right way is to renname y to spot and keep y in the integral :-)

[Torsten]

You integrate a function with respect to a variable y. All expressions that depend on the integration variable y must be defined as functions of y in your code. So each time you use the integration variable y in an expression, this expression must be defined as a function of y. Rename the variable y in the list of inputs to the function "Price_SV" into yconst and replace y by yconst in your code where you want to refer to this constant. Then change the other expressions still containing the integration variable y to functions.

[Simon Christensen]

Done, thanks :-). Do you know how I can calibrate the model to market data (with lsqnonlin)? I have lots of attempts, but nothing seems to work

[Torsten]

What are the parameters you want to fit ? What are the simulated data that you compute with your model and that you want to compare with the market data ?

lsqnonlin doesn't want you to return the sum of squared differences, but the differences themselves as a vector.

And if it were the sum of squared differences, your definition is wrong. Instead of

sum(data-Price_SV).^2

you had to define it as

sum((data-Price_SV).^2)

[Simon Christensen]

The parameters are Sigma, Kappa, Theta, eta, VIX, (all 1 x 1 elements) and I have data on vT (1 x1 matrix) , vt (1x1 matrix) K (10 x1 matrix, strike prices), r (1 x 1 matrix), y (1 x1 matrix) and then I have data on prices for options (10 x1 ), and I need to calibrate the prices i.e., minimize the squared difference between P_data and Price_SV, where P_data at data on market prices and Price_SV are prices from the model :-)

[Simon Christensen]

I have this five times, so 5 different sets of containing values on vT, vt, K, S, P_data, which need to be used to calibration (not sure if I shold calibrate to each set separately or all sets at once)

[Simon Christensen]

Really hope you can give a hint or help, worked with this for a very, very long time.

[Rena Berman]

(Answers Dev) Restored edit

[Dyuman Joshi]

@Rena Berman, can anything be done about the comments that were edited away as well?

[Rena Berman]

(Answers Dev) @Dyuman Joshi, I restored the comments.

[Dyuman Joshi]

Thank you, Rena!

Answer 1

[Torsten]

Assuming that Price_SV returns a vector of the same size as P_data, your call would look somehow like

fun =@(p)Price_SV(p(1), p(2), p(3), p(4), p(5), vT, vt, K, r, yconst) - P_data;
p0 = [0.61, 3.21, 0.19, 0.7, 1];
p = lsqnonlin(fun,p0)

or maybe

fun = @(p)arrayfun(@(K,P_data)Price_SV(p(1), p(2), p(3), p(4), p(5), vT, vt, K, r, yconst) - P_data,K,P_data)

if each element of K(i) gives the value of Price_SV(i) that is to be compared with P_data(i).

[Simon Christensen]

Thanks :-) assuming I use the first one, how do I get Price_SV to be a 10 x 1 matrix since my P-data are that size?

And also, every time I optimize I get 'Local minimum found. Optimization completed because the size of the gradient is less than1e-4 times the value of the function tolerance', the found values almost perfectly matches initial values, which is kinda annoying (maybe that because I only used one data point)

[Simon Christensen]

And also, do you know how I choose suitable values for initial values? Seems to have a rather large impact here :-)

[Torsten]

Thanks :-) assuming I use the first one, how do I get Price_SV to be a 10 x 1 matrix since my P-data are that size?

Within the function, make a loop over the number of P_data. Or use the second definition of fun.

And also, every time I optimize I get 'Local minimum found. Optimization completed because the size of the gradient is less than1e-4 times the value of the function tolerance', the found values almost perfectly matches initial values, which is kinda annoying (maybe that because I only used one data point)

Check in advance what your function returns by calling it with different parameters . I always got 0 as output for Price_SV which doesn't seem reasonable.

do you know how I choose suitable values for initial values?

I don't know about the background of your problem. Thus if somebody can find good initial guesses, it's you. Or try "MultiStart" which chooses a range of random values within specified ranges.

But first you should take care to make your function work properly.

[Simon Christensen]

Actually, I also observed that. Almost whenever ycon > K, the Price_SV gives zero, which doesn't make any sense. Hmmm. At the moment I'm not sure why that the case, but sure it has something do with Bessel function. Also, the returned values from estimation are way too high compared with other empirical findings, unfortunately.

[Torsten]

Almost whenever ycon > K, the Price_SV gives zero, which doesn't make any sense.

Seems you work with the ycon instead of y in the integration. max(ycon-K,0) gives 0 if ycon < K, thus you integrate the 0-function. This changes if y ranges from 0 to Inf.

[Simon Christensen]

Hmm, I will have to think about for a while, I'm not sure at the moment. I have written this code:

%Price_SV(Sigma, Kappa, Theta, eta, VIX, vT, vt, K, r, y)

%https://se.mathworks.com/help/fininst/calibrate-option-pricing-model-using-heston-model.html

K = [m2(:,3,1)] ;

Prices = [m2(:,2,1)]

%Construct objective function

for i = 1:length(K)

fun = @(Param) Prices(i) - Price_SV(Param(1), Param(2), Param(3), Param(4), Param(5), vT(2), vt(2), K(i), r(2), Spot(2));

end

p0 = [0.5 0.5 0.5 0.5 0.5];

lb = [0 0 0 0 0]

p = lsqnonlin(fun,p0, lb)

One last question please: If I want to add all sets of data (the 5 sets mentioned above), instead of only one set as shown in the code above, do you know have I could add them to the fun function?

So the p_data (which is Prices = [m2(:,2,1)] above), is called Prices1 = [m1(:,2,1)], Prices2 = [m2(:,2,1) and so on.

[Simon Christensen]

Above all, I must say really appreciate yor help! Writing a project in financial engineering, and almost no previous exposure to Matlab!

[Torsten]

for i = 1:length(K)
    fun = @(Param) Prices(i) - Price_SV(Param(1), Param(2), Param(3), Param(4), Param(5), vT(2), vt(2), K(i), r(2), Spot(2));
end

This will not work - if you want to use a loop, it has to be inside the function "Price_SV".

Concatenate all the data you have to two long vectors K and Prices and use

fun = @(Param) arrayfun(@(K,Prices) Prices - Price_SV(Param(1), Param(2), Param(3), Param(4), Param(5), vT(2), vt(2), K, r(2), Spot(2)),K,Prices);

This is the easiest solution.

See here for an example:

K = [1 4 7];
Prices = [1 4 9];
fun = @(Param)arrayfun(@(K,Prices) Prices - Price_SV(Param(1), K),K,Prices);
fun(2)

ans = 1×3

0 6 14

function Call_SV = Price_SV(p,K)
  Call_SV = p - K;
end

[Simon Christensen]

Hmmm, i tried this one

%Price_SV(Sigma, Kappa, Theta, eta, VIX, vT, vt, K, r)

%https://se.mathworks.com/help/fininst/calibrate-option-pricing-model-using-heston-model.html

rho_J =0

Mu_V =0

Mu_S =0

Lambda =0

S = [Spot(2), Spot(3) Spot(4), Spot(5), Spot(6)]

%concatenate prices, real and theo:

Prices_true= [m2(:,2,1); m3(:,2,1); m4(:,2,1); m5(:,2,1); m6(:,3,1)]

K_true= [m2(:,3,1); m3(:,3,1); m4(:,3,1); m5(:,3,1); m6(:,3,1)]

M2 = @(Param) SV_Heston(Param(1), Param(2), Param(3), Param(4), Param(5), rho_J, Mu_V, Mu_S, Lambda, S(1), K_true(1), r(2), vT(2), vt(2), 0.25);

M3 = @(Param) SV_Heston(Param(1), Param(2), Param(3), Param(4), Param(5), rho_J, Mu_V, Mu_S, Lambda, S(2), K_true(3), r(3), vT(3), vt(3), 0.25);

M4 = @(Param) SV_Heston(Param(1), Param(2), Param(3), Param(4), Param(5), rho_J, Mu_V, Mu_S, Lambda, S(3), K_true(4), r(4), vT(4), vt(4), 0.25);

M5 = @(Param) SV_Heston(Param(1), Param(2), Param(3), Param(4), Param(5), rho_J, Mu_V, Mu_S, Lambda, S(4), K_true(5), r(5), vT(5), vt(5), 0.25);

M6 = @(Param) [SV_Heston(Param(1), Param(2), Param(3), Param(4), Param(5), rho_J, Mu_V, Mu_S, Lambda, S(5), K_true(6), r(6), vT(6), vt(6), 0.25);

M = {M2, M3, M4, M5, M6}

%Construct objective function

fun = @(Param) arrayfun(@(K_true,Prices_true) Prices_true - M K_true,Prices_true);

p0 = [3 0.25 0.25 0.25 0.095];

lb = [0 0 0 0 0]

sol = lsqnonlin(fun, p0, lb)

But unfortunately, once it is a function handle datatype, it is not possible to subtract. Do you know how I can proceed?

[Torsten]

%Price_SV(Sigma, Kappa, Theta, eta, VIX, vT, vt, K, r)
%https://se.mathworks.com/help/fininst/calibrate-option-pricing-model-using-heston-model.html
rho_J  =0
Mu_V   =0
Mu_S   =0
Lambda =0
S      = [Spot(2), Spot(3) Spot(4), Spot(5), Spot(6)]
%concatenate prices, real and theo:
Prices_true= [m2(:,2,1); m3(:,2,1); m4(:,2,1); m5(:,2,1); m6(:,2,1)]
K_true= [m2(:,3,1); m3(:,3,1); m4(:,3,1); m5(:,3,1); m6(:,3,1)]
S_array = [S(1)*ones(size(m2(:,2,1)));...
           S(2)*ones(size(m3(:,2,1)));...
           S(3)*ones(size(m4(:,2,1)));...
           S(4)*ones(size(m5(:,2,1)));...
           S(5)*ones(size(m6(:,2,1)))];
r_array = [r(2)*ones(size(m2(:,2,1)));...
           r(3)*ones(size(m3(:,2,1)));...
           r(4)*ones(size(m4(:,2,1)));...
           r(5)*ones(size(m5(:,2,1)));...
           r(6)*ones(size(m6(:,2,1)))];
vT_array = [vT(2)*ones(size(m2(:,2,1)));...
           vT(3)*ones(size(m3(:,2,1)));...
           vT(4)*ones(size(m4(:,2,1)));...
           vT(5)*ones(size(m5(:,2,1)));...
           vT(6)*ones(size(m6(:,2,1)))];
vt_array = [vt(2)*ones(size(m2(:,2,1)));...
           vt(3)*ones(size(m3(:,2,1)));...
           vt(4)*ones(size(m4(:,2,1)));...
           vt(5)*ones(size(m5(:,2,1)));...
           vt(6)*ones(size(m6(:,2,1)))];
fun = @(Param)arrayfun(@(K_true,Prices_true,S_array,r_array,vT_array,vt_array)Prices_true-SV_Heston(Param(1), Param(2), Param(3), Param(4), Param(5), rho_J, Mu_V, Mu_S, Lambda, S_array, K_true, r_array, vT_array, vt_array, 0.25),K_true,Prices_true,S_array,r_array,vT_array,vt_array);
p0 = [3 0.25 0.25 0.25 0.095];
lb = [0 0 0 0 0]
sol = lsqnonlin(fun, p0, lb)