*! version 1.0.0  07jul2005 

mata

MVN_MAXDIM = 20
MVN_MAXPTS = 50

g_reps = J(MVN_MAXPTS,1,5)
g_p = J(MVN_MAXPTS,1,.)
g_eta = J(MVN_MAXPTS,1,NULL)
i = 0

g_p[++i] = 101
g_eta[i] = &(44,31,27,40,2,6,14,26,16,26,26,2,26,14,14,17,36,17,2)

g_p[++i] = 151
g_eta[i] = &(56,68,74,29,23,75,75,69,69,69,35,75,35,69,35,35,52,61,52)

g_p[++i] = 211
g_eta[i] = &(64,51,24,70,5,105,99,81,81,81,81,22,22,105,102,60,84,63,40)

g_p[++i] = 251
g_eta[i] = &(70,66,123,11,64,50,125,5,2,125,125,117,5,117,117,117,63,125,117)

g_p[++i] = 307
g_eta[i] = &(119,75,48,133,97,41,153,2,108,20,20,20,20,2,100,153,41,153,39)

g_p[++i] = 353
g_eta[i] = &(149,121,69,117,86,23,176,23,23,23,23,46,46,176,46,23,23,23,46)

g_p[++i] = 401
g_eta[i] = &(155,33,95,199,10,109,2,2,2,174,174,53,75,111,22,22,111,112,2)

g_p[++i] = 449
g_eta[i] = &(165,24,122,149,60,66,2,2,2,2,2,2,15,219,15,2,15,15,15)

g_p[++i] = 503
g_eta[i] = &(186,124,219,220,165,173,251,2,2,2,2,2,93,239,132,251,173,173,157)

g_p[++i] = 547
g_eta[i] = &(209,244,92,37,115,77,190,2,2,2,259,259,227,47,47,47,273,9,81)

g_p[++i] = 601
g_eta[i] = &(223,284,181,72,60,148,148,300,148,148,148,170,251,134,134,134,134,134,295)

g_p[++i] = 653
g_eta[i] = &(191,244,60,258,3,125,125,2,326,104,270,104,78,114,142,23,129,133,133)

g_p[++i] = 701
g_eta[i] = &(271,215,298,331,3,76,285,2,350,151,19,340,285,145,145,145,315,158,156)

g_p[++i] = 751
g_eta[i] = &(311,141,133,152,251,155,319,375,319,375,319,122,86,224,236,35,35,2,375)

g_p[++i] = 809
g_eta[i] = &(300,166,239,306,364,364,35,263,283,307,20,20,348,320,204,283,20,364,384)

g_p[++i] = 853
g_eta[i] = &(372,150,132,275,252,114,166,2,167,202,124,124,291,291,37,291,2,426,426)

g_p[++i] = 907
g_eta[i] = &(347,273,105,147,303,83,247,453,453,224,224,101,101,304,449,449,181,451,451)

g_p[++i] = 953
g_eta[i] = &(283,80,398,44,330,201,461,476,2,246,461,246,246,303,135,120,90,90,461)

g_p[++i] = 1009
g_eta[i] = &(390,445,207,334,337,291,477,2,504,2,504,2,2,503,337,337,226,55,2)

g_p[++i] = 1511
g_eta[i] = &(559,584,97,600,506,62,214,352,2,2,2,2,517,517,550,250,176,176,176)

g_p[++i] = 2003
g_eta[i] = &(830,421,854,410,18,881,92,762,2,1001,268,268,268,268,268,268,90,90,268)

g_p[++i] = 3001
g_eta[i] = &(1140,772,174,890,1459,1096,1085,627,627,2,2,245,245,245,49,49,1437,1437,5)

g_p[++i] = 4001
g_eta[i] = &(1478,722,113,62,132,1999,10,10,1016,200,200,20,2,2,508,759,445,787,788)

g_p[++i] = 5003
g_eta[i] = &(1850,1611,1659,1618,1173,513,1668,205,1611,2501,2501,2501,755,810,2238,2397,2397,457,104)

g_p[++i] = 6007
g_eta[i] = &(2488,2831,2610,2660,2006,2895,1314,827,1202,2004,2,1202,2004,2727,2727,2727,1885,814,2681)

