Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

mvtcdfqmc error when mvncdf used within fsolve

Asked by Mark Whirdy on 13 Feb 2013

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;

0 Comments

Mark Whirdy

2 Answers

Answer by Mark Whirdy on 13 Feb 2013
Accepted answer

0 Comments

Mark Whirdy
Answer by Tom Lane on 13 Feb 2013

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.

2 Comments

Mark Whirdy on 13 Feb 2013

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 on 15 Feb 2013

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

Tom Lane

Contact us