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

Hi all,

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 = 0, alpha1 = 1 and alpha2=1. Which is not a realistic solution.

function [L_t,sig_A] = Hull_Clean(Impl_Vol1, Impl_Vol2, tau1, tau2, Px_Last1, Px_Last2, STRIKE1, STRIKE2,t,T,r) %Impl_Vol1, Impl_Vol2, tau1, tau2, Px_Last1, Px_Last2, STRIKE1, STRIKE2,  
%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); 


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