mata:

mata clear

real scalar integrate(real scalar lower, real scalar upper, real vector arg) {
  if (((lower == .) & (upper == .)) | ((lower == 0) & (upper == .)) |
        ((lower != .) & (upper != .))) {
    return(integrate_main(lower, upper, arg))
  }
  else if ((lower == .) & (upper != .)) {
    return(integrate_main(0, upper, arg) + integrate_main(0, ., arg))
  }
  else if ((lower != 0) & (upper == .)) {
    return(integrate_main(lower, 0, arg) + integrate_main(0, ., arg))
  }
  else {
    return(integrate_main(lower, upper, arg))
  }  
}

real matrix integrate_main(real lower, real upper, real vector arg) {
  if ((lower != .) & (upper != .)) {
    rw  = ( .9997137268,  .9984919506,  .9962951347,  .9931249370,
	        .9889843952,  .9838775407,  .9778093585,  .9707857758,
			.9628136543,  .9539007829,  .9440558701,  .9332885350,
			.9216092981,  .9090295710,  .8955616450,  .8812186794,
			.8660146885,  .8499645279,  .8330838799,  .8153892383,
			.7968978924,  .7776279096,  .7575981185,  .7368280898,
			.7153381176,  .6931491994,  .6702830156,  .6467619085,
			.6226088602,  .5978474702,  .5725019326,  .5465970121,
			.5201580199,  .4932107892,  .4657816498,  .4378974022,
			.4095852917,  .3808729816,  .3517885264,  .3223603439,
			.2926171880,  .2625881204,  .2323024818,  .2017898641,
			.1710800805,  .1402031372,  .1091892036,  .0780685828,
			.0468716824,  .0156289844, -.0156289844, -.0468716824,
		   -.0780685828, -.1091892036, -.1402031372, -.1710800805,
		   -.2017898641, -.2323024818, -.2625881204, -.2926171880,
           -.3223603439, -.3517885264, -.3808729816, -.4095852917,
		   -.4378974022, -.4657816498, -.4932107892, -.5201580199,
		   -.5465970121, -.5725019326, -.5978474702, -.6226088602,
           -.6467619085, -.6702830156, -.6931491994, -.7153381176,
		   -.7368280898, -.7575981185, -.7776279096, -.7968978924,
		   -.8153892383, -.8330838799, -.8499645279, -.8660146885,
           -.8812186794, -.8955616450, -.9090295710, -.9216092981,
		   -.9332885350, -.9440558701, -.9539007829, -.9628136543,
		   -.9707857758, -.9778093585, -.9838775407, -.9889843952,
           -.9931249370, -.9962951347, -.9984919506, -.9997137268 \
			.0007346345,  .0017093927,  .0026839254,  .0036559612,
			.0046244501,  .0055884280,  .0065469485,  .0074990733,
            .0084438715,  .0093804197,  .0103078026,  .0112251140,
			.0121314577,  .0130259479,  .0139077107,  .0147758845,
			.0156296211,  .0164680862,  .0172904606,  .0180959407,
            .0188837396,  .0196530875,  .0204032326,  .0211334421,
			.0218430024,  .0225312203,  .0231974232,  .0238409603,
			.0244612027,  .0250575445,  .0256294029,  .0261762192,
            .0266974592,  .0271926134,  .0276611982,  .0281027557,
			.0285168543,  .0289030896,  .0292610841,  .0295904881,
			.0298909796,  .0301622651,  .0304040795,  .0306161866,
            .0307983790,  .0309504789,  .0310723374,  .0311638357,
			.0312248843,  .0312554235,  .0312554235,  .0312248843,
			.0311638357,  .0310723374,  .0309504789,  .0307983790,
            .0306161866,  .0304040795,  .0301622651,  .0298909796,
			.0295904881,  .0292610841,  .0289030896,  .0285168543,
			.0281027557,  .0276611982,  .0271926134,  .0266974592,
            .0261762192,  .0256294029,  .0250575445,  .0244612027,
			.0238409603,  .0231974232,  .0225312203,  .0218430024,
			.0211334421,  .0204032326,  .0196530875,  .0188837396,
            .0180959407,  .0172904606,  .0164680862,  .0156296211,
			.0147758845,  .0139077107,  .0130259479,  .0121314577,
			.0112251140,  .0103078026,  .0093804197,  .0084438715,
            .0074990733,  .0065469485,  .0055884280,  .0046244501,
			.0036559612,  .0026839254,  .0017093927,  .0007346345)
    sum = rw[2, ]:*sf((upper - lower)/2*rw[1, ] :+ (upper + lower)/2, arg)
    return((upper - lower)/2*quadrowsum(sum))
  }
  else if (lower == 0 & upper == .) {
    rw  = (374.9841128000, 355.2613119000, 339.4351019000, 325.6912634000,
	       313.3295340000, 301.9858553000, 291.4401336000, 281.5463283000,
	       272.2011700000, 263.3281685000, 254.8686293000, 246.7762410000, 
           239.0136298000, 231.5500680000, 224.3598948000, 217.4213933000,
		   210.7159729000, 204.2275596000, 197.9421331000, 191.8473694000,
		   185.9323602000, 180.1873909000, 174.6037612000, 169.1736398000,
           163.8899456000, 158.7462485000, 153.7366875000, 148.8559014000,
		   144.0989700000, 139.4613646000, 134.9389050000, 130.5277232000,
		   126.2242308000, 122.0250921000, 117.9271991000, 113.9276512000,
           110.0237356000, 106.2129115000, 102.4927949000,  98.8611460500,
		    95.3158573500,  91.8549432600,  88.4765308200,  85.1788512100,
 	 	    81.9602322200,  78.8190913200,  75.7539294700,  72.7633254300,
            69.8459306400,  67.0004645200,  64.2257101200,  61.5205102900,
		    58.8837639800,  56.3144229900,  53.8114889200,  51.3740103300,
		    49.0010802100,  46.6918335400,  44.4454451100,  42.2611274800,
            40.1381290600,  38.0757323700,  36.0732524100,  34.1300351200,
		    32.2454560000,  30.4189187700,  28.6498541500,  26.9377187300,
		    25.2819938700,  23.6821847600,  22.1378194500,  20.6484479700,
            19.2136415800,  17.8329919500,  16.5061104700,  15.2326276000,
		    14.0121922500,  12.8444711800,  11.7291484900,  10.6659250900,
		     9.6545182440,   8.6946611140,   7.7861023780,   6.9286058290,
             6.1219500310,   5.3659279860,   4.6603468360,   4.0050275820,
		     3.3998048270,   2.8445265430,   2.3390538500,   1.8832608260,
		     1.4770343300,   1.1202738350,    .8128912841,    .5548109376,
              .3459691810,    .1863141021,    .0758036120,    .0143861470 \
		     7.596410e-96,   3.146290e-94,   9.411750e-96,   7.09970e-106,
		     3.82990e-102,   2.953620e-97,   2.133670e-97,   2.228220e-97,
             1.326430e-95,   5.932280e-99,   1.291530e-99,   1.15800e-100,
		     1.657940e-95,   4.53820e-100,   1.815580e-90,   1.437780e-93,
		     2.020360e-91,   1.279760e-88,   6.704800e-86,   2.884880e-83,
             1.037900e-80,   3.152470e-78,   8.153750e-76,   1.809860e-73,
		     3.471850e-71,   5.792630e-69,   8.454960e-67,   1.085370e-64,
		     1.231380e-62,   1.240240e-60,   1.113570e-58,   8.947340e-57,
             6.456270e-55,   4.197750e-53,   2.466810e-51,   1.313980e-49,
		     6.361270e-48,   2.806040e-46,   1.130480e-44,   4.168860e-43,
		     1.410140e-41,   4.383790e-40,   1.254830e-38,   3.313060e-37,
             8.081630e-36,   1.824200e-34,   3.815860e-33,   7.407430e-32,
		     1.336210e-30,   2.242650e-29,   3.506270e-28,   5.112370e-27,
		     6.959200e-26,   8.853230e-25,   1.053590e-23,   1.174020e-22,
             1.226010e-21,   1.200850e-20,   1.104090e-19,   9.536250e-19,
		     7.743090e-18,   5.914430e-17,   4.252610e-16,   2.880120e-15,
		     1.838350e-14,   1.106500e-13,   6.283520e-13,   3.368210e-12,
             1.705060e-11,   8.154800e-11,   3.686330e-10,   1.575600e-09,
		     6.369710e-09,   2.436430e-08,   8.820060e-08,   3.022640e-07,
		     9.808340e-07,   3.014270e-06,   8.774310e-06,    .0000241958,
              .0000632109,    .0001564521,    .0003668548,    .0008148716,
		      .0017143197,    .0034149800,    .0064389510,    .0114854424,
		      .0193678281,    .0308463086,    .0463401336,    .0655510093,
              .0870966385,    .1083141121,    .1254070908,    .1340433397,
		      .1303566130,    .1121151033,    .0796767462,    .0363926059)
    sum = rw[2, ]:*exp(rw[1, ]):*sf(rw[1, ], arg)
    return(quadrowsum(sum))
  }
  else if (lower == . & upper == .) {
    rw  = ( 13.4064873400,  12.82379975000,  12.3429642200,  11.9150619400,
	        11.5214154000,  11.15240439000,  10.8022607500,  10.4671854200,
		    10.1445099400,   9.83226980800,   9.5289658230,   9.2334208900,
             8.9446892170,   8.66199616800,   8.3846969400,   8.1122473110,
		     7.8441823840,   7.58010080800,   7.3196528220,   7.0625310600,
		     6.8084633530,   6.55720703200,   6.3085443610,   6.0622788330,
             5.8182321350,   5.57624164900,   5.3361583600,   5.0978451050,
		     4.8611750920,   4.62603063600,   4.3923020790,   4.1598868550,
		     3.9286886830,   3.69861685900,   3.4695856360,   3.2415136800,
             3.0143235800,   2.78794142400,   2.5622964020,   2.3373204640,
		     2.1129479960,   1.88911553700,   1.6657615090,   1.4428259700,
		     1.2202503910,    .99797743610,    .7759507615,    .5541148236,
              .3324146923,    .11079587240,   -.1107958724,   -.3324146923,
		     -.5541148236,   -.77595076150,   -.9979774361,  -1.2202503910,
		    -1.4428259700,   -1.6657615090,  -1.8891155370,  -2.1129479960,
            -2.3373204640,   -2.5622964020,  -2.7879414240,  -3.0143235800,
		    -3.2415136800,   -3.4695856360,  -3.6986168590,  -3.9286886830,
		    -4.1598868550,   -4.3923020790,  -4.6260306360,  -4.8611750920,
            -5.0978451050,   -5.3361583600,  -5.5762416490,  -5.8182321350,
		    -6.0622788330,   -6.3085443610,  -6.5572070320,  -6.8084633530,
		    -7.0625310600,   -7.3196528220,  -7.5801008080,  -7.8441823840,
            -8.1122473110,   -8.3846969400,  -8.6619961680,  -8.9446892170,
		    -9.2334208900,   -9.5289658230,  -9.8322698080, -10.1445099400,
		   -10.4671854200,  -10.8022607500, -11.1524043900, -11.5214154000,
           -11.9150619400,  -12.3429642200, -12.8237997500, -13.4064873400 \
		     5.908070e-79,    1.972860e-72,   3.083030e-67,   9.019220e-63,
		     8.518880e-59,    3.459480e-55,   7.191530e-52,   8.597560e-49,
             6.420730e-46,    3.185220e-43,   1.100470e-40,   2.748780e-38,
		     5.116230e-36,    7.274570e-34,   8.067430e-32,   7.101810e-30,
		     5.037790e-28,    2.917350e-26,   1.394840e-24,   5.561030e-23,
             1.865000e-21,    5.302320e-20,   1.286830e-18,   2.682490e-17,
		     4.829840e-16,    7.548900e-15,   1.028870e-13,   1.227880e-12,
		     1.287900e-11,    1.191300e-10,   9.747920e-10,   7.075860e-09,
             4.568130e-08,    2.629100e-07,   1.351800e-06,   6.221520e-06,
		      .0000256762,     .0000951716,    .0003172920,    .0009526922,
		      .0025792733,     .0063030003,    .0139156652,    .0277791274,
              .0501758127,     .0820518274,    .1215379868,    .1631300305,
		      .1984628503,     .2188926296,    .2188926296,    .1984628503,
		      .1631300305,     .1215379868,    .0820518274,    .0501758127,
              .0277791274,     .0139156652,    .0063030003,    .0025792733,
		      .0009526922,     .0003172920,    .0000951716,    .0000256762,
		     6.221520e-06,    1.351800e-06,   2.629100e-07,   4.568130e-08,
             7.075860e-09,    9.747920e-10,   1.191300e-10,   1.287900e-11,
		     1.227880e-12,    1.028870e-13,   7.548900e-15,   4.829840e-16,
		     2.682490e-17,    1.286830e-18,   5.302320e-20,   1.865000e-21,
             5.561030e-23,    1.394840e-24,   2.917350e-26,   5.037790e-28,
		     7.101810e-30,    8.067430e-32,   7.274570e-34,   5.116230e-36,
		     2.748780e-38,    1.100470e-40,   3.185220e-43,   6.420730e-46,
             8.597560e-49,    7.191530e-52,   3.459480e-55,   8.518880e-59,
		     9.019220e-63,    3.083030e-67,   1.972860e-72,   5.908070e-79)
    sum = rw[2, ]:*exp(rw[1, ]:^2):*sf(rw[1, ], arg)
    return(quadrowsum(sum))
  }
}