g_p[++i] = 7001
g_eta[i] = &(1952,3061,2661,376,1583,2354,3378,2866,809,1813,3500,3375,15,885,885,1867,15,15,1867)

g_p[++i] = 8009
g_eta[i] = &(2430,3677,1982,2178,2050,2595,506,2803,2663,3927,2,4004,4004,4004,4004,887,887,887,887)

g_p[++i] = 9001
g_eta[i] = &(3798,3642,2003,3824,1714,1333,4499,4135,4183,4228,455,4500,2114,4228,1631,1631,2044,2044,1631)

g_p[++i] = 10007
g_eta[i] = &(4129,544,3227,4305,3720,4961,3335,5,4604,3683,3683,2,3683,2527,2527,4953,1286,1069,2)

g_p[++i] = 11003
g_eta[i] = &(4618,5197,3887,2747,1795,3760,2430,3,4610,3677,1179,2,1179,3677,2,2,3677,4669,4669)

g_p[++i] = 12007
g_eta[i] = &(4565,4564,5912,418,793,3987,4003,3621,4669,489,205,6003,2,3428,1533,2939,2465,5998,6003)

g_p[++i] = 13001
g_eta[i] = &(5071,5930,3065,686,4106,5101,6368,4334,3708,3708,3708,2,3708,3708,3708,3708,3708,3187,3187)

g_p[++i] = 14009
g_eta[i] = &(5351,1673,6187,1715,362,2887,1292,4642,1633,2297,205,7004,6902,6902,6902,205,6902,1189,1190)

g_p[++i] = 15013
g_eta[i] = &(6222,5619,4914,4035,690,216,607,4783,190,1977,1591,7506,7506,5883,231,231,4978,190,5785)

g_p[++i] = 16001
g_eta[i] = &(5911,4833,3365,5927,4632,4436,10,5334,40,3008,1008,20,20,20,8000,8000,8000,8000,126)

g_p[++i] = 17011
g_eta[i] = &(6511,1156,2212,6680,7062,3968,4693,3,6057,2968,402,2,2,2,2,8505,5262,3589,5262)

g_p[++i] = 18013
g_eta[i] = &(6611,6669,1805,2166,5688,8057,1065,4311,2400,5770,5012,2,2,2,6109,8219,8219,8900,6352)

g_p[++i] = 19001
g_eta[i] = &(7218,5875,7967,4308,6218,3343,7068,3,7827,5349,5229,6886,2,9500,4027,7619,3344,5947,5946)

g_p[++i] = 20011
g_eta[i] = &(7607,2759,7248,9432,4951,8593,181,6670,173,10,5064,5064,10005,6544,6544,6544,6544,6544,6544)

g_p[++i] = 25013
g_eta[i] = &(6974,9949,12352,4695,3910,10004,3080,9389,5083,7127,7852,7852,2,2,12506,7117,3233,3233,7117)

g_p[++i] = 30011
g_eta[i] = &(12574,2326,3544,13604,7616,5139,14021,15004,123,5460,8118,1116,2,15005,12599,3311,12599,12599,12599)

g_p[++i] = 35023
g_eta[i] = &(12957,2591,10158,8313,10441,7182,16930,11675,3,15087,16871,3907,17511,17511,2,2,4882,4882,165)

g_p[++i] = 40111
g_eta[i] = &(11080,5189,18518,1400,9834,17057,3971,14216,13370,11631,11631,11631,5426,2,2,16125,11631,2894,2894)

g_p[++i] = 45007
g_eta[i] = &(17381,21219,4250,20549,21444,16201,1868,9058,21422,11276,1856,1856,20803,22503,22503,22503,2,2,2)

g_p[++i] = 50021
g_eta[i] = &(19389,7086,23277,8315,2441,1163,20008,11363,16674,15162,17976,15778,15777,2,2,2,17976,15409,12881)

g_p[++i] = 60013
g_eta[i] = &(23210,16047,17895,17057,23969,1064,5569,29172,20004,10043,25857,20921,20921,2,2,25285,6368,6368,6368)

