非线性方程组的求解系统(包含二元累积正态分布函数)

Question

我在解决 4 个未知的方程组时遇到了一些困难。这些方程是非线性的，包含二元正态累积分布函数。我尝试了多个数值优化包(fmincon，fminsearch，fsolve,..(，但它们都返回L = 1，sigma_A = 0，alpha1 = 1和alpha2 = 1的"简单解"。这不是一个现实的解决方案。我希望 L 在 0.1 到 0.95 的范围内，siga_A为非零。

提前非常感谢

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

Answer 1

• 使用lsqnonlin，因为函数输出是标量

• 设sigma_A为自由变量

sigma_A_lb = -inf; sigma_A_ub = inf

代码如下

% 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) ;

---

结果

VALUES =
    |-------------------------------------|
    |   L    |  sigma_A| alpha1  | alpha2 | 
    |-------------------------------------|
    |0.9757  | 0.0052  | 0.9979  |  0.9963|
    |-------------------------------------|
fval =
   8.5170e-10

---

---

检查注释中给出的条件

% 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];

结果

VALUES =
    |-------------------------------------|
    |   L    |  sigma_A| alpha1  | alpha2 | 
    |-------------------------------------|
    |0.8770  | 0.0265  | 0.9895  |  0.9813|
    |-------------------------------------|
fval =
   2.2257e-08