mvtcdfqmc error when mvncdf used within fsolve

4 views (last 30 days)
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;

Accepted Answer

Mark Whirdy
Mark Whirdy on 13 Feb 2013

More Answers (1)

Tom Lane
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
Tom Lane
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.

Sign in to comment.

Categories

Find more on Systems of Nonlinear Equations in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!