g_p[++i] = 70001
g_eta[i] = &(26747,11787,21248,8839,8487,6029,22882,16887,11667,5088,2849,10422,29567,2,35000,2,2,11734,11734)

g_p[++i] = 80021
g_eta[i] = &(33165,28441,37125,14123,10030,28239,35585,37801,30043,434,434,37828,7949,40010,40010,40010,25016,25016,13698)

g_p[++i] = 90001
g_eta[i] = &(34438,9150,11149,13128,28255,28715,27052,26857,44999,18460,33507,3674,3674,4404,45000,45000,45000,45000,45000)

g_p[++i] = 100003
g_eta[i] = &(38763,13758,31528,48778,15608,23828,2961,17206,50000,47636,26982,39335,44974,44974,50001,2,2,40739,17584)

if (i < MVN_MAXPTS) {
	errprintf("Note: resetting the maximum number of lattice generators from %g to %g\n", 
		MVN_MAXPTS, i)

	MVN_MAXPTS = i
	g_reps = g_reps[1::MVN_MAXPTS,1]
	g_p = g_p[1::MVN_MAXPTS,1]
	g_eta = g_eta[1::MVN_MAXPTS,1]
}
unlink("mvnlattice.mmat")
fh = fopen("mvnlattice.mmat", "w")
fputmatrix(fh, g_p)
fputmatrix(fh, g_eta)
fputmatrix(fh, g_reps)
fclose(fh)

mata drop g_p
mata drop g_eta
mata drop g_reps
mata drop i
mata drop fh
mata drop MVN_MAXDIM
mata drop MVN_MAXPTS


real rowvector mvnlattice(real vector a,    /* lower bound on integration 	*/
			  real vector b,    /* upper bound on integeration	*/
			  real matrix R,    /* covariance or correlation matrix	*/
			| real scalar maxn, /* max # of integration points	*/
			  real scalar eps)  /* desired error			*/
{
	real scalar p, n, fh
	real colvector a1, b1, ip
	real vector N, reps
	real matrix R1
	real rowvector result
	pointer(real vector) vector eta

	n = rows(R)
	if (cols(R) != n) {
		errprintf("R is not square\n")
		exit(3205)
	}
	if (missing(R)) {
		errprintf("R has missing values\n")
		exit(3351)
	}
	if (!issymmetric(R)) {
		errprintf("R is not symmetric; try R = (R+R')/2 to ensure symmetry\n")
		exit(3200)
	}
	if (length(a)!=n) {
		errprintf("lower bound vector, a, and R do not have the same dimension\n")
		exit(3200)
	}
	if (length(b)!=n) {
		errprintf("upper bound vector, b, and R do not have the same dimension\n")
		exit(3200)
	}
	if (n == 1) {
		real scalar s
	
		if (R <= 0) {
			errprintf("variance must be positive\n")
			exit(3300)
		}
		if (a > b) {
			errprintf("lower bound, a, must be greater than the upper bound, b\n")
			exit(3300)
		}
		s = sqrt(R)
		result = (normal(b/s)-normal(a/s),0,0)
	}
	else {
		stata("qui findfile mvnlattice.mmat")
		fh = fopen(st_global("r(fn)"), "r")
		N = fgetmatrix(fh)
		eta = fgetmatrix(fh)
		reps = fgetmatrix(fh)
		fclose(fh)
		
		pragma unset ip
		if (cols(a)==1) a1 = a
		else a1 = a'

		if (cols(b)==1) b1 = b
		else b1 = b'
		R1 = R
		p = mvnlattice_pivchol(a1,b1,R1,ip)
		if (missing(eps)) {
			eps = 1.0e-4
		}
		if (missing(maxn)) {
			maxn = 10000
		}
		if (p < 0) {
			result = mvnlattice_(a1, b1, R1, maxn, eps, N, eta, reps)	
		}
		else {
			result = J(1,3,.)
			result[1] = p
		}
	}
	return(result)
} 

