/// Derivative of cdf version 14 mata: mata clear mata: mata set matastrict on mata: real vector mvnxpb(real colvector b, real colvector mean, real matrix V, | real scalar dfdx) { real vector ret real vector x x = b - mean if (args()==3) ret = _mvnxpb(x, V) else if(args()==4) ret = partiald(x, V) return(ret) } /// multivariate normal CDF by genz biv.conditioning real scalar _mvnxpb(real colvector b, real matrix r) { // define variables real scalar n, ep, i, k, p, pb, ckk, dem, de, s, im, cii, u, w, v, am, bm, _ha, _hb, t real matrix c, D real colvector ap, bp, d, y, ai, bi, E, ab, bb, cb, sg, _hc1, _hc2, _hhc1, ina, inb n=rows(r) c = r //ap = a bp = b d = sqrt(diagonal(c)) for (i=1; i<=n; i++){ //ap[i,1] = ap[i,1]/d[i,1] bp[i,1] = bp[i,1]/d[i,1] } c = corr(c) y = J(n,1,0) p = 1 pb = 1 D = I(n) for (k=1; k<=n; k++){ im = k ckk = 0 dem = 1 s = 0 /// find minimum outer value for integration for (i=k; i<=n; i++){ cii = sqrt(c[i,i]) if (i>1){ if (k==1) s = s if (k>1) s = c[|i,1\i,(k-1)|]*y[|1\(k-1)|] } //ai = (ap[i]-s)/cii bi = (bp[i]-s)/cii de = normal(bi)-0 if (de<=dem){ ckk = cii dem = de //am = ai bm = bi im = i } } if (im>k){ //_ha = ap[im,1] //ap[im,1] = ap[k,1] //ap[k,1] = _ha _hb = bp[im,1] bp[im,1] = bp[k,1] bp[k,1] = _hb if (k>1){ _hc1 = c[|im,1\im,(k-1)|] c[|im,1\im,(k-1)|] = c[|k,1\k,(k-1)|] c[|k,1\k,(k-1)|] = _hc1 } if ((im+1)<=n){ _hc2 = c[|(im+1),im\n,im|] c[|(im+1),im\n,im|] = c[|(im+1),k\n,k|] c[|(im+1),k\n,k|] = _hc2 } if ((k+1)<=(im-1)){ t = c[|(k+1),k\(im-1),k|] c[|(k+1),k\(im-1),k|] = c[|im,(k+1)\im,(im-1)|]' c[|im,(k+1)\im,(im-1)|] = t' c[im,im] = c[k,k] } else if ((k+1)>(im-1)){ c[im,im] = c[k,k] } } if (ckk>0){ c[k,k] = ckk if ((k+1)<=n){ c[|k,(k+1)\k,n|] = J(1,(n-k),0) for (i=(k+1); i<=n; i++){ c[i,k] = c[i,k]/ckk c[|i,(k+1)\i,i|] = c[|i,(k+1)\i,i|] - c[i,k]*c[|(k+1),k\i,k|]' } } if (abs(dem)>0){ y[k] = (0-normalden(bm))/dem } } p = p*dem ////bivariate product if (mod(k,2)==0){ u = c[(k-1),(k-1)] v = c[k,k] w = c[k,(k-1)] c[|(k-1),(k-1)\n,k|] = c[|(k-1),(k-1)\n,k|]*(1/u, 0\-w/(u*v), 1/v) //ab = ap[|(k-1),1\k,1|] bb = bp[|(k-1),1\k,1|] cb = J(rows(bb),1,0) if (k>2){ cb = c[|(k-1),1\k,(k-2)|]*y[|1,1\(k-2),1|] } sg = (u^2, u*w\u*w, (w^2+v^2)) D[|(k-1),(k-1)\k,k|] = sg //ina = ab - cb inb = bb - cb E = bvnmmg(inb,sg) pb = pb*E[3,1] y[|(k-1),1\k,1|] = (E[1,1], E[2,1])' } } if (mod(n,2)==1) pb = pb*dem return(pb) /// function end } /// partial derivative of Phi(ss;vv,L) real colvector partiald(real colvector s, real matrix L) { real scalar n, i real colvector part, s_no_i, colll, v_i, at_vec real matrix L_no_i n = rows(L) part = J(n,1,.) for (i=1; i<=n; i++){ L_no_i = elimmat(L,i) s_no_i = elimvec(s,i) colll = collvec(L,i) v_i = colll*L[i,i]^(-1)*s[i] at_vec = s_no_i - v_i part[i] = normalden(s[i],0,sqrt(L[i,i]))*_mvnxpb(at_vec,L_no_i) } return(part) } /// Eliminate t-th row and column of matrix M real matrix elimmat(real matrix M, real scalar t) { real scalar n real matrix m1, m n = rows(M) if (t>1 & t1 & t1 & t