0

Stuck on numerically evaluating this integral and only have experience in doing simple quadrature schemes so bear with me if this looks amateurish. The nested integral is below using Gauss-Hermite(-inf,+inf) and a variant of Gauss Laguerre scheme (0,+inf) - over a gaussian and gamma(v/2,2) densities. I get to the point where the I m nearly finished the integration, intermediate steps look good, but I m stuck on how to combine the weights to evaluate the overall integral. I'd be really grateful for suggestions in modifing the code/other idea to code a better quadrature scheme that solves the problem. \begin{equation} \int^{\infty}{-\infty}\int^{\infty}{0} \prod_{i=1}^n\Phi \left(\frac{\sqrt{w/v}\,C{i}-a{i}Z}{\sqrt{1a {i^2}}\right)f{z}(Z)f{w}(W)dwdz \end{equation}

% script defines nodes, weights, parameters then calls one main and one subfunction

rho=0.3; nfirms=10; h=repmat(0.1,[1,nfirms]); T=1; R=0.4; v=8; alpha=v/2;GaussPts=15;

% Quadrature nodes - gaussian and gamma(v/2) from Miranda and Fackler CompEcon
    % toolbox
[x_norm,w_norm] = qnwnorm(GaussPts,0,1);
[x_gamma,w_gamma] = qnwgamma(GaussPts,alpha);
L_mat=zeros(nfirms+1,GaussPts);

for i=1:1:GaussPts;
     L_mat(:,i) = TC_gamma(x_norm(i,:),x_gamma(i,:),h,rho,T,v,nfirms);
end; 

w_norm_mat= repmat(w_norm',nfirms+1,1);
w_gamma_mat = repmat(w_gamma',nfirms+1,1);
% need to weight L_mat by the gaussian and chi-sq i.e, (gamma v/2,2)?
ucl = L_mat.*w_norm;%?? HERE
ucl2 = sum(ucl.*w_gamma2,2);% ?? HERE


function [out] = TC_gamma(x_norm,x_gamma,h,rho,T,v,nfirms)
% calls subfunction feeds into recursion

qki= Vec_CondPTC_gamma(x_norm,x_gamma,h,rho,T,v)' ;

fpdf=zeros(nfirms+1,nfirms+1);
% start at the first point on the tree
 fpdf(1,1)=1; 
    for i=2:nfirms+1 ;
     fpdf(1,i)=fpdf(1,i-1)*(1-qki(:,i-1));
       for j=2:nfirms+1;
       fpdf(j,i)=fpdf(j,i-1)*(1-qki(:,i-1))+fpdf(j-1,i-1)*qki(:,i-1);
     end
   fpdf(i,i)=fpdf(i-1,i-1)*qki(:,i-1);
  end
out=fpdf(:,end);
end% of function TC_gamma


function qki= Vec_CondPTC_gamma(x_norm,x_gamma,h,rho,T,v) 
PD = (1-exp(-kron(h,T)));DB = tinv(PD,v);

a=rho.^0.5; sqrt1_a2 = sqrt(1-sum(a.*a,2));

aM = gtimes(a, x_norm'); Sqrt_W=gamcdf(x_gamma,v/2,2).^0.5;

DB_times_W= gtimes(DB,Sqrt_W); DB_minus_aM = gminus(DB_times_W',aM);

qki=normcdf(grdivide(DB_minus_aM,sqrt1_a2));
end% of function Vec_CondPTC
RVNorman
  • 123
  • 1
  • 8
  • Your code is hard to read; can you please indent it? Also, why can you not just use the `quad` function in matlab? – Isaac Jul 25 '12 at 19:59
  • Hi Issac,apologies indents added to improve readability; Three reasons 1. I ve not seen/used the inbuilt functions for integration of two differing densities (N(0,1) & Chi-sq here), 2 quad cant easily handle infinite limits like this case; 3 dblquad - needs areas to be mapped to a rectangle, not sure how to do this here. Open to suggestions. – RVNorman Jul 25 '12 at 20:54
  • The approach seems sound mathematically, I don't know enough matlab to help, but you'll have 2 nested loops, and you just have to multiply the weights. – Alexandre C. Jul 25 '12 at 21:26
  • Hi Alexandre, could you indicate how you would do this in C or someother language? It may be an issue with my scaling of the weights but simple element by element multiplication doesnt give the required results. – RVNorman Aug 03 '12 at 11:16
  • LaTeX markup seems broken; I tried to render it over on [math.se] (to insert screenshot here), but could not figure out what you mean by `\sqrt{1a {i^2}` –  Mar 02 '15 at 19:58
  • Not sure whats happening with Latex markup. In the code ai=rho.^0.5 and ai2=(1-ai.*ai).^0.5. Its the square root of the correlation to the common factor eqns 5 and 6 in this [link](https://www.business.unsw.edu.au/About-Site/Schools-Site/risk-actuarial-site/Documents/Y.%20Goegebeur,%20T.%20Hoedemakers%20and%20J.%20Tistaert%20-%20Synthetic%20CDO%20Pricing%20using%20the%20Student%20t%20Factor%20Model%20with%20Random%20Recovery.pdf) – RVNorman Mar 03 '15 at 20:35

0 Answers0