/* pivot R on interval ranges (b, a) and factor */
real scalar mvnlattice_pivchol(real colvector a, real colvector b, real matrix R, 
	real colvector ip)
{
	real scalar  n
	real colvector d, inf, i, j
	real matrix D

	real scalar MVN_INFINITY

	MVN_INFINITY = 1000.0 

	/* assume dimensions are compatable */
	m = cols(R)

	ip = 1::m
	inf = J(m,1,0)
	d = J(m,1,0)

	mina = -MVN_INFINITY
	maxb = MVN_INFINITY

	eps = 100.0*st_numscalar("c(epsdouble)")
	s = sqrt(diagonal(R))
	if (missing(s) | sum(s==0.0)) {
		errprintf("R is not positive definite\n")
		exit(3353)
	}
	i = (b:>=.)
	mvnlattice_substitute(inf, i, 2)
	j = mvnlattice_select(ip,!i)
	b[j,1] = b[j]:/s[j]

	i = (a:>=.)
	j = mvnlattice_select(ip,i)
	inf[j,1] = inf[j] :+ 1
	j = mvnlattice_select(ip,!i)
	a[j,1] = a[j]:/s[j]

	i = (inf:==0)
	j = mvnlattice_select(ip,i)
	d[j,1] = b[j]-a[j]
	if (sum(d:<0)) {
		errprintf("lower bound, a, must be greater than the upper bound, b\n")
		exit(3300)
	}
	if (sum(d[j]:<=eps)) {
		/* interval width of zero; distribution is zero */
		return(0)
	}
	j = mvnlattice_select(ip,(inf:==1))
	d[j,1] = b[j]:-mina
	j = mvnlattice_select(ip,(inf:==2))
	d[j,1] = maxb:-a[j]
	mvnlattice_substitute(d, (inf:==3), 2*MVN_INFINITY)
	ip[.] = order(d,1)
	D = diag(1:/s)
	R = D*R*D

	n = sum(inf:<3)
	i = 1::n
	a = (a[ip])[i]
	b = (b[ip])[i]
	R = (R[ip,ip])[|1,1\n,n|]
	_cholesky(R)
	if (missing(R)) {
		errprintf("R is not positive definite\n")
		exit(3353)
	}
	return(-1)
}

real matrix mvnlattice_select(real matrix x, real vector s)
{
        real scalar i, j, k, l, n, c, r
        real vector res
                                     
	n = length(s)
	r = rows(x)
	c = cols(x)
	if (r*c != n) _error(3200)
	
	k = sum(s)
	res = J(k,1,.)
	if (k > 0) {
		l = 0
		for (k=1; k<=n; k++) {
			if (s[k] > 0) {
				i = floor((k-1)/c)+1
				j = k-(i-1)*c
				res[++l,1] = x[i,j]
			}
		}
	}
        return(res)
}

void mvnlattice_substitute(real matrix x, real vector s, real vector val)
{
        real scalar i, j, k, l, n, c, r, m
                                     
	n = length(s)
	r = rows(x)
	c = cols(x)
	if (r*c != n) _error(3200)
	
	k = sum(n)
	if (k == 0) return
	m = length(val)

	l = 0
	for (k=1; k<=n; k++) {
		if (s[k] > 0) {
			i = floor((k-1)/c)+1
			j = k-(i-1)*c
			l = mod(++l,m)+1
			x[i,j] = val[l]
		}
	}
}

void mvnlattice_transform(real scalar j, real colvector xj, real scalar aj0, 
	real scalar bj0, real matrix R, real matrix e, real colvector q)
{
	real scalar k, im
	real colvector aj, bj, na, nb

	im = 0
	n = length(xj)
	if (aj0 >= .) {
		im = 1
	}
	else {
		aj = J(n,1,aj0)
	}
	if (bj0 >= .) {
		im = im + 2
	}
	else {
		bj = J(n,1,bj0)
	}
	if (im == 1) {
		for (k=1; k<j; k++) {
			bj = bj :- R[j,k]:*e[.,k]
		}
		bj = bj:/R[j,j]
		nb = normal(bj)
		e[.,j] = invnormal(xj:*nb)
		q[.] = q:*nb
	}
	else if (im == 2) {
		for (k=1; k<j; k++) {
			aj = aj :- R[j,k]:*e[.,k]
		}
		aj = aj:/R[j,j]
		na = normal(aj)
		e[.,j] = invnormal((1:-xj):*na:+xj)
		q[.] = q:*(1:-na)
	}
	else if (im == 3) {
		q[.] = J(n,1,1)
	}
	else {
		for (k=1; k<j; k++) {
			aj = aj :- R[j,k]:*e[.,k]
			bj = bj :- R[j,k]:*e[.,k]
		}
		aj = aj:/R[j,j]
		bj = bj:/R[j,j]
		nb = normal(bj)
		na = normal(aj)
		e[.,j] = invnormal((1:-xj):*na:+xj:*nb)
		q[.] = q:*(nb:-na)
	}
}

