MATLAB Answers

## mvtcdfqmc error when mvncdf used within fsolve

Asked by Mark Whirdy

### Mark Whirdy (view profile)

on 13 Feb 2013
Accepted Answer by Mark Whirdy

### Mark Whirdy (view profile)

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;```

## 2 Answers

### Mark Whirdy (view profile)

Answer by Mark Whirdy

on 13 Feb 2013
Accepted answer

### Tom Lane (view profile)

Answer by Tom Lane

### Tom Lane (view profile)

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.

Mark Whirdy

### Mark Whirdy (view profile)

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

### Tom Lane (view profile)

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.

#### Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply today