mvtcdfqmc error when mvncdf used within fsolve

3 次查看（过去 30 天）
Mark Whirdy
Mark Whirdy 2013-2-13
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;

请先登录，再回答此问题。

采纳的回答

Mark Whirdy
Mark Whirdy 2013-2-13

更多回答（1 个）

Tom Lane
Tom Lane 2013-2-13
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 个评论
Mark Whirdy
Mark Whirdy 2013-2-13
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 2013-2-15
Sometimes people use abs(H) in place of H in order to force H to be positive. This may be worth a try.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Biotech and Pharmaceutical 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by