*! Date    : 17 Aug 2021
*! Version : 1.6
*! Authors : Michael J Grayling & Adrian P Mander

/*
  16/04/15 v1.0 Basic version complete.
  01/06/15 v1.1 Extended to use variable re-ordering.
  20/09/15 v1.2 Removed dependency on integrate() being pre-installed and
                tailored integration to our particular requirements.
  21/09/15 v1.3 Made small changes for improved efficiency.
  08/09/17 v1.4 Changed so df = 0 or df = . corresponds to mvnormal. Added
                support for non-integer df.
  03/01/19 v1.5 Minor changes for speed. Converted delta and sigma to be
                optional with internal defaults.
  17/08/21 v1.6	Fixed a small printing error.
*/

program define mvt, rclass
version 15.0
syntax , LOWer(numlist miss) UPPer(numlist miss) [DELta(numlist) ///
         Sigma(string) df(real 1) INTegrator(string) SHIfts(integer 12) ///
		 SAMples(integer 1000) ALPha(real 3)]

///// Check input variables ////////////////////////////////////////////////////

if ("`integrator'" == "") {
  local integrator "pmvnormal"
}
local lenlower:list sizeof lower
local lenupper:list sizeof upper
local lendelta:list sizeof delta
if (`lenlower' > 100) {
  di "{error}Only multivariate t distributions of dimension up to 100 are supported."
  exit(198)
}
if (`lenlower' != `lenupper') {
  di "{error}Vector of lower limits (lower) and vector of upper limits (upper) must be of equal length."
  exit(198)
}
forvalues i = 1/`lenlower' {
  local loweri:word `i' of `lower'
  local upperi:word `i' of `upper'
  if ((`loweri' != .) & (`upperi' != .)) {
    if (`loweri' >= `upperi') {
      di "{error}Each lower integration limit (in lower) must be strictly less than the corresponding upper limit (in upper)."
      exit(198)
    }
  }
}
if ((`lendelta' != 0) & (`lenlower' != `lendelta')) {
  di "{error}Vector of lower limits (lower) and vector of non-centrality parameters (delta) must be of equal length."
  exit(198)
}
if ("`sigma'" != "") {
  if (colsof(`sigma') != rowsof(`sigma')) {
    di "{error}Scale matrix Sigma (sigma) must be square."
    exit(198)
  }
  if (`lenlower' != colsof(`sigma')) {
    di "{error}Vector of lower limits (lower) must be the same length as the dimension of scale matrix Sigma (sigma)."
    exit(198)
  }
  cap mat chol = cholesky(`sigma')
  if (_rc > 0) {
    di "{error}Scale matrix Sigma (sigma) must be symmetric positive-definite."
    exit(198)
  }
  mat sigma    = (`sigma')
}
else {
  mat sigma    = I(`lenlower')
}
if (`df' < 0) {
  di "{error}Degrees of freedom (df) must be greater than or equal to zero."
  exit(198)
}
if (("`integrator'" != "mvnormal") & ("`integrator'" != "pmvnormal")) {
  di "{error}Choice of integrator (integrator) must be mvnormal or pmvnormal."
  exit(198)
}
if (`shifts' < 1) {
  di "{error}Number of shifts of the Quasi-Monte-Carlo integration algorithm to use (shifts) must be a strictly positive integer."
  exit(198)
}
if (`samples' < 1) {
  di "{error}Number of samples to use in each shift of the Quasi-Monte-Carlo integration algorithm (samples) must be a strictly positive integer."
  exit(198)
}
if (`alpha' <= 0) {
  di "{error}Chosen Monte-Carlo confidence factor (alpha) must be strictly positive."
  exit(198)
}

///// Perform main computations ////////////////////////////////////////////////

local matalower ""
foreach l of local lower {
  if "`matalower'" == "" local matalower "`l'"
  else local matalower "`matalower',`l'"
}
mat lower   = (`matalower')
local mataupper ""
foreach l of local upper {
  if "`mataupper'" == "" local mataupper "`l'"
  else local mataupper "`mataupper',`l'"
}
mat upper   = (`mataupper')
if (`lendelta' != 0) {
  local matadelta ""
  foreach l of local delta {
    if "`matadelta'" == "" local matadelta "`l'"
    else local matadelta "`matadelta',`l'"
  }
  mat delta = (`matadelta')
}
else {
  mat delta = J(1, `lenlower', 0)
}
mata: mvt_void(`df', "`integrator'", `shifts', `samples', `alpha')

///// Output ///////////////////////////////////////////////////////////////////

mat returns            = returns
return scalar error    = returns[2, 1]
return scalar integral = returns[1, 1]
di "{txt}integral = {res}" returns[1, 1] "{txt}"
di "{txt}error    = {res}" returns[2, 1] "{txt}"
end

///// Mata /////////////////////////////////////////////////////////////////////

mata:

void mvt_void(df, integrator, shifts, samples, alpha) {
  lower   = st_matrix("lower")
  upper   = st_matrix("upper")
  delta   = st_matrix("delta")
  Sigma   = st_matrix("sigma")
  st_matrix("returns", mvt_mata(lower, upper, delta, Sigma, df, integrator,
                                shifts, samples, alpha))
}

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 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) {
  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 {
    k                        = rows(Sigma)
    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 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, 1]                   = sqrt(Sigma[1, 1])
	  C[2::k, 1]                = Sigma[2::k, 1]/C[1, 1]
	}
	else {
	  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]), J(1, k - 1, 0)))
    }
	else {
	  d                         = J(samples, k, 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, k - 1, 0))
	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 vector sf(real rowvector s, real vector arg) {
  return(s:*(1 :+ (s:^2)/(arg[1] + arg[2] - 1)):^(-(arg[1] + arg[2])/2))
}

end