real colvector invmvnormal_mata(real scalar p, real vector mean,
                                real matrix Sigma, string tail,
								real scalar max_iter, real scalar tolerance,
								string integrator, real scalar shifts,
								real scalar samples) {		
					
  if ((p == 0) | (p == 1)) {
    return((. \ 0 \ 0 \ 0 \ 0))
  }
  else {
    if (rows(Sigma) == 1) {
      if (tail == "upper") {
        p          = 1 - p
      }
      else if (tail == "both") {
        p          = 0.5 + 0.5*p
      }
      return((invnormal(p) + mean \ 0 \ 0 \ 0 \ 0))
    }
    else {
	  if (tail == "both") {
	    a          = 10^-6
	  }
	  else {
	    a          = -10
	  }
	  b            = 10
	  fa           = invmvnormal_mata_int(a, p, mean, Sigma, tail, integrator,
	                                      shifts, samples)
	  if (fa == 0) {
	    return((a \ 0 \ 0 \ 0 \ 0))
	  }
	  fb           = invmvnormal_mata_int(b, p, mean, Sigma, tail, integrator,
	                                      shifts, samples)
	  if (fb == 0) {
	    return((b \ 0 \ 0 \ 0 \ 0))
	  }
	  if (tail == "both") {
	    while ((fa > 0) & (a >= 10^-20)) {
	      a        = a/2
	      fa       = invmvnormal_mata_int(a, p, mean, Sigma, tail, integrator,
	                                      shifts, samples)
	    }
	    while ((fb < 0) & (b <= 10^6)) {
	      b        = 2*b
	      fb       = invmvnormal_mata_int(b, p, mean, Sigma, tail, integrator,
	                                      shifts, samples)
	    }
	  }
	  else {
	    while ((((fa > 0) & (fb > 0)) | ((fa < 0) & (fb < 0))) & (a >= -10^6)) {
	      a        = 2*a
	      b        = 2*b
	      fa       = invmvnormal_mata_int(a, p, mean, Sigma, tail, integrator,
	                                      shifts, samples)
	      fb       = invmvnormal_mata_int(b, p, mean, Sigma, tail, integrator,
	                                      shifts, samples)
	    }
	  }
	  if (fa == 0) {
	    return((a \ 0 \ 0 \ 0 \ 0))
	  }
	  else if (fb == 0) {
	    return((b \ 0 \ 0 \ 0 \ 0))
	  }
	  else if (((fa < 0) & (fb > 0)) | ((fa > 0) & (fb < 0))) {
        half_tol   = 0.5*tolerance
		c          = b
        fc         = fb
        for (iter = 1; iter <= max_iter; iter++) {
          if (((fb > 0) & (fc > 0)) | ((fb < 0) & (fc < 0))) {
            c      = a
	        fc     = fa
	        d      = b - a
	        e      = d
	      }
	      if (abs(fc) < abs(fb)) {
	        a      = b
	        b      = c
	        c      = a
	        fa     = fb
	        fb     = fc
	        fc     = fa
	      }
	      tol1     = 6e-8*abs(b) + half_tol
	      xm       = 0.5*(c - b)
		  if ((abs(xm) <= tol1) | (fb == 0)) {
	        return((b \ abs(0.5*(b - a)) \ 0 \ fb \ iter))
	      }
	      if ((abs(e) >= tol1) & (abs(fa) > abs(fb))) {
	        s      = fb/fa
	        if (a == c) {
	          pi   = 2*xm*s
		      q    = 1 - s
	        }
	        else {
	          q    = fa/fc
		      r    = fb/fc
		      pi   = s*(2*xm*q*(q - r) - (b - a)*(r - 1))
              q    = (q - 1)*(r - 1)*(s - 1)
	        }
	        if (pi > 0) {
	          q    = -q
	        }
	        pi     = abs(pi)
	        if (2*pi < min((3*xm*q - abs(tol1*q), abs(e*q)))) {
	          e    = d
		      d    = pi/q
	        }
	        else {
	          e    = d = xm
	        }
	      }
	      else {
	        e      = d = xm
	      }
	      a        = b
	      fa       = fb
		  if (abs(d) > tol1) {
	        b      = b + d
	      }
	      else {
	        b      = b + tol1*sign(xm)
	      }
	      fb       = invmvnormal_mata_int(b, p, mean, Sigma, tail, integrator,
	                                      shifts, samples)
        }
		return((b \ abs(0.5*(b - a)) \ 1 \ fb \ iter))
	  }
	  else {
	    return((. \ . \ 2 \ . \ .))
	  }
    }
  }
}

real scalar invmvnormal_mata_int(real scalar q, real scalar p, real vector mean,
                                 real matrix Sigma, string tail,
								 string integrator, real scalar shifts,
								 real scalar samples) {
  k          = rows(Sigma)
  if (tail == "lower") {
    a        = J(1, k, .)
    b        = J(1, k, q)
  }
  else if (tail == "upper") {
    a        = J(1, k, q)
    b        = J(1, k, .)
  }
  else {
    a        = J(1, k, -q)
    b        = J(1, k, q)
  }
  if (integrator != "mvnormal") {
    return(pmvnormal_mata(a, b, mean, Sigma, shifts, samples, 3)[1] - p)
  }
  else {
    for (i = 1; i <= k; i++) {
	  if (a[i] == .) {
	    a[i] = -8e307
	  }
	  if (b[i] == .) {
	    b[i] = 8e307
	  }
	}
	return(mvnormalcv(a, b, mean, vech(Sigma)') - p)
  }
}

real colvector invmvt_mata(real scalar p, real vector delta, real matrix Sigma,
                           real scalar df, string tail, real scalar max_iter,
						   real scalar tolerance, string integrator,
						   real scalar shifts, real scalar samples) {		
  if ((p == 0) | (p == 1)) {
    return((. \ 0 \ 0 \ 0 \ 0))
  }
  else {
    k              = rows(Sigma)
    if (k == 1) {
      if (tail == "upper") {
        p          = 1 - p
      }
      else if (tail == "both") {
        p          = 0.5 + 0.5*p
      }
	  if ((df == 0) | (df == .)) {
        return((invnormal(p) + delta \ 0 \ 0 \ 0 \ 0))
	  }
	  else {
	    return((invt(df, p) + delta \ 0 \ 0 \ 0 \ 0))
	  }
    }
    else {
	  if (tail == "both") {
	    a          = 10^-6
	  }
	  else {
	    a          = -10
	  }
	  b            = 10
	  fa           = invmvt_mata_int(a, p, delta, Sigma, df, tail, integrator,
	                                 shifts, samples)
	  if (fa == 0) {
	    return((a \ 0 \ 0 \ 0 \ 0))
	  }
	  fb           = invmvt_mata_int(b, p, delta, Sigma, df, tail, integrator,
	                                 shifts, samples)
	  if (fb == 0) {
	    return((b \ 0 \ 0 \ 0 \ 0))
	  }
	  if (tail == "both") {
	    while ((fa > 0) & (a >= 10^-20)) {
	      a        = a/2
	      fa       = invmvt_mata_int(a, p, delta, Sigma, df, tail, integrator,
	                                 shifts, samples)
	    }
	    while ((fb < 0) & (b <= 10^6)) {
	      b        = 2*b
	      fb       = invmvt_mata_int(b, p, delta, Sigma, df, tail, integrator,
	                                 shifts, samples)
	    }
	  }
	  else {
	    while ((((fa > 0) & (fb > 0)) | ((fa < 0) & (fb < 0))) & (a >= -10^6)) {
	      a        = 2*a
	      b        = 2*b
	      fa       = invmvt_mata_int(a, p, delta, Sigma, df, tail, integrator,
	                                 shifts, samples)
	      fb       = invmvt_mata_int(b, p, delta, Sigma, df, tail, integrator,
	                                 shifts, samples)
	    }
	  }
	  if (fa == 0) {
	    return((a \ 0 \ 0 \ 0 \ 0))
	  }
	  else if (fb == 0) {
	    return((b \ 0 \ 0 \ 0 \ 0))
	  }
	  else if (((fa < 0) & (fb > 0)) | ((fa > 0) & (fb < 0))) {
        half_tol   = 0.5*tolerance
		c          = b
        fc         = fb
        for (iter = 1; iter <= max_iter; iter++) {
		  if (((fb > 0) & (fc > 0)) | ((fb < 0) & (fc < 0))) {
            c      = a
	        fc     = fa
	        d      = b - a
	        e      = d
	      }
	      if (abs(fc) < abs(fb)) {
	        a      = b
	        b      = c
	        c      = a
	        fa     = fb
	        fb     = fc
	        fc     = fa
	      }
	      tol1     = 6e-8*abs(b) + half_tol
	      xm       = 0.5*(c - b)
		  if ((abs(xm) <= tol1) | (fb == 0)) {
	        return((b \ abs(0.5*(b - a)) \ 0 \ fb \ iter))
	      }
	      if ((abs(e) >= tol1) & (abs(fa) > abs(fb))) {
	        s      = fb/fa
	        if (a == c) {
	          pi   = 2*xm*s
		      q    = 1 - s
	        }
	        else {
	          q    = fa/fc
		      r    = fb/fc
		      pi   = s*(2*xm*q*(q - r) - (b - a)*(r - 1))
              q    = (q - 1)*(r - 1)*(s - 1)
	        }
	        if (pi > 0) {
	          q    = -q
	        }
	        pi     = abs(pi)
	        if (2*pi < min((3*xm*q - abs(tol1*q), abs(e*q)))) {
	          e    = d
		      d    = pi/q
	        }
	        else {
	          e    = d = xm
	        }
	      }
	      else {
	        e      = d = xm
	      }
	      a        = b
	      fa       = fb
		  if (abs(d) > tol1) {
	        b      = b + d
	      }
	      else {
	        b      = b + tol1*sign(xm)
	      }
	      fb       = invmvt_mata_int(b, p, delta, Sigma, df, tail, integrator,
	                                 shifts, samples)
        }
		return((b \ abs(0.5*(b - a)) \ 1 \ fb \ iter))
	  }
	  else {
	    return((. \ . \ 2 \ . \ .))
	  }
    }
  }
}

real scalar invmvt_mata_int(real scalar q, real scalar p, real vector delta,
                            real matrix Sigma, real scalar df, string tail,
							string integrator, real scalar shifts,
							real scalar samples) {
  k   = rows(Sigma)
  if (tail == "lower") {
    a = J(1, k, .)
    b = J(1, k, q)
  }
  else if (tail == "upper") {
    a = J(1, k, q)
    b = J(1, k, .)
  }
  else {
    a = J(1, k, -q)
    b = J(1, k, q)
  }
  return(mvt_mata(a, b, delta, Sigma, df, integrator, shifts, samples,
                  3)[1] - p)
}

real colvector invtmvnormal_mata(real scalar p, real vector lowert,
                                 real vector uppert, real vector mean,
                                 real matrix Sigma, string tail,
								 real scalar max_iter, real scalar tolerance,
								 string integrator, real scalar shifts,
								 real scalar samples) {		
					
  if (p == 0) {
    return((min(lowert) \ 0 \ 0 \ 0 \ 0))
  }
  else if (p == 1) {
    return((max(uppert) \ 0 \ 0 \ 0 \ 0))
  }
  else {
    k               = rows(Sigma)
	if (integrator != "mvnormal") {
	  denominator   = pmvnormal_mata(lowert, uppert, mean, Sigma, shifts,
	                                 samples, 3)[1]
	}
	else {
	  for (i = 1; i <= k; i++) {
	    if (lowert[i] == .) {
	      lowert[i] = -8e307
	    }
	    if (uppert[i] == .) {
	      uppert[i] = 8e307
	    }
	  }
	  denominator   = mvnormalcv(a, b, mean, vech(Sigma)')
	}
	a               = min(lowert)
	if (a == .) {
	  a             = -8e307
	}
	b               = max(uppert)
	if (b == .) {
	  b             = 8e307
	}
	fa              = invtmvnormal_mata_int(a, p, lowert, uppert, mean, Sigma,
	                                        tail, integrator, shifts, samples,
											denominator)
    fb              = invtmvnormal_mata_int(b, p, lowert, uppert, mean, Sigma,
	                                        tail, integrator, shifts, samples,
											denominator)
	if (((fa < 0) & (fb > 0)) | ((fa > 0) & (fb < 0))) {
      half_tol    = 0.5*tolerance
      c           = b
      fc          = fb
      for (iter = 1; iter <= max_iter; iter++) {
        if (((fb > 0) & (fc > 0)) | ((fb < 0) & (fc < 0))) {
          c       = a
	      fc      = fa
	      d       = b - a
	      e       = d
	    }
	    if (abs(fc) < abs(fb)) {
	      a       = b
	      b       = c
	      c       = a
	      fa      = fb
	      fb      = fc
	      fc      = fa
	    }
	    tol1      = 6e-8*abs(b) + half_tol
	    xm        = 0.5*(c - b)
		if ((abs(xm) <= tol1) | (fb == 0)) {
	      return((b \ abs(0.5*(b - a)) \ 0 \ fb \ iter))
	    }
	    if ((abs(e) >= tol1) & (abs(fa) > abs(fb))) {
	      s       = fb/fa
	      if (a == c) {
	        pi    = 2*xm*s
		    q     = 1 - s
	      }
	      else {
	        q     = fa/fc
		    r     = fb/fc
		    pi    = s*(2*xm*q*(q - r) - (b - a)*(r - 1))
            q     = (q - 1)*(r - 1)*(s - 1)
	      }
	      if (pi > 0) {
	        q     = -q
	      }
	      pi      = abs(pi)
	      if (2*pi < min((3*xm*q - abs(tol1*q), abs(e*q)))) {
	        e     = d
		    d     = pi/q
	      }
	      else {
	        e     = d = xm
	      }
	    }
	    else {
	      e       = d = xm
	    }
	    a         = b
	    fa        = fb
		if (abs(d) > tol1) {
	      b       = b + d
	    }
	    else {
	      b       = b + tol1*sign(xm)
	    }
	    fb        = invtmvnormal_mata_int(b, p, lowert, uppert, mean, Sigma,
		                                  tail, integrator, shifts, samples,
										  denominator)
      }
      return((b \ abs(0.5*(b - a)) \ 1 \ fb \ iter))
	}
	else {
	  return((. \ . \ 2 \ . \ .))
	}
  }
}

real scalar invtmvnormal_mata_int(real scalar q, real scalar p,
                                  real vector lowert, real vector uppert,
								  real vector mean, real matrix Sigma,
								  string tail, string integrator,
								  real scalar shifts, real scalar samples, 
								  real scalar denominator) {
  
  k             = rows(Sigma)
  if (tail == "lower") {
    a           = J(1, k, .)
    b           = J(1, k, q)
  }
  else if (tail == "upper") {
    a           = J(1, k, q)
    b           = J(1, k, .)
  }
  else {
    a           = J(1, k, -q)
    b           = J(1, k, q)
  }
  aact          = J(1, k, 0)
  bact          = J(1, k, 0)
  for (i = 1; i <= k; i++) {
    aact[i]     = max((a[i], lowert[i]))
	bact[i]     = min((b[i], uppert[i]))
  }
  if (integrator != "mvnormal") {
    return(pmvnormal_mata(aact, bact, mean, Sigma, shifts, samples,
	                      3)[1]/denominator - p)
  }
  else {
    for (i = 1; i <= k; i++) {
	  if (aact[i] == .) {
	    aact[i] = -8e307
	  }
	  if (b[i] == .) {
	    bact[i] = 8e307
	  }
	}
	return(mvnormalcv(aact, bact, mean, vech(Sigma)')/denominator - p)
  }
}

real colvector invtmvt_mata(real scalar p, real vector lowert,
                            real vector uppert, real vector delta,
                            real matrix Sigma, real scalar df, string tail,
							real scalar max_iter, real scalar tolerance,
							string integrator, real scalar shifts,
							real scalar samples) {		
					
  if (p == 0) {
    return((min(lowert) \ 0 \ 0 \ 0 \ 0))
  }
  else if (p == 1) {
    return((max(uppert) \ 0 \ 0 \ 0 \ 0))
  }
  else {
    k               = rows(Sigma)
	denominator     = mvt_mata(lowert, uppert, delta, Sigma, df, integrator,
	                           shifts, samples, 3)[1]
	a               = min(lowert)
	if (a == .) {
	  a             = -8e307
	}
	b               = max(uppert)
	if (b == .) {
	  b             = 8e307
	}
	fa              = invtmvt_mata_int(a, p, lowert, uppert, delta, Sigma, df,
	                                   tail, integrator, shifts, samples,
									   denominator)
    fb              = invtmvt_mata_int(b, p, lowert, uppert, delta, Sigma, df,
	                                   tail, integrator, shifts, samples,
									   denominator)
	if (((fa < 0) & (fb > 0)) | ((fa > 0) & (fb < 0))) {
      half_tol    = 0.5*tolerance
      c           = b
      fc          = fb
      for (iter = 1; iter <= max_iter; iter++) {
        if (((fb > 0) & (fc > 0)) | ((fb < 0) & (fc < 0))) {
          c       = a
	      fc      = fa
	      d       = b - a
	      e       = d
	    }
	    if (abs(fc) < abs(fb)) {
	      a       = b
	      b       = c
	      c       = a
	      fa      = fb
	      fb      = fc
	      fc      = fa
	    }
	    tol1      = 6e-8*abs(b) + half_tol
	    xm        = 0.5*(c - b)
		if ((abs(xm) <= tol1) | (fb == 0)) {
	      return((b \ abs(0.5*(b - a)) \ 0 \ fb \ iter))
	    }
	    if ((abs(e) >= tol1) & (abs(fa) > abs(fb))) {
	      s       = fb/fa
	      if (a == c) {
	        pi    = 2*xm*s
		    q     = 1 - s
	      }
	      else {
	        q     = fa/fc
		    r     = fb/fc
		    pi    = s*(2*xm*q*(q - r) - (b - a)*(r - 1))
            q     = (q - 1)*(r - 1)*(s - 1)
	      }
	      if (pi > 0) {
	        q     = -q
	      }
	      pi      = abs(pi)
	      if (2*pi < min((3*xm*q - abs(tol1*q), abs(e*q)))) {
	        e     = d
		    d     = pi/q
	      }
	      else {
	        e     = d = xm
	      }
	    }
	    else {
	      e       = d = xm
	    }
	    a         = b
	    fa        = fb
		if (abs(d) > tol1) {
	      b       = b + d
	    }
	    else {
	      b       = b + tol1*sign(xm)
	    }
	    fb        = invtmvt_mata_int(b, p, lowert, uppert, delta, Sigma, df,
	                                 tail, integrator, shifts, samples,
									 denominator)
      }
      return((b \ abs(0.5*(b - a)) \ 1 \ fb \ iter))
	}
	else {
	  return((. \ . \ 2 \ . \ .))
	}
  }
}

real scalar invtmvt_mata_int(real scalar q, real scalar p, real vector lowert,
                             real vector uppert, real vector delta,
                             real matrix Sigma, real scalar df, string tail,
							 string integrator, real scalar shifts,
							 real scalar samples, real scalar denominator) {
  k          = rows(Sigma)
  if (tail == "lower") {
    a        = J(1, k, .)
    b        = J(1, k, q)
  }
  else if (tail == "upper") {
    a        = J(1, k, q)
    b        = J(1, k, .)
  }
  else {
    a        = J(1, k, -q)
    b        = J(1, k, q)
  }
  aact       = J(1, k, 0)
  bact       = J(1, k, 0)
  for (i = 1; i <= k; i++) {
    aact[i]  = max((a[i], lowert[i]))
	bact[i]  = min((b[i], uppert[i]))
  }
  return(mvt_mata(aact, bact, delta, Sigma, df, integrator, shifts, samples,
                  3)[1]/denominator - p)
}

real colvector mvnormalden_mata(real vector x, real vector mean,
                                real matrix Sigma,| real matrix C) {
  if (args() < 4) {
    C         = cholesky(Sigma)
  }
  log_density = -sum(log(diagonal(C))) -
                  0.5*(rows(C)*log(2*pi()) + colsum(lusolve(C, (x - mean)'):^2))
  return((exp(log_density) \ log_density))
}

real colvector mvt_mata(real vector lower, real vector upper, real vector delta,
                        real matrix Sigma, real scalar df, string integrator,
						real scalar shifts, real scalar samples,
						real scalar alpha) {
  k                          = rows(Sigma)
  if ((df == 0) | (df == .)) {
	if (integrator == "mvnormal") {
	  for (i = 1; i <= k; i++) {
	    if (lower[i] == .) {
	      lower[i]           = -8e307
	    }
	    if (upper[i] == .) {
	      upper[i]           = 8e307
	    }
	  }
	  return((mvnormalcv(lower, upper, delta, vech(Sigma)') \ .))
	}
	else {
	  return(pmvnormal_mata(lower, upper, delta, Sigma, shifts, samples, alpha))
	}
  }
  else {
    if (k == 1) {
      if ((lower == .) & (upper == .)) {
	    I                    = 1
	  }
	  else if ((lower != .) & (upper == .)) {
	    I                    = 1 - t(df, lower - delta)
	  }
	  else if ((lower == .) & (upper != .)) {
	    I                    = t(df, upper - delta)
	  }
	  else {
	    I                    = t(df, upper - delta) - t(df, lower - delta)
	  }
	  E                      = 0
    }
    else {
      a                      = lower - delta
      b                      = upper - delta
      C                      = J(k, k, 0)
      u                      = y = J(1, k - 1, 0)
      atilde                 = (a[1], J(1, k - 2, 0))
      btilde                 = (b[1], J(1, k - 2, 0))
	  sqrt_Sigma11           = sqrt(Sigma[1, 1])
      for (i = 1; i <= k - 1; i++) {
        args = J(1, k - i + 1, 0)
        for (j = 1; j <= k - i + 1; j++){
          s                  = j + i - 1
          if (i > 1) {
		    y2               = sum(y[1::(i - 1)]:^2)
			Cy               = sum(C[s, 1::(i - 1)]:*y[1::(i - 1)])
			C2               = sum(C[s, 1::(i - 1)]:^2)
            if ((a[s] != .) & (b[s] != .)) {
              args[j]        =
			    t(df + i - 1,
			      sqrt((df + i - 1)/(df + y2))*((b[s] - Cy)/
				                                  sqrt(Sigma[s, s] - C2))) -
				t(df + i - 1,
				  sqrt((df + i - 1)/(df + y2))*((a[s] - Cy)/
						                          sqrt(Sigma[s, s] - C2)))
            }
            else if ((a[s] == .) & (b[s] != .)) {
              args[j]        =
			    t(df + i - 1,
				  sqrt((df + i - 1)/(df + y2))*((b[s] - Cy)/
                                                  sqrt(Sigma[s, s] - C2)))
            }
            else if ((b[s] == .) & (a[s] != .)) {
              args[j]        = 1 -
			    t(df + i - 1,
				  sqrt((df + i - 1)/(df + y2))*((a[s] - Cy)/
                                                  sqrt(Sigma[s, s] - C2)))
            }
            else if ((a[s] == .) & (b[s] == .)) {
              args[j]        = 1
            }
          } 
          else {
            if ((a[s] != .) & (b[s] != .)) {
              args[j]        = t(df, b[s]/sqrt_Sigma11) -
			                     t(df, a[s]/sqrt_Sigma11)
            }
            else if ((a[s] == .) & (b[s] != .)) {
              args[j]        = t(df, b[s]/sqrt_Sigma11)
            }
            else if ((b[s] == .) & (a[s] != .)) {
              args[j]        = 1 - t(df, a[s]/sqrt_Sigma11)
            }
            else if ((a[s] == .) & (b[s] == .)) {
              args[j]        = 1
            }
          }
        }
		ii                   = ww =.
        minindex(args, 1, ii, ww)
        m                    = i - 1 + ii[1]
        if (i != m) {
	      tempa              = a
          tempb              = b
          tempa[i]           = a[m]
          tempa[m]           = a[i]
          a                  = tempa
          tempb[i]           = b[m]
          tempb[m]           = b[i]
          b                  = tempb
          tempSigma          = Sigma
          tempSigma[i, 1::k] = Sigma[m, 1::k]
          tempSigma[m, 1::k] = Sigma[i, 1::k]
          Sigma              = tempSigma
          Sigma[1::k, i]     = tempSigma[1::k, m]
          Sigma[1::k, m]     = tempSigma[1::k, i]
	    }
        if (i > 1) {
          if (i != m) {
		    tempC            = C
            tempC[i, 1::k]   = C[m, 1::k]
            tempC[m, 1::k]   = C[i, 1::k]
            C                = tempC
            C[1::k, i]       = tempC[1::k, m]
            C[1::k, m]       = tempC[1::k, i]
		  }
          C[i, i]            = sqrt(Sigma[i, i] - sum(C[i, 1::(i - 1)]:^2))
          for (s = i + 1; s <= k; s++){
            C[s, i]          =
			  (Sigma[s, i] - sum(C[i, 1::(i - 1)]:*C[s, 1::(i - 1)]))/C[i, i]
          }
		  Cy                 = sum(C[i, 1::(i - 1)]:*y[1::(i - 1)])
          atilde[i]          = (a[i] - Cy)/C[i, i]
          btilde[i]          = (b[i] - Cy)/C[i, i]
        } else {
          C[i, i]            = sqrt(Sigma[i, i])
          C[2::k, i]         = Sigma[2::k, i]/C[i, i]
        }
        arg                  = (df, i)
        u[i]                 =
		  (gamma((df + 1)/2)/(gamma((df + i - 1)/2)*((df + i - 1)*pi())^0.5))*
            integrate(atilde[i], btilde[i], arg)/
                     (t(df + i - 1, btilde[i]) - t(df + i - 1, atilde[i]))
        if (i == 1) {
          y[i]               = u[i]
        }
        else {
          y[i]               =
		    u[i]*sqrt((df + sum(y[1::(i - 1)]:^2))/(df + i - 1))
        }
      }
      C[k, k]                = sqrt(Sigma[k, k] - sum(C[k, 1::(k - 1)]:^2))
	  sqrt_primes            = sqrt((2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37,
	                               41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83,
								   89, 97, 101, 103, 107, 109, 113, 127, 131,
								   137, 139, 149, 151, 157, 163, 167, 173, 179,
								   181, 191, 193, 197, 199, 211, 223, 227, 229,
								   233, 239, 241, 251, 257, 263, 269, 271, 277,
								   281, 283, 293, 307, 311, 313, 317, 331, 337,
								   347, 349, 353, 359, 367, 373, 379, 383, 389,
								   397, 401, 409, 419, 421, 431, 433, 439, 443,
							       449, 457, 461, 463, 467, 479, 487, 491, 499,
								   503, 509, 521, 523, 541)[1::(k - 1)])
	  I                      = V = 0
      d                      = J(1, k, 0)
      e                      = J(1, k, 1)
      f                      = J(1, k, 0)
      if (a[1] != .) {
         d[1]                = t(df, a[1]/C[1, 1])
      }
      if (b[1] != .) {
         e[1]                = t(df, b[1]/C[1, 1])
      }
      f                      = (e[1] - d[1], J(1, k - 1, 0))
	  for (i = 1; i <= shifts; i++) {
		Ii                   = 0
        Delta                = runiform(1, k - 1)
        for (j = 1; j <= samples; j++) {
		  u[1]               =
		    invt(df, d[1] + abs(2*mod(j*sqrt_primes[1] + Delta[1], 1) - 1)*
			                  (e[1] - d[1])) 
          y[1]               = u[1]
          for (l = 2; l <= k; l++) {
		    if ((a[l] != .) & (b[l] != .)) {
			  y2             = sum(y[1::(l - 1)]:^2)
			  Cy             = sum(C[l, 1::(l - 1)]:*y[1::(l - 1)])
		      d[l]           = t(df + l - 1,
			                     sqrt((df + l - 1)/(df + y2))*
                                   ((a[l] - Cy)/C[l, l]))
			  e[l]           = t(df + l - 1,
			                     sqrt((df + l - 1)/(df + y2))*
                                   ((b[l] - Cy)/C[l, l]))
		    }
		    else if ((a[l] != .) & (b[l] == .)) {
              d[l]           =
			    t(df + l - 1,
				  sqrt((df + l - 1)/(df + sum(y[1::(l - 1)]:^2)))*
                    ((a[l] - sum(C[l, 1::(l - 1)]:*y[1::(l - 1)]))/C[l, l]))
		    }
            else if ((a[l] == .) & (b[l] != .)) {
              e[l]           =
			    t(df + l - 1,
				  sqrt((df + l - 1)/(df + sum(y[1::(l - 1)]:^2)))*
                    ((b[l] - sum(C[l, 1::(l - 1)]:*y[1::(l - 1)]))/C[l, l]))
            }
		    f[l]             = (e[l] - d[l])*f[l - 1]
            if (l < k) {
              u[l]           =
			    invt(df + l - 1,
				     d[l] + abs(2*mod(j*sqrt_primes[l] + Delta[l], 1) - 1)*
				              (e[l] - d[l]))
              y[l]           =
			    u[l]*sqrt((df + sum(y[1::(l - 1)]:^2))/(df + l - 1))
            }
          }
          Ii                 = Ii + (f[k] - Ii)/j
        }
        del                  = (Ii - I)/i
        I                    = I + del
        V                    = (i - 2)*V/i + del^2
        E                    = alpha*sqrt(V)
      }
    }
  }
  return((I \ E))
}

real colvector mvtden_mata(real vector x, real vector delta, real matrix Sigma,
                           real scalar df,| real matrix C) {
  if (args() < 5) {
    C           = cholesky(Sigma)
  }
  if ((df == 0) | (df == .)) {
    log_density = -sum(log(diagonal(C))) -
                    0.5*(rows(C)*log(2*pi()) +
					       colsum(lusolve(C, (x - delta)'):^2))
  }
  else {
    k           = rows(C)
    log_density = lngamma((df + k)/2) -
	                (lngamma(df/2) + sum(log(diagonal(C))) +
                      (k/2)*log(pi()*df)) -
					0.5*(df + k)*
					  log(1 + (colsum(lusolve(C, (x - delta)'):^2)/df))
  }
  return((exp(log_density) \ log_density))
}

real colvector pmvnormal_mata(real vector lower, real vector upper, 
                              real vector mean, real matrix Sigma,
							  real scalar shifts, real scalar samples,
							  real scalar alpha) {
  k                             = rows(Sigma)
  if (k == 1) {
    if ((lower == .) & (upper == .)) {
	  I                         = 1
	}
	else if ((lower != .) & (upper == .)) {
	  I                         = 1 - normal((lower - mean)/sqrt(Sigma))
	}
	else if ((lower == .) & (upper != .)) {
	  I                         = normal((upper - mean)/sqrt(Sigma))
	}
	else {
	  sqrt_Sigma                = sqrt(Sigma)
      I                         = normal((upper - mean)/sqrt_Sigma) -
	                                normal((lower - mean)/sqrt_Sigma)
	}
	E                           = 0
  }
  else if (k == 2) {
    if (lower[1] == .) lower[1] = -8e307
	if (lower[2] == .) lower[2] = -8e307
	if (upper[1] == .) upper[1] = 8e307
	if (upper[2] == .) upper[2] = 8e307
	sqrt_Sigma                  = sqrt((Sigma[1, 1], Sigma[2, 2]))
	a                           = (lower - mean):/(sqrt_Sigma[1], sqrt_Sigma[2])
	b                           = (upper - mean):/(sqrt_Sigma[1], sqrt_Sigma[2])
	r                           = Sigma[1, 2]/(sqrt_Sigma[1]*sqrt_Sigma[2])
	I                           =
	  binormal(b[1], b[2], r) + binormal(a[1], a[2], r) -
	    binormal(b[1], a[2], r) - binormal(a[1], b[2], r)
	E                           = 0
  }
  else {
    a                           = lower - mean
    b                           = upper - mean
	C                           = J(k, k, 0)
	zero_k_min_2                = J(1, k - 2, 0)
	zero_k_min_1                = (zero_k_min_2, 0)
    y                           = zero_k_min_1
    atilde                      = (a[1], zero_k_min_2)
    btilde                      = (b[1], zero_k_min_2)
	sqrt_Sigma11                = sqrt(Sigma[1, 1])
    args                        = J(1, k, 1)
	for (j = 1; j <= k; j++) {
      if ((a[j] != .) & (b[j] != .)) {
        args[j]                 = normal(b[j]/sqrt_Sigma11) -
			                        normal(a[j]/sqrt_Sigma11) 
      }
      else if ((a[j] == .) & (b[j] != .)) {
        args[j]                 = normal(b[j]/sqrt_Sigma11)
      }
      else if ((b[j] == .) & (a[j] != .)) {
        args[j]                 = 1 - normal(a[j]/sqrt_Sigma11)
      }
    }
	ii                          = ww = .
    minindex(args, 1, ii, ww)
    if (ii[1] != 1) {
	  tempa                     = a
      tempb                     = b
      tempa[1]                  = a[ii[1]]
      tempa[ii[1]]              = a[1]
      a                         = tempa
      tempb[1]                  = b[ii[1]]
      tempb[ii[1]]              = b[1]
      b                         = tempb
      tempSigma                 = Sigma
      tempSigma[1, ]            = Sigma[ii[1], ]
      tempSigma[ii[1], ]        = Sigma[1, ]
      Sigma                     = tempSigma
      Sigma[, 1]                = tempSigma[, ii[1]]
      Sigma[, ii[1]]            = tempSigma[, 1]
	}
    C[, 1]                      = (sqrt_Sigma11 \ Sigma[2::k, 1]/sqrt_Sigma11)
    if (atilde[1] != btilde[1]) {
	  y[1]                      = (normalden(atilde[1]) - normalden(btilde[1]))/
                                    (normal(btilde[1]) - normal(atilde[1]))
	}
	for (i = 2; i <= k - 1; i++) {
      args                      = J(1, k - i + 1, 1)
	  i_vec                     = 1::(i - 1)
      for (j = 1; j <= k - i + 1; j++) {
        s                       = j + i - 1
		if ((a[s] != .) & (b[s] != .)) {
		    Cy                  = sum(C[s, i_vec]:*y[i_vec])
		    denom               = sqrt(Sigma[s, s] - sum(C[s, i_vec]:^2))
			args[j]             = normal((b[s] - Cy)/denom) -
                                    normal((a[s] - Cy)/denom)
          }
          else if ((a[s] == .) & (b[s] != .)) {
	        args[j]             =
			  normal((b[s] - sum(C[s, i_vec]:*y[i_vec]))/
			           sqrt(Sigma[s, s] - sum(C[s, i_vec]:^2))) 
          }
          else if ((b[s] == .) & (a[s] != .)) {
  	        args[j]             =
			  1 - normal((a[s] - sum(C[s, i_vec]:*y[i_vec]))/
			               sqrt(Sigma[s, s] - sum(C[s, i_vec]:^2)))
          }
      }
	  ii                        = ww = .
      minindex(args, 1, ii, ww)
      m                         = i - 1 + ii[1]
	  if (i != m) {
	    tempa                   = a
        tempb                   = b
        tempa[i]                = a[m]
        tempa[m]                = a[i]
        a                       = tempa
        tempb[i]                = b[m]
        tempb[m]                = b[i]
        b                       = tempb
        tempSigma               = Sigma
        tempSigma[i, ]          = Sigma[m, ]
        tempSigma[m, ]          = Sigma[i, ]
        Sigma                   = tempSigma
        Sigma[, i]              = tempSigma[, m]
        Sigma[, m]              = tempSigma[, i]
		tempC                   = C
        tempC[i, ]              = C[m, ]
        tempC[m, ]              = C[i, ]
        C                       = tempC
        C[, i]                  = tempC[, m]
        C[, m]                  = tempC[, i]
	  }
      C[i, i]                   = sqrt(Sigma[i, i] - sum(C[i, i_vec]:^2))
      i_vec2                    = (i + 1)::k
	  C[i_vec2, i]              =
	    (Sigma[i_vec2, i] - rowsum(J(k - i, 1, C[i, i_vec]):*C[i_vec2, i_vec]))/
		  C[i, i]
	  Cy                        = sum(C[i, i_vec]:*y[i_vec])
      atilde[i]                 = (a[i] - Cy)/C[i, i]
      btilde[i]                 = (b[i] - Cy)/C[i, i]
	  if (atilde[i] != btilde[i]) {
	    y[i]                    = (normalden(atilde[i]) - normalden(btilde[i]))/
                                    (normal(btilde[i]) - normal(atilde[i]))
	  }
    }
	C[k, k]                     = sqrt(Sigma[k, k] - sum(C[k, 1::(k - 1)]:^2))
    I                           = V = 0
    if (a[1] != .) {
      d                         = J(samples, 1, (normal(a[1]/C[1, 1]), zero_k_min_1))
    }
	else {
	  d                         = J(samples, 1, (zero_k_min_1, 0))
	}
    if (b[1] != .) {
      e                         = J(samples, 1, (normal(b[1]/C[1, 1]), J(1, k - 1, 1)))
    }
	else {
	  e                         = J(samples, 1, J(1, k, 1))
	}
    f                           = (e[, 1] - d[, 1], J(samples, 1, zero_k_min_1))
	y                           = J(samples, k - 1, 0)
	Delta                       = runiform(shifts, k - 1)
	samples_sqrt_primes         =
	  (1::samples)*sqrt((2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47,
	                     53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107,
						 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167,
						 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229,
						 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283,
						 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359,
						 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431,
						 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491,
						 499, 503, 509, 521, 523, 541)[1::(k - 1)])
	Ii                          = J(1, shifts, 0)
	for (i = 1; i <= shifts; i++) {
	  for (l = 2; l <= k; l++) {
		l_vec                 = 1::(l - 1)
		y[, l - 1]            =
		  invnormal(d[, l - 1] + abs(2*mod(samples_sqrt_primes[, l - 1] :+
       				                     Delta[i, l - 1], 1) :- 1):*
								   (e[, l - 1] - d[, l - 1]))
		if ((a[l] != .) & (b[l] != .)) {
		  Cy                  = rowsum(J(samples, 1, C[l, l_vec]):*y[, l_vec])
		  d[, l]                = normal((a[l] :- Cy)/C[l, l])
	      e[, l]                = normal((b[l] :- Cy)/C[l, l])
		  f[, l]                = (e[, l] :- d[, l]):*f[, l - 1]
		}
		else if ((a[l] != .) & (b[l] == .)) {
          d[, l]                =
	        normal((a[l] :- rowsum(J(samples, 1, C[l, l_vec]):*y[, l_vec]))/C[l, l])
		  f[, l]                = (1 :- d[, l]):*f[, l - 1]
		}
        else if ((a[l] == .) & (b[l] != .)) {
          e[, l]                =
	        normal((b[l] :- rowsum(J(samples, 1, C[l, l_vec]):*y[, l_vec]))/C[l, l])
		  f[, l]                = e[, l]:*f[, l - 1]
        }
		else {
		  f[, l]                = f[, l - 1]
		}
      }
	  for (j = 1; j <= samples; j++) {
	    Ii[i]                   = Ii[i] + (f[j, k] - Ii[i])/j
	  }
	  del                       = (Ii[i] - I)/i
	  I                         = I + del
	  V                         = (i - 2)*V/i + del^2
      E                         = alpha*sqrt(V)
	}
  }
  return((I \ E))
}

real matrix rmvnormal_mata(real scalar n, real vector mean,
                           real matrix Sigma,| real matrix C) {
  if (args() < 3) {
    C = cholesky(Sigma)
  }
  return((C*rnormal(rows(C), n, 0, 1))' + J(n, 1, mean))
}

real matrix rmvt_mata(real scalar n, real vector delta, real matrix Sigma,
                      real scalar df, | real matrix C) {
  if (args() < 5) {
    C = cholesky(Sigma)
  }
  if ((df == 0) | (df == .)) {
    return((C*rnormal(rows(C), n, 0, 1))' + J(n, 1, delta))
  }
  else {
    return((C*rnormal(rows(C), n, 0, 1))':/sqrt(rchi2(n, 1, df)/df) +
			 J(n, 1, delta))
  }
}

real matrix rtmvnormal_mata(real scalar n, real vector lowert,
                            real vector uppert, real vector mean,
							real matrix Sigma,| real matrix C) {
  if (args() < 6) {
    C              = cholesky(Sigma)
  }
  k                = rows(C)
  rmatrix          = J(n, k, 0)
  i                = 1
  while (i <= n) {
    gensample      = (C*rnormal(k, 1, 0, 1))' + mean
	check          = 1
	for (j = 1; j <= k; j++) {
	  if ((gensample[j] < lowert[j]) | (gensample[j] > uppert[j])) {
	    check      = 0
		break
	  }
	}
	if (check == 1) {
	  rmatrix[i, ] = gensample
	  i            = i + 1
	}
  }
  return(rmatrix)
}

real matrix rtmvt_mata(real scalar n, real vector lowert, real vector uppert,
                       real vector delta, real matrix Sigma,
					   real scalar df,| real matrix C) {
  if (args() < 7) {
    C              = cholesky(Sigma)
  }
  k                = rows(C)
  rmatrix          = J(n, k, 0)
  i                = 1
  while (i <= n) {
    if ((df > 0) & (df != .)) {
      gensample    = (C*rnormal(k, 1, 0, 1))':/sqrt(rchi2(1, 1, df)) + delta
	}
	else {
	  gensample    = (C*rnormal(k, 1, 0, 1))' + delta
	}
    check          = 1
	for (j = 1; j <= k; j++) {
	  if ((gensample[j] < lowert[j]) | (gensample[j] > uppert[j])) {
	    check      = 0
		break
	  }
	}
	if (check == 1) {
	  rmatrix[i, ] = gensample
	  i            = i + 1
	}
  }
  return(rmatrix)
}

real vector sf(real rowvector s, real vector arg) {
  return(s:*(1 :+ (s:^2)/(arg[1] + arg[2] - 1)):^(-(arg[1] + arg[2])/2))
}

real scalar tmvnormal_mata(real vector lower, real vector upper,
                           real vector lowert, real vector uppert,
						   real vector mean, real matrix Sigma,
						   string integrator, real scalar shifts,
						   real scalar samples) {
  k                 = rows(Sigma)
  loweract          = J(1, k, 0)
  upperact          = J(1, k, 0)
  for (i = 1; i <= k; i++) {
    loweract[i]     = max((lower[i], lowert[i]))
	upperact[i]     = min((upper[i], uppert[i]))
  }
  if (integrator == "mvnormal") {
    for (i = 1; i <= k; i++) {
	  if (loweract[i] == .) {
	    loweract[i] = -8e307
	  }
	  if (upperact[i] == .) {
	    upperact[i] = 8e307
	  }
	  if (lowert[i] == .) {
	    lowert[i]   = -8e307
	  }
	  if (uppert[i] == .) {
	    uppert[i]   = 8e307
	  }
	}
    I               = mvnormalcv(loweract, upperact, mean, vech(Sigma)')/
	                    mvnormalcv(lowert, uppert, mean, vech(Sigma)')
  } 	
  else {
    I               = pmvnormal_mata(loweract, upperact, mean, Sigma, shifts,
	                                 samples, 3)[1]/
	                    pmvnormal_mata(lowert, uppert, mean, Sigma, shifts, 
						               samples, 3)[1]
  }
  return(I)
}

real colvector tmvnormalden_mata(real vector x, real vector lowert,
                                 real vector uppert, real vector mean,
				   			     real matrix Sigma, string integrator,
							     real scalar shifts,
							     real scalar samples,| real matrix C) {
  if (args() < 9) {
    C               = cholesky(Sigma)
  }
  k                 = rows(Sigma)
  check             = 1
  for (i = 1; i <= k; i++) {
    if ((x[i] < lowert[i]) | (x[i] > uppert[i])) {
	  check         = 0
	  break
	}
  }
  if (check == 0) {
    return((0 \ .))
  }
  else {
    if (integrator == "pmvnormal") {
	  denominator   = pmvnormal_mata(lowert, uppert, mean, Sigma, shifts,
	                                 samples, 3)[1]
	}
	else {
	  for (i = 1; i <= k; i++) {
	    if (lowert[i] == .) {
		  lowert[i] = -8e307
		}
		if (uppert[i] == .) {
		  uppert[i] = 8e307
		}
	  }
	  denominator   = mvnormalcv(lowert, uppert, mean, vech(Sigma)')
	}
	log_density     = -sum(log(diagonal(C))) -
	                    0.5*(k*log(2*pi()) +
					           0.5*colsum(lusolve(C, (x - mean)'):^2)) -
					    log(denominator)
    return((exp(log_density) \ log_density))
  }
}

real scalar tmvt_mata(real vector lower, real vector upper, real vector lowert,
                      real vector uppert, real vector delta, real matrix Sigma,
					  real scalar df, string integrator, real scalar shifts,
					  real scalar samples) {
  k             = rows(Sigma)
  loweract      = J(1, k, 0)
  upperact      = J(1, k, 0)
  for (i = 1; i <= k; i++) {
    loweract[i] = max((lower[i], lowert[i]))
	upperact[i] = min((upper[i], uppert[i]))
  }
  I             = mvt_mata(lower, upper, delta, Sigma, df, integrator, shifts,
                           samples, 3)[1]/
		            mvt_mata(lowert, uppert, delta, Sigma, df, integrator,
					         shifts, samples, 3)[1]
  return(I)
}

real colvector tmvtden_mata(real vector x, real vector lowert,
                            real vector uppert, real vector delta,
							real matrix Sigma, real scalar df,
						    string integrator, real scalar shifts,
						    real scalar samples,| real matrix C) {
  if (args() < 10) {
    C             = cholesky(Sigma)
  }
  if ((df > 0) & (df != .)) {
    k             = rows(C)
    check         = 1
    for (i = 1; i <= k; i++) {
      if (x[i] < lowert[i] | x[i] > uppert[i]) {
	    check     = 0
	    break
	  }
    }
    if (check == 0) {
      return((0 \ .))
    }
    else {
	  denominator = mvt_mata(lowert, uppert, delta, Sigma, df, integrator,
	                         shifts, samples, 3)[1]
      log_density = lngamma((df + k)/2) -
	                  (lngamma(df/2) + sum(log(diagonal(C))) +
					    (k/2)*log(pi()*df)) -
					  0.5*(df + k)*
					  log(1 + (colsum(lusolve(C, (x - delta)'):^2)/df)) -
					  log(denominator[1])
	  return((exp(log_density) \ log_density))
    }
  }
  else {
    return(tmvnormalden_mata(x, lowert, uppert, delta, Sigma, integrator,
	                         shifts, samples, C))
  }
}

end