mvtcdfqmc error when mvncdf used within fsolve

Save the function below & call it a couple of times by

E_t = 2000;
sig_E = 0.5;
K = [10*ones(5,1);1010];
T = [1:6]';
t = 0;
r = 0.1;
[A_t,sig_A,D_t,s,A_T] = ic_calcGeskeStructModel(E_t,sig_E,K,T,t,r);

Sometimes I get the error from the quasi-monte-carlo algo mvtcdfqmc, (but only sometimes). Problem is due to fsolve chosing negative values of x which causes a log(-x)=complex and erfc() to throw an error at the complex number. Seems I need a constrained non-linear equation-solver (to constrain all x>0), any suggestions?

" ??? Error using ==> erfc Input must be real and full.

Error in ==> mvtcdfqmc>normp at 185 p = 0.5 * erfc(-z ./ sqrt(2));

.... etc

"

********************************************************************* *********************************************************************

    function [A_t,sig_A,D_t,s,A_T] = ic_calcGeskeStructModel(E_t,sig_E,K,T,t,r)
    %
    %
    %   Input:
    %       E_t     Equity Value of Firm          [1x1]
    %       sig_E   Equity Implied Volatility     [1x1]
    %       X       Cashflows (Coupon&Principal)  [nx1]
    %       T       Time of Cashflows (years)     [nx1]
    %       t       Current time                  [1x1]  
    %       r       Spot Rates                    [nx1]
    %
    %   Output:
    %       A_t     Value of Firm's Assets        [1x1]
    %       sig_A   Volatility of Firm's Assets   [1x1]
    %       D_t     Value of Firm's Debt          [1x1]
    %       s       Credit Spread                 [1x1]  
    %       A_T     Debt Default Barrier Vector   [nx1]
    %
    %
    %
    %   Example:
    %       E_t = 2000;
    %       sig_E = 0.5;
    %       K = [10*ones(5,1);1010];
    %       T = [1:6]';
    %       t = 0;
    %       r = 0.1;
    %       [A_t,sig_A,D_t,s,A_T] = ic_calcGeskeStructModel(E_t,sig_E,K,T,t,r);
    %
    % References:
    %   [1] Geske R., 1979, "The Valuation of Compound Options", 
    %       Journal of Financial Economics 7, pp. 63-81.
    %       http://bbs.cenet.org.cn/uploadImages/20035291315398755.pdf
    %   [2] Geske R., 1977, "The Valuation of Corporate Liabilities as Compound Options", 
    %       Journal of Financial and Quantitative Analysis, Vol. 12, No. 4, November, pp. 541-552.
    %       http://www.defaultrisk.com/pa_price_25.htm
    %   [3] Thomassen L. and Van Wouwe M., 2001, "The n-fold Compound Option",
    %       Working Paper 2001-041, UA Faculty of Applied Economics, Antwerp.
    %       http://ideas.repec.org/p/ant/wpaper/2001041.html
    %

    %%CREATE CORRELATION MATRIX

    % rho = sqrt(T_i/T_j) i<j
    rho = T(:,ones(size(T,1),1)); % Repmat the Cashflow Time Vector
    rho = triu(rho./rho') + triu(rho./rho',+1)'; % Create Symmetrical Matrix 
    rho = rho.^0.5;

CREATE BLACK-SCHOLES-PARAMETER VECTORIZED ANONYMOUS FUNCTIONS

    % A_t   [1x1]
    % sig_A [1x1]
    % H     [mx1]; m counts along the cashflows [1 <= m <= n]
    % T     [mx1]
    % t     [mx1]
    % r     [mx1]

    d_p = @(A_t,sig_A,H,T,r,t)((1./(sig_A.*sqrt(T-t))).*(log(A_t./H) + (r + 0.5*sig_A^2).*(T-t)));
    d_m = @(A_t,sig_A,H,T,r,t)((1./(sig_A.*sqrt(T-t))).*(log(A_t./H) + (r - 0.5*sig_A^2).*(T-t)));

    %%SOLVE FOR ASSET VALUE, VOLATILITY & DEBT-BARRIERS [A_t, sig_A, H[n]]

    tol = 1e-3;
    maxfunevals = 1e5;
    o = optimset('MaxFunEvals',maxfunevals,'TolFun',tol); % Default values are two fine for tractable solution
    p = 2 + size(K,1);
    g = @(x)objfun(x,d_p,d_m,T,K,t,r,E_t,sig_E,rho,p,o);
    x0 = [E_t+sum(K./(1+r)); sig_E*E_t/(E_t+sum(K./(1+r))); K]; 

    options = optimset('Display','iter','Algorithm','levenberg-marquardt'); 
    [x,fval] = fsolve(g,x0,options); % Solve using 

    A_t = x(1);
    sig_A = x(2);
    A_T = x(3:p);

    %%SOLVE FOR FIRM DEBT VALUE [D_t]

    D_t = A_t - E_t; % Balance Sheet

    %%SOLVE FOR CREDIT SPREAD [s]

    % Define a Credit Spread [s] as
    % D_t = K1*exp(-(r1+s)*(T1-t)) + K2*exp(-(r2+s)*(T2-t))

    spread = @(s)(K'*exp(-(r+s)*(T-t)) - D_t);
    s = fzero(spread,0);

    %
    %
    %
    %%OBJECTIVE FUNCTION FOR SIMULTANEOUS SOLUTION OF ASSET VAL, VOL, & DEBT-BARRIER-VECTOR

    function F = objfun(x,d_p,d_m,T,K,t,r,E_t,sig_E,rho,p,o)

    % x(1)   = A_t
    % x(2)   = sig_A
    % x(3:p) = H

    F = [...
        K(1:end-1) + (K(2:end).*exp(-r.*diff(T)).*normcdf(d_m(x(1),x(2),x(4:p),T(2:end),r,T(1:end-1)),0,1)) + x(4:p).*normcdf(d_p(x(1),-x(2),x(4:p),T(2:end),r,T(1:end-1)),0,1) - x(4:p);
        x(1).*mvncdf(d_p(x(1),x(2),x(3:p),T,r,t)',[],rho,o) - sum(K.*mvncdf(trilinf(d_p(x(1),x(2),x(3:p,ones(6,1)),T(:,ones(6,1)),r,t)'),[],rho,o)) - K(end).*exp(-r*T(end))*normcdf(d_m(x(1),x(2),x(p),T(end),r,t),0,1) - E_t;
        x(1).*x(2).*mvncdf(d_p(x(1),x(2),x(3:p),T,r,t)',[],rho,o) - E_t.*sig_E;...
        ];

    function Y = trilinf(X)

    Y = tril(X);
    Y(Y==0) = Inf;

Answer 1

[Mark Whirdy]

Got it

< http://www.mathworks.com/matlabcentral/answers/7549>

Answer 2

[Tom Lane]

When I try this, it's because H becomes negative and so log(A_t./H) becomes complex.

You may find it helpful to set a breakpoint at the line that defines d_p, and specify that it is inside the anonymous function, and make it conditional on any(H(:)<0). That could help you figure out what is going on.

[Mark Whirdy]

Hi Tom, thanks.

Indeed all members of x should be positive (in context) and where fsolve pushes them negative, the log(-x[i]) produces a complex number which erfc can't handle.

Seems I need CONSTRAINED (all x>0) nonlinear equation solver, fmincon requires a scalar function though & F is vector.

Any other workarounds?

[Tom Lane]

Sometimes people use abs(H) in place of H in order to force H to be positive. This may be worth a try.