Solving system on non-linear equations (containing bivariate cumulative normal distribution function)

Question

I am experiencing some struggles in solving a system of 4 equations with 4 unknows. The equations are non-linear and contain a bivariate normal cumulative distribution function. I tried multiple numerical optimization packages (fmincon, fminsearch, fsolve,..) but they all return the 'simple solution' of L = 1, sigma_A = 0, alpha1 = 1 and alpha2=1. Which is not a realistic solution. I would expect L to be in the range of 0.1 to 0.95, and siga_A to be nonzero.

Many thanks in advance

function [L_t,sig_A,s,p,DTD] = Hull_Clean(Impl_Vol1, Impl_Vol2, tau1, tau2, Px_Last1, Px_Last2, STRIKE1, STRIKE2,t,T,r) 



%Impl_Vol1=0.21; Impl_Vol2=0.31; tau1=30/252; tau2=30/252; Px_Last1=76.16; Px_Last2=76.16; STRIKE1=70; STRIKE2=65;t=0;T=5;r=0.03;

%% RESHAPE VARIABLES TO 3D TO FACILITATE MATRIXWISE BIVAR NEWTON-REPHSON 

[n,m] = size(Impl_Vol1);
nm = n*m;

Impl_Vol1 = reshape(Impl_Vol1,1,1,nm); Impl_Vol2 = reshape(Impl_Vol2,1,1,nm); STRIKE1 = reshape(STRIKE1,1,1,nm); STRIKE2 = reshape(STRIKE2,1,1,nm);
tau1 = reshape(tau1,1,1,nm); tau2 = reshape(tau2,1,1,nm); Px_Last1 = reshape(Px_Last1,1,1,nm); Px_Last2 = reshape(Px_Last2,1,1,nm);



%% Calculate known variables CONTROLEER ALLES NOG, OOK MET GESKE 1979 ARTIKEL OP CORRECTHEID FORMULES!!!!

Kappa1 = STRIKE1 .* exp(-r .* tau1) ./ Px_Last1; % GEBRUIK ALTIJD STRIKE EN Px_Last per stock of hele bedrijf  
Kappa2 = STRIKE2 .* exp(-r .* tau2) ./ Px_Last2;

d1_star1 = (-log(Kappa1) ./ (Impl_Vol1 .* sqrt(tau1))) + (0.5 .* Impl_Vol1 .* sqrt(tau1)); d2_star1 = d1_star1 - (Impl_Vol1 .* sqrt(tau1));
d1_star2 = (-log(Kappa2) ./ (Impl_Vol2 .* sqrt(tau2))) + (0.5 .* Impl_Vol2 .* sqrt(tau2)); d2_star2 = d1_star2 - (Impl_Vol2 .* sqrt(tau2));

FNN1 = (Kappa1 .* fcnN(-d2_star1)) - fcnN(-d1_star1);
FNN2 = (Kappa2 .* fcnN(-d2_star2)) - fcnN(-d1_star2);

%% Define functions of unknown variables

d1 = @(L_t,sig_A,C)((1./(sig_A(:,:,C).*sqrt(T(:,:,C)-t))).*(-log(L_t(:,:,C)) + ((0.5.*sig_A(:,:,C).^2).*(T(:,:,C)-t))));
d2 = @(L_t,sig_A,C)((1./(sig_A(:,:,C).*sqrt(T(:,:,C)-t))).*(-log(L_t(:,:,C)) - ((0.5.*sig_A(:,:,C).^2).*(T(:,:,C)-t))));

a1 = @(alpha, sig_A, tau, C)((1./(sig_A(:,:,C).*sqrt(tau(:,:,C)-t))).*(-log(alpha(:,:,C)) + ((0.5.*sig_A(:,:,C).^2).*(tau(:,:,C)-t))));  % MOETEN HIER ZEKER - en + zo? ?CHECK
a2 = @(alpha, sig_A, tau, C)((1./(sig_A(:,:,C).*sqrt(tau(:,:,C)-t))).*(-log(alpha(:,:,C)) - ((0.5.*sig_A(:,:,C).^2).*(tau(:,:,C)-t))));