real colvector mvnlattice_mod1(real colvector a) 
{
	return(a:-floor(a))
}

real rowvector mvnlattice_(real colvector a, real colvector b, real matrix R, 
	real scalar maxn, real scalar eps, real vector N, 
	pointer (real vector) vector eta, real vector reps)
{
	real scalar i, j, k, l, n, nl, nk, nn
	real scalar etak, p, pl, pk, d1, v, vl, s, aj, bj
	real colvector q1, q2, z, kv
	real rowvector u, z0
	real matrix e1, e2

	real scalar MVN_MAXPTS, MVN_BLOCKSZ

	MVN_BLOCKSZ = 50000
	MVN_MAXPTS = length(N)
	
	m = rows(R)
	z0 = J(m,1,0)
	u = J(m,1,0)
	z = J(0,1,0)
	e1 = e2 = J(0,m,0)
	q1 = q2 = J(0,1,0)
	
	p = v = 0
	s = 1
	n = 0
	i = 0 
	do {
		k = mod(i,MVN_MAXPTS)+1
		nl = N[k]
		etak = (*eta[k])[m-1]
		z0[1] = 1
	        for (j=2; j<=m; j++) {
			z0[j] = mod(z0[j-1]*etak,nl)
		}	
		pl = 0
		vl = 0
		nb = floor(nl/MVN_BLOCKSZ)+1

		for (l=1; l<=reps[k]; l++) { 
			/* use random shift for error analysis */
			/* Sec. 4.6 of Sloan and Joe (1994)    */
			u[.] = uniform(m,1) 
			k2 = 0
			pk = 0
			nn = 0
			
			for (ib=1; ib<=nb; ib++) {
				k1 = k2 + 1
				k2 = min((k2+MVN_BLOCKSZ,nl))
				kv = (k1::k2):/nl
				nk = length(kv)
				if (length(z) != nk) {
					e1 = e2 = J(nk,m,0)
					z = J(nk,1,0)
					q1 = q2 = J(nk,1,1)
				}
				else {
					q1[.] = q2[.] = J(nk,1,1)
				} 
				for (j=1; j<=m; j++) {
					/* random shift and periodizing transform (Genz, 2002) */
					z[.] = abs(2:*mvnlattice_mod1(kv:*z0[j]:+u[j]):-1)
					aj = a[j]
					bj = b[j]
					/* transform over a unit hypercube (Genz 1992) */
					mvnlattice_transform(j, z, aj, bj, R, e1, q1)
					z[.] = 1:-z
					mvnlattice_transform(j, z, aj, bj, R, e2, q2)
				}
				pk = pk+sum(q1)+sum(q2)
				nn = nn + 2*nk
			}
			pk = pk/nn
			pl = pl + (pk-pl)/l
			if (l > 1) {
				d1 = l-1
				vl = d1*vl/l + (pk-pl)^2/d1
			}
			n = n + 2*nl
			if (n>maxn && l>1) {
				break
			}
		}
		vl = vl/(reps[k]-1)/reps[k];
		if (++i == 1) {
			v = vl
			p = pl
		}
		else if (l > 1) {
			p = p + v*(pl-p)/(v+vl)
			v = v*vl/(v+vl)
		}
		s = 10*sqrt(v)
	/* Chebyshev inequality for error bound 	*/
	/* Pr(|eps| < k*s) >= 1-1/k^2, k=10, Pr>=.99	*/
	} while(n<=maxn && s>=eps)
	if (n>maxn & s>eps) {
		printf("{txt}note: estimated absolute error %g using %g points\n", s, n)
	}
	return((p,s,n))
}

end