*! version 1.0 25jun2007, Henrik Stovring version 9.0 mata: // create integration nodes on interval (a; b) real matrix MCatnode(real matrix x, real matrix ab) { real scalar M real colvector delta real matrix xtr M = cols(x) delta = (ab[.,2] - ab[.,1]) :/ M xtr = J(rows(x), 2*M, .) for (j=1; j<=M; j++) { xtr[., j] = ab[.,1] + delta :* (j - 1) + x[.,j] :* delta xtr[., M+j] = ab[.,1] + delta :* (j - 1) + (- x[.,j] :+ 1) :* delta } return(xtr) } // create integration weights for MCatnode real matrix MCatwt(real matrix x, real matrix ab) { real scalar M M = cols(x) return((ab[.,2] - ab[.,1]) :/ (2 * M)) } // create integration points on interval (exp(-infty), exp(-a)), // ie. (0,exp(-a)) real matrix tlMCatnode(real matrix x, real colvector a) { real scalar M real matrix xtr M = cols(x) xtr = J(rows(x), 2*M, .) for (j=1; j<=M; j++) { xtr[., j] = a - log((x[.,j] :+ (j - 1)) / M) xtr[., M+j] = a - log((-x[.,j] :+ j) / M) } return(xtr) } // create integration weights for tlMCatnode real matrix tlMCatwt(real matrix x, | real colvector a) { pragma unused a real scalar M M = cols(x) uitr = J(rows(x), 2*M, .) for (j=1; j<=M; j++) { uitr[., j] = (x[.,j] :+ (j - 1)):^(-1) :/ 2 uitr[., M+j] = (-x[.,j] :+ j):^(-1) :/ 2 } return(uitr) } // create integration result from evaluated nodes and weights real colvector intres(real matrix fx, real matrix wt) { return(quadrowsum(fx :* wt)) } end