d1_tau = @(L_t, alpha, sig_A, tau, C)((1./(sig_A(:,:,C).*sqrt(T(:,:,C)-tau(:,:,C)))).*(-log(L_t(:,:,C) ./ alpha(:,:,C)) + ((0.5.*sig_A(:,:,C).^2).*(T(:,:,C)-tau(:,:,C)))));
d2_tau = @(L_t, alpha, sig_A, tau, C)((1./(sig_A(:,:,C).*sqrt(T(:,:,C)-tau(:,:,C)))).*(-log(L_t(:,:,C) ./ alpha(:,:,C)) - ((0.5.*sig_A(:,:,C).^2).*(T(:,:,C)-tau(:,:,C)))));


%% System of nonlinear equations 

Eq1 = @(L_t, sig_A, alpha, tau, C)(((alpha.*fcnN(d1_tau(L_t, alpha, sig_A, tau, C))) - (L_t(:,:,C).*fcnN(d2_tau(L_t, alpha, sig_A, tau, C)) )...
    )./(fcnN(d1(L_t,sig_A,C))-(L_t(:,:,C).*fcnN(d2(L_t,sig_A,C)) ))); % Of -1 kappa hier zodat hij == 0 kan solven? 

Eq1_1 = @(L_t, sig_A, alpha1, C)((((alpha1(:,:,C).*fcnN(d1_tau(L_t, alpha1, sig_A, tau1, C))) - (L_t(:,:,C).*fcnN(d2_tau(L_t, alpha1, sig_A, tau1, C)) )...
    )./(fcnN(d1(L_t,sig_A,C))-(L_t(:,:,C).*fcnN(d2(L_t,sig_A,C)) ))) - Kappa1(:,:,C)) ; 

Eq1_2 = @(L_t, sig_A, alpha2, C)((((alpha2(:,:,C).*fcnN(d1_tau(L_t, alpha2, sig_A, tau2, C))) - (L_t(:,:,C).*fcnN(d2_tau(L_t, alpha2, sig_A, tau2, C)) )...
    )./(fcnN(d1(L_t,sig_A,C))-(L_t(:,:,C).*fcnN(d2(L_t,sig_A,C)) ))) - Kappa2(:,:,C)) ; 

Eq2_1 = @(L_t, sig_A, alpha1, C)(((L_t(:,:,C)*fcnM(-a2(alpha1, sig_A, tau1, C),d2(L_t,sig_A,C), - sqrt(tau1(:,:,C) ./ T(:,:,C)) ))  - ...
    fcnM(-a1(alpha1, sig_A, tau1, C), d1(L_t,sig_A,C), - sqrt(tau1(:,:,C) ./T(:,:,C))) + ...
    (Kappa1(:,:,C) .* fcnN(-a2(alpha1, sig_A, tau1, C)) .* (fcnN(d1(L_t, sig_A, C)) - (L_t(:,:,C).*fcnN(d2(L_t, sig_A,C))) ))) - ...
    (FNN1(:,:,C) .* (fcnN(d1(L_t, sig_A, C)) - (L_t(:,:,C).*fcnN(d2(L_t, sig_A,C))) ) )) ; 

Eq2_2 = @(L_t, sig_A, alpha2, C)(((L_t(:,:,C)*fcnM(-a2(alpha2, sig_A, tau2, C),d2(L_t,sig_A,C), - sqrt(tau2(:,:,C) ./ T(:,:,C)) ))  - ...
    fcnM(-a1(alpha2, sig_A, tau2, C), d1(L_t,sig_A,C), - sqrt(tau2(:,:,C) ./T(:,:,C))) + ...
    (Kappa1(:,:,C) .* fcnN(-a2(alpha2, sig_A, tau2, C)) .* (fcnN(d1(L_t, sig_A, C)) - (L_t(:,:,C).*fcnN(d2(L_t, sig_A,C))) ))) - ...
    (FNN2(:,:,C) .* (fcnN(d1(L_t, sig_A, C)) - (L_t(:,:,C).*fcnN(d2(L_t, sig_A,C))) ) )) ; 



