nested quadrature in matlab using two quadrature rules
===

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

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