matlab doesn't calibrate via lsqnonlin when initial values depart from the actual parameters

Hi guys,

I have written a program that calibrates two parameters from the option prices. I have simulated 10 option prices for known parameters and now try to check my calibration function.

If I set initial values of the paramters close to the actual one (the difference is not more that 0.4) then it calibrates the right parameters in 4 seconds. When I set the initial values slightly different from the initial (e.g. I set it equal to 2 when actual is 1), or when I change the upper and lower bounds and run the calibration, matlab calibration doesn't converge - it is always busy and I have to kill the process. ( I even left it for 1 night, the result was still "busy")

Please help me, changing the accuracy and number of function evaluation doesn't help!

[Matt J]

Let's have a look at the code.

[Daniel]

%%%%%%%%%%Calibration part  %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%initial parameters and lower\upper bounds %%

x0=[0.2,-1.5];
lb = [0.001,-2];
ub = [1.2,0];

%%%%%%%%%%reading market data from from xls file %%%%%%%%%%%%

mktdata=xlsread('marketdata');
K = mktdata(2:11,1);
T = mktdata(1,2);
mktprice = mktdata(2:11,2);
S=50;

%%%%%%%%%%%%%Calibration - Determining the parameters %%%%%%%%%

options=optimset('TolFun',1e-3);
x = lsqnonlin(@LSQDiff,x0,lb,ub,options);
toc;

Then follows the code for each function:

    function [ Diff ] = LSQDiff( x )

%%%%%%%%%%%%%%%%%%%Allocate memory %%%%%%%%%%%%%%%%%
Diff=zeros(length(K),length(T));

%%%%%%%%%%Calculate the matrix of differences %%%%%%%%%%%%%%%
    for i=1:length(T)
     Diff(:,i)= mktprice(:,i)-Putprice(S,K,T(i),x(1),x(2));
    end
reshape(Diff,length(T)*length(K),1);

next one

function [Put]=Putprice(S,K,T,c,beta)
%%%%%%%declare the parameters for previously calibrated a(t) and b(t)
a1=0;  a2=10;  b1=0;  b2=0.02;  r1=0;  r2=0.05;  q1=0;  q2=0;  syms u;

%%%%%%%%%%%%%%%%%specify r(t) and q(t)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
r=r1*u+r2;
q=q1*u+q2;
rfunc=@(u)r1*u+r2;
qfunc=@(u)q1*u+q2;

%%%%%%%%%%%%%%%%%specify a(t) and b(t) functions %%%%%%%%%%%%%%%%%%%%%
a=a1*u+a2;
b=b1*u+b2;
alpha=r-q+b;

bfunc=@(u) b1*u+b2;
taufunc=@(u)a^2*exp(-2*abs(beta)*int(alpha,u,0,u));
rbfunc=@(u) r1*u+r2+b1*u+b2;
alphafunc=@(u) r1*u+r2+b1*u+b2-q1*u-q2;

res=taufunc(u);
resfunc=matlabFunction(res);

%%%%%%%%%%%%%%%%%compute integral for tau!%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tau=integral(@(u)resfunc(u),0,T);

%%%%%%%%%%%%%%%%%%continue with the pricing %%%%%%%%%%%%%%%%%%%%%%%%%
nuplus=(c+0.5)/abs(beta);
delta=2*(nuplus+1);
x=1/abs(beta)*S^abs(beta);
k=1/abs(beta).*(K.^abs(beta)).*exp(-abs(beta).*integral(alphafunc,0,T));

%%%%%%%%%%%%%%Clean Price of the option %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Put_clean=exp(-integral(rbfunc,0,T)).*K.*(x^2/tau)^(1/(2*abs(beta))).*Phiminus(-0.5/abs(beta),k.^2/tau,delta,x^2/tau)...
    -exp(-integral(qfunc,0,T))*S.*Phiminus(0,k.^2/tau,delta,x^2/tau);

%%%%%%%%%%%%%Default part of the option %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  Def=1-exp(-integral(bfunc,0,T))*(x^2/tau)^(1/(2*abs(beta)))*momentum(-1/(2*abs(beta)),delta,x^2/tau);
  Put_def=K.*exp(-integral(rfunc,0,T)).*Def;
  Put=Put_clean+Put_def;
  end

Last one

function y=Phiminus(p,k,delta,alpha)
temp1=2^p*exp(-alpha/2);
tol=10^(-60);
%if I set here smaller tolarenace value, my calibrated parameters significantly differ from the real ones (0.6 insteadof 1). But matlab does it in 4 seconds, so it should be ok.
n=0;
x=1;
%f=zeros(10,length(k));
    while(x>tol)
    f(n+1,:)=temp1*(alpha/2)^n*gamma(p+delta/2+n).*...
gammainc(k/2,p+delta/2+n)./(gamma(n+1)*gamma(delta/2+n));
    x=f(n+1,1);
    n=n+1;
    end
z=sum(f);
y=z';

end

Hope it's readable and not too long. I would really appreciate any help.

[Matt J]

OK. Well, given it's complexity, it might be easier to help if you write the function you're trying to minimize in mathematical form, rather than code form.

For now, I'll just remark that non-differentiable operations like abs(beta) are probably "illegal", since lsqnonlin uses smooth algorithms.

[Daniel]

http://postimg.org/image/8bhnhe5h9/

this is the function that has two parameters beta and с (с is inside the delta)

r(u) and b(u) are constants

Answer 1

[Matt J]

It sounds like you have a highly ill-conditioned problem which could be slowing convergence. The non-differentiability menationed in my Comment above might also be a factor.

Because it's only a 2-variable problem, you can probably investigate what's going on by making a surf() plot of the objective function (the squared residual) and see if it's shaped funny.

It might also help to terminate the algorithm early, using MaxIterations or similar, and seeing where it's getting stuck. Then, look for this location on the surf plot and see if the local shape of the surface is behaving oddly there.

[Matt J]

One other comment. You are constraining the unknown parameter x(2), or "beta" to be negative, lying between lb(2)=-2 and ub(2)=0.

This is strange, because the objective function only ever uses abs(beta). Thus, it makes me wonder why you don't just get rid of the abs() and set lb(2)=0 and ub(2)=2. That would take care of the differentiability issue.

[Matt J]

One other comment (a question really). How did you arrive at TolFun=1e-3? It is not at all obvious to me from looking at your code why this would be a better order of magnitude for TolFun than the default 1e-6.

[Daniel]

The problem is that when I run the code it simply doesn't stop. The status is always 'busy'. Even when I change the number of iterations to five, it still doesn't give any result. I played with the tolerance level just to obtain any result, but it doesn't work neither for high nor for low tolerance level.

I also tried to run the optimization with fmincon - situation is the same, no result. I have to kill the process manually.

[Matt J]

Have you tried executing

 x0=[0.2,-1.5];
 LSQDiff(x0)

separately to see if it can even evaluate the first point bug free? If it hangs, you should set a breakpoint in the first line of LSQDiff and then use DBSTEP to step line by line until you find the precise line within the function where it starts to hang.

[Daniel]

You are right. For some values that are far from the initial ones LSQDiff hangs.

[Matt J]

Did you locate which line?

[Daniel]

it stops here

%%%%%%%%%%%%%Default part of the option %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Def=1-exp(-integral(bfunc,0,T))*...(x^2/tau)^(1/(2*abs(beta)))*momentum(-1/(2*abs(beta)),delta,x^2/tau);

The problem is in hypergeom function, for some values of the parameters it gives NaN

That's actually very bad for the calibration, because I don't know which upper and lower bounds and initial values I should set...