%% Solve system on non linearr equiations 
opts = optimset('tolfun',0,'tolx',0,'maxfun',Inf);
opts=optimset('Algorithm','Levenberg-Marquardt');
x0 = [0.1, 0.5 ,0.8,0.9];

fun = @(x)[Eq1_1(x(1), x(2), x(3), true); Eq1_2(x(1), x(2), x(4), true); Eq2_1(x(1), x(2), x(3), true); Eq2_2(x(1), x(2), x(4), true)];

[VALUES,fval] = fsolve(fun, x0, opts)
L_t = VALUES(1); sig_A = VALUES(2); alpha1 = VALUES(3); alpha2 = VALUES(4); 




%% SOLVE FOR CREDIT SPREAD [s]
d1 = @(L_t,sig_A,C)((1/(sig_A *sqrt(T -t)))*(-log(L_t ) + ((0.5*sig_A^2)*(T -t))));
d2 = @(L_t,sig_A,C)((1/(sig_A *sqrt(T -t)))*(-log(L_t ) - ((0.5*sig_A^2)*(T -t))));

s = - log(fcnN(d2(L_t, sig_A, true)) + (fcnN(-d1(L_t, sig_A, true)) / L_t)) / (T - t);

p = fcnN(-d2(L_t,sig_A,true)); % P(A_t < K) = N(-d_m)

DTD = d2(L_t,sig_A,true);

end


%
%% SUBFUNCTIONS

function p=fcnN(x)
p=0.5*(1.+erf(x/sqrt(2)));
end
%
function p=fcnn(x)
p=exp(-0.5*x^2)/sqrt(2*pi);
end

function Y = inv3d(X)
    Y = -X;
    Y(2,2,:) = X(1,1,:);
    Y(1,1,:) = X(2,2,:);
    detMat = 1/(X(1,1,:)*X(2,2,:) - X(1,2,:)*X(2,1,:));
    detMat = detMat(ones(1,2),ones(2,1),:);
    Y = detMat*Y;
end

function p=fcnM(a,b,rho)
X = [a;b];
mu = [0;0];
sigma = [1, rho; rho, 1];

p = mvncdf(X,mu,sigma);
end

function p=fcnM1(a,b,rho)
    if(a <= 0 && b <= 0 && rho <= 0)

    aprime = a/(sqrt(2*(1-(rho^2))));
    bprime = b/(sqrt(2*(1-(rho^2))));
    A = [0.3253030 0.4211071 0.1334425 0.006374323];
    B = [0.1337764 0.6243247 1.3425378 2.2626645];

    F = 'exp(aprime*(2*x - aprime)+ (bprime*(2*y - bprime)) + (2*rho *(x - aprime)*(y-bprime)))'; 
    t = 0;

    for i=1:4
        for j=1:4
            x = B(i);
            y = B(j);
            t = t + A(i)*A(j)*eval(F);

        end
    end

    p = (sqrt(1-rho^2)/pi) * t;

elseif (a * b * rho <= 0)

        if (a <=0 && b >=0 && rho >=0)
            p = normcdf(a) - fcnM1(a,-b,-rho);
        elseif (a >=0 && b <=0 && rho >=0)
                p = normcdf(b) - fcnM1(-a,b,-rho);
        elseif (a >=0 && b >=0 && rho <=0) %modified here at 1:45 AM
                p = normcdf(a) + normcdf(b) - 1 + fcnM1(-a,-b,rho);
        end
elseif  a*b*rho > 0;
    %Could not use the In-Built function sign(x) because it is +1 if x>=0
    %not just x>0 as in Matlab.

    if(a >= 0), 
        asign =1 ;
    else
        asign = -1;
    end

    if(b >= 0), 
        bsign =1 ;
    else
        bsign = -1;
    end

    rho1 = (rho*a - b)*asign/(sqrt(a^2 - (2*rho*a*b) + b^2));
    rho2 = (rho*b - a)*bsign/(sqrt(a^2 - (2*rho*a*b) + b^2));
    delta = (1-(asign*bsign))/4;

        p = fcnM1(a,0,rho1) + fcnM1(b,0,rho2) - delta ;
    end
end

Did you try to include bounds to your parameters (`lb` and `ub` in `fmincon`)? Or any other constraints on your parameters/equations? Or different initial guesses? — rinkert, Jun 06 '19 at 11:43
yes, I did. For example if is set lb = [0, 0.02, 0, 0] and ub = [ 0.98, inf, inf, inf], it will just return something like [0.98, 0.02, 1, 1] — Python_Donald, Jun 06 '19 at 12:19
And what if you use random initial conditions (that are kind of reasonable perhaps), or use a [genetic algorithm](https://nl.mathworks.com/help/gads/ga.html) to do so for you. — rinkert, Jun 06 '19 at 12:22
Unfortunately I do not have the Global Optimization Toolbox. Hence I cannot use the ga() function — Python_Donald, Jun 06 '19 at 12:45
Alternatively, you could do a grid search, find 'good' candidates, and use these as starting points for the optimizer. This also gives you some idea of the influence of parameters on the cost function, though this is really easy to visualize with 4 parameters. — rinkert, Jun 06 '19 at 12:55
I tried all sorts of starting values, this hardly changes the outcome of the optimization. — Python_Donald, Jun 06 '19 at 13:13

Adam · Answer 1 · 2019-06-06T18:54:38.843

Use lsqnonlin since your function output is a scalar
Let sigma_A be a free variable

sigma_A_lb = -inf; sigma_A_ub = inf

The code is as follow

% Given initials 
Impl_Vol1=0.21; Impl_Vol2=0.31; tau1=30/252; tau2=30/252; Px_Last1=76.16; Px_Last2=76.16; STRIKE1=70; STRIKE2=65;t=0;T=5;r=0.03;

% Change this section

x0 = [0.1, 0.02 ,0.8,0.9];
lb = [0.1, -inf, 0, 0]; 
ub = [0.98, inf, inf, inf];
fun = @(x)[Eq1_1(x(1), x(2), x(3), true); Eq1_2(x(1), x(2), x(4), true); Eq2_1(x(1), x(2), x(3), true); Eq2_2(x(1), x(2), x(4), true)];

options = optimoptions('lsqnonlin','FunctionTolerance',1e-10);
[VALUES,fval] = lsqnonlin(fun,x0,lb,ub, options) ;

Result

VALUES =
    |-------------------------------------|
    |   L    |  sigma_A| alpha1  | alpha2 | 
    |-------------------------------------|
    |0.9757  | 0.0052  | 0.9979  |  0.9963|
    |-------------------------------------|

fval =

   8.5170e-10

Checking with conditions given in comments

% Initial conditions given in comments

x0 = [0.1, 0.1 ,0.8,0.9];
lb = [0.1, -inf, 0, 0];
ub = [0.98, inf, inf, inf];

Result

VALUES =
    |-------------------------------------|
    |   L    |  sigma_A| alpha1  | alpha2 | 
    |-------------------------------------|
    |0.8770  | 0.0265  | 0.9895  |  0.9813|
    |-------------------------------------|

fval =

   2.2257e-08

Hi, thanks for your answer. However, I think that this solution just depends on your lower bound of sig_A in this case. [0.8241 0.0100 0.9619 0.9930] is the solution for x0 = [0.1, 0.1 ,0.8,0.9]; lb = [0.1, 0.01, 0, 0]; ub = [0.98, inf, inf, inf]. So sigma_A will just be equal to your lower bound. — Python_Donald, Jun 06 '19 at 14:49
Just let `sigma_A` be free, then the solution obtained will depend only on the initial guess. Check the new update answer — Adam, Jun 06 '19 at 17:05

Solving system on non-linear equations (containing bivariate cumulative normal distribution function)

1 Answers1