version 12.1 mata: mata clear real scalar mvnormalden_mata(real vector x, real vector mean, real matrix Sigma) { // Perform checks on input variables if (rows(Sigma) ~= cols(Sigma)) { "{error} Covariance matrix Sigma (Sigma) must be square." exit(198) } if (length(x) ~= length(mean)) { "{error} Vector of quantiles (x) and mean vector (mean) must be of equal length." exit(198) } if (length(mean) ~= cols(Sigma)) { "{error} Mean vector (mean) must be the same length as the dimension of covariance matrix Sigma (sigma)." exit(198) } // Initialise all required variables C = cholesky(Sigma) k = rows(C) // Compute density tmp = lusolve(C, (x :- mean)') rss = colsum(tmp:^2) logdensity = -sum(log(diagonal(C))) - 0.5*k*log(2*pi()) - 0.5*rss density = exp(logdensity) // Return result return(density) } real scalar mvtden_mata(real vector x, real vector delta, real matrix Sigma, real scalar df) { // Perform checks on input variables if (rows(Sigma) ~= cols(Sigma)) { "{error} Scale matrix Sigma (Sigma) must be square." exit(198) } if (length(x) ~= length(delta)) { "{error} Vector of quantiles (x) and vector of non-centrality parameters (delta) must be of equal length." exit(198) } if (length(delta) ~= cols(Sigma)) { "{error} Vector of non-centrality parameters (delta) must be the same length as the dimension of scale matrix Sigma (Sigma)." exit(198) } if ((df < 1) | (mod(df, 1) ~= 0)) { "{error} Degree of freedom (df) must be a positive integer." exit(198) } // Initialise all required variables C = cholesky(Sigma) k = rows(C) // Compute density tmp = lusolve(C, (x :- delta)') rss = colsum(tmp:^2) logdensity = lngamma((df + k)/2) - (lngamma(df/2) + sum(log(diagonal(C))) + (k/2)*log(pi()*df)) - 0.5*(df + k)*log(1 + (rss/df)) density = exp(logdensity) // Return result return(density) } real matrix rmvnormal_mata(real scalar n, real vector mean, real matrix Sigma) { // Perform checks on input variables if (rows(Sigma) ~= cols(Sigma)) { "{error} Covariance matrix Sigma (Sigma) must be square." exit(198) } if (length(mean) ~= cols(Sigma)) { "{error} Mean vector (mean) must be the same length as the dimension of covariance matrix Sigma (sigma)." exit(198) } if ((n < 1) | (mod(n, 1) ~= 0)) { "{error} Number of random vectors to generate (n) must be a strictly positive integer." exit(198) } // Initialise all required variables C = cholesky(Sigma) k = rows(C) // Initialise a matrix with each row the vector mean meanmat = J(n, k, 0) for (i = 1; i <= n; i++) { meanmat[i, 1::k] = mean } // Create a matrix of random N(0,1) variables z = J(k, n, 0) for (i = 1; i <= k; i++) { for (j = 1; j <= n; j++) { z[i, j] = rnormal(1, 1, 0, 1) } } // Create a matrix of random MVN vectors with distirbution MVN(mean, Sigma) rmvnormal = (C*z)' + meanmat // Return result return(rmvnormal) } real matrix rmvt_mata(real scalar n, real vector delta, real matrix Sigma, real scalar df) { // Perform checks on input variables if (rows(Sigma) ~= cols(Sigma)) { "{error} Scale matrix Sigma (Sigma) must be square." exit(198) } if (length(delta) ~= cols(Sigma)) { "{error} Vector of non-centrality parameters (delta) must be the same length as the dimension of scale matrix Sigma (Sigma)." exit(198) } if ((n < 1) | (mod(n, 1) ~= 0)) { "{error} Number of random vectors to generate (n) must be a strictly positive integer." exit(198) } if ((df < 1) | (mod(df, 1) ~= 0)) { "{error} Degree of freedom (df) must be a positive integer." exit(198) } // Initialise all required variables C = cholesky(Sigma) k = rows(C) // Initialise a matrix with each row the vector delta deltamat = J(n, k, 0) for (i = 1; i <= n; i++) { deltamat[i, 1::k] = delta } // Create a matrix of random normal variables and a vector of random chi2 // variables z = J(k, n, 0) rchi2 = J(n, 1, 0) for (i = 1; i <= n; i++) { for (j = 1; j <= k; j++) { z[j, i] = rnormal(1, 1, 0, 1) } rchi2[i] = rchi2(1, 1, df) } // Create a matrix of random multivariate normal vectors with distribution // MVN(0, Sigma) rmvnormal = (C*z)' // Create a matrix of the desired multivariate t vectors with distribution // MVT(delta, Sigma) rmvt = rmvnormal:/sqrt(rchi2/df) + deltamat // Return result return(rmvt) } real vector mvnormal_mata(real vector lower, real vector upper, real vector mean, real matrix Sigma,| real scalar M, real scalar N, real scalar alpha) { // Check if optional arguments have been specified if (args() < 7) { M = 12 N = 1000 alpha = 3 } // Perform checks on input variables if (cols(Sigma) ~= rows(Sigma)) { "{error} Covariance matrix Sigma (Sigma) must be square." exit(198) } if (length(lower) ~= length(upper)) { "{error} Vector of lower limits (lower) and vector of upper limits (upper) must be of equal length." exit(198) } if (length(lower) ~= length(mean)) { "{error} Mean vector (mean) and vector of lower limits (lower) must be of equal length." exit(198) } if (length(mean) ~= cols(Sigma)) { "{error} Mean vector (mean) must be the same length as the dimension of covariance matrix Sigma (Sigma)." exit(198) } if (length(mean) > 100) { "{error} Only multivariate normal distributions of dimension up to 100 are supported." exit(198) } for (i = 1; i <= length(lower); i++) { if ((lower[i] ~= .) & (upper[i] ~= .)) { if (lower[i] > upper[i]) { "{error} Each lower integration limit (in lower) must be strictly less than the corresponding upper limit (in upper)." exit(198) } } } if ((M < 1) | (mod(M, 1) ~= 0)) { "{error} Number of shifts of the Quasi-Monte-Carlo integration algorithm to use (M) must be a strictly positive integer." exit(198) } if ((N < 1) | (mod(N, 1) ~= 0)) { "{error} Number of samples to use in each shift of the Quasi-Monte-Carlo integration algorithm (N) must be a strictly positive integer." exit(198) } if (alpha <= 0) { "{error} Chosen Monte-Carlo confidence factor (alpha) must be strictly positive." exit(198) } // Initialise all required variables k = rows(Sigma) // If we're dealing with the univariate or bivariate normal case use the // standard functions if (k == 1) { if (lower == .) { lower = -8e+307 } if (upper == .) { upper = 8e+307 } I = normal((upper - mean)/sqrt(Sigma)) - normal((lower - mean)/sqrt(Sigma)) E = 0 } else if (k == 2) { for (i = 1; i <= 2; i++) { if (lower[i] == .) { lower[i] = -8e+307 } if (upper[i] == .) { upper[i] = 8e+307 } } I = binormal((upper[1] - mean[1])/sqrt(Sigma[1, 1]), (upper[2] - mean[2])/sqrt(Sigma[2, 2]), Sigma[1,2]/sqrt(Sigma[1,1]*Sigma[2,2])) + binormal((lower[1] - mean[1])/sqrt(Sigma[1, 1]), (lower[2] - mean[2])/sqrt(Sigma[2, 2]), Sigma[1,2]/sqrt(Sigma[1,1]*Sigma[2,2])) - binormal((upper[1] - mean[1])/sqrt(Sigma[1, 1]), (lower[2] - mean[2])/sqrt(Sigma[2, 2]), Sigma[1,2]/sqrt(Sigma[1,1]*Sigma[2,2])) - binormal((lower[1] - mean[1])/sqrt(Sigma[1, 1]), (upper[2] - mean[2])/sqrt(Sigma[2, 2]), Sigma[1,2]/sqrt(Sigma[1,1]*Sigma[2,2])) E = 0 } // If we're dealing with an actual multivariate case then use the // Quasi-Monte-Carlo Randomised-Lattice Separation-Of-Variables method else if (k > 2) { // Algorithm is for the case with 0 means, so adjust lower and upper // and then re-order variables for maximum efficiency a = lower - mean b = upper - mean C = J(k, k, 0) y = J(1, k - 1, 0) atilde = J(1, k - 1, 0) btilde = J(1, k - 1, 0) atilde[1] = a[1] btilde[1] = b[1] // Loop over each column of Sigma for (i = 1; i <= k - 1; i++) { // Determine variate with minimum expectation args = J(1, k - i + 1, 0) for (j = 1; j <= k - i + 1; j++){ s = j + i - 1 if (i > 1) { if ((a[s] ~= .) & (b[s] ~= .)) { args[j] = normal((b[s] - sum(C[s, 1::(i - 1)]:*y[1::(i - 1)]))/ sqrt(Sigma[s, s] - sum(C[s, 1::(i - 1)]:^2))) - normal((a[s] - sum(C[s, 1::(i - 1)]:*y[1::(i - 1)]))/ sqrt(Sigma[s, s] - sum(C[s, 1::(i - 1)]:^2))) } else if ((a[s] == .) & (b[s] ~= .)) { args[j] = normal((b[s] - sum(C[s, 1::(i - 1)]:*y[1::(i - 1)]))/ sqrt(Sigma[s, s] - sum(C[s, 1::(i - 1)]:^2))) } else if ((b[s] == .) & (a[s] ~= .)) { args[j] = 1 - normal((a[s] - sum(C[s, 1::(i - 1)]:*y[1::(i - 1)]))/ sqrt(Sigma[s, s] - sum(C[s, 1::(i - 1)]:^2))) } else if ((a[s] == .) & (b[s] == .)) { args[j] = 1 } } else { if ((a[s] ~= .) & (b[s] ~= .)) { args[j] = normal(b[s]/sqrt(Sigma[1, 1])) - normal(a[s]/sqrt(Sigma[1, 1])) } else if ((a[s] == .) & (b[s] ~= .)) { args[j] = normal(b[s]/sqrt(Sigma[1, 1])) } else if ((b[s] == .) & (a[s] ~= .)) { args[j] = 1 - normal(a[s]/sqrt(Sigma[1, 1])) } else if ((a[s] == .) & (b[s] == .)) { args[j] = 1 } } } minindex(args, 1, ii, ww) m = i - 1 + ii[1] // Change elements m and i of a, b, Sigma and C 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){ 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] // Compute next column of C and next value of atilda and btilda 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] } atilde[i] = (a[i] - sum(C[i, 1::(i - 1)]:*y[1::(i - 1)]))/C[i, i] btilde[i] = (b[i] - sum(C[i, 1::(i - 1)]:*y[1::(i - 1)]))/C[i, i] } else { C[i, i] = sqrt(Sigma[i, i]) C[2::k, i] = Sigma[2::k, i]/C[i, i] } // Compute next value of y using normalden y[i] = (normalden(atilde[i]) - normalden(btilde[i]))/ (normal(btilde[i]) - normal(atilde[i])) } // Set the final element of C C[k, k] = sqrt(Sigma[k, k] - sum(C[k, 1::(k - 1)]:^2)) // List of the first 100 primes to use in the integration algorithm p = (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) // Initialise all required variables I = 0 V = 0 sqp = p[1::(k - 1)]:^0.5 d = J(1, k, 0) e = J(1, k, 1) f = J(1, k, 0) // First elements of d, e and f are always the same if (a[1] ~= .) { d[1] = normal(a[1]/C[1, 1]) } if (b[1] ~= .) { e[1] = normal(b[1]/C[1, 1]) } f[1] = e[1] - d[1] // Perform M shifts of the Quasi-Monte-Carlo integration algorithm for (i = 1; i <= M; i++) { Ii = 0 // We require k - 1 random uniform numbers Delta = runiform(1, k - 1) // Use N samples in each shift for (j = 1; j <= N; j++) { // Loop to compute other values of d, e and f for (l = 2; l <= k; l++) { y[l - 1] = invnormal(d[l - 1] + abs(2*(mod(j*sqp[l - 1] + Delta[l - 1], 1)) - 1)* (e[l - 1] - d[l - 1])) if (a[l] ~= .) { d[l] = normal((a[l] - sum(C[l, 1::(l - 1)]:*y[1::(l - 1)]))/C[l, l]) } if (b[l] ~= .) { e[l] = normal((b[l] - sum(C[l, 1::(l - 1)]:*y[1::(l - 1)]))/C[l, l]) } f[l] = (e[l] - d[l])*f[l - 1] } Ii = Ii + (f[k] - Ii)/j } // Update the values of the variables delta = (Ii - I)/i I = I + delta V = (i - 2)*V/i + delta^2 E = alpha*sqrt(V) } } // Return result return((I, E)) } real vector mvt_mata(real vector lower, real vector upper, real vector delta, real matrix Sigma, real scalar df,| real scalar M, real scalar N, real scalar alpha) { // Check if optional arguments have been specified if (args() < 8) { M = 12 N = 1000 alpha = 3 } // Perform checks on input variables if (cols(Sigma) ~= rows(Sigma)) { "{error} Scale matrix Sigma (Sigma) must be square." exit(198) } if (length(lower) ~= length(upper)) { "{error} Vector of lower limits (lower) and vector of upper limits (upper) must be of equal length." exit(198) } if (length(lower) ~= length(delta)) { "{error} Vector of lower limits (lower) and vector of non-centrality parameters (delta) must be of equal length." exit(198) } if (length(delta) ~= cols(Sigma)) { "{error} Vector of non-centrality parameters (delta) must be the same length as the dimension of scale matrix Sigma (sigma)." exit(198) } if (length(delta) > 100) { "{error} Only multivariate t distributions of dimension up to 100 are supported." exit(198) } for (i = 1; i <= length(lower); i++) { if ((lower[i] ~= .) & (upper[i] ~= .)) { if (lower[i] > upper[i]) { "{error} Each lower integration limit (in lower) must be strictly less than the corresponding upper limit (in upper)." exit(198) } } } if ((df < 1) | (mod(df, 1) ~= 0)) { "{error} Degree of freedom (df) must be a positive integer." exit(198) } if ((M < 1) | (mod(M, 1) ~= 0)) { "{error} Number of shifts of the Quasi-Monte-Carlo integration algorithm to use (M) must be a strictly positive integer." exit(198) } if ((N < 1) | (mod(N, 1) ~= 0)) { "{error} Number of samples to use in each shift of the Quasi-Monte-Carlo integration algorithm (N) must be a strictly positive integer." exit(198) } if (alpha <= 0) { "{error} Chosen Monte-Carlo confidence factor (alpha) must be strictly positive." exit(198) } // Initialise all required variables k = rows(Sigma) // If we're dealing with a univariate t variable use standard functions if (k == 1) { I = (1 - ttail(df, upper - delta)) - (1 - ttail(df, lower - delta)) E = 0 } // If we're dealing with an actual multivariate case then use the // Quasi-Monte-Carlo Randomised-Lattice Separation-Of-Variables method else { // Algorithm is for the case with 0 non-centrality parameters, so adjust // lower and upper and then re-order variables for maximum efficiency a = lower - delta b = upper - delta C = J(k, k, 0) u = J(1, k - 1, 0) y = J(1, k - 1, 0) atilde = J(1, k - 1, 0) btilde = J(1, k - 1, 0) atilde[1] = a[1] btilde[1] = b[1] // Loop over each column of Sigma for (i = 1; i <= k - 1; i++) { // Determine variate with minimum expectation args = J(1, k - i + 1, 0) for (j = 1; j <= k - i + 1; j++){ s = j + i - 1 if (i > 1) { if ((a[s] ~= .) & (b[s] ~= .)) { args[j] = (1 - ttail(df + i - 1, sqrt((df + i - 1)/ (df + sum(y[1::(i - 1)]:^2)))* ((b[s] - sum(C[s, 1::(i - 1)]:*y[1::(i - 1)]))/ sqrt(Sigma[s, s] - sum(C[s, 1::(i - 1)]:^2))))) - (1 - ttail(df + i - 1, sqrt((df + i - 1)/ (df + sum(y[1::(i - 1)]:^2)))* ((a[s] - sum(C[s, 1::(i - 1)]:*y[1::(i - 1)]))/ sqrt(Sigma[s, s] - sum(C[s, 1::(i - 1)]:^2))))) } else if ((a[s] == .) & (b[s] ~= .)) { args[j] = 1 - ttail(df + i - 1, sqrt((df + i - 1)/ (df + sum(y[1::(i - 1)]:^2)))* ((b[s] - sum(C[s, 1::(i - 1)]:*y[1::(i - 1)]))/ sqrt(Sigma[s, s] - sum(C[s, 1::(i - 1)]:^2)))) } else if ((b[s] == .) & (a[s] ~= .)) { args[j] = 1 - (1 - ttail(df + i - 1, sqrt((df + i - 1)/ (df + sum(y[1::(i - 1)]:^2)))* ((a[s] - sum(C[s, 1::(i - 1)]:*y[1::(i - 1)]))/ sqrt(Sigma[s, s] - sum(C[s, 1::(i - 1)]:^2))))) } else if ((a[s] == .) & (b[s] == .)) { args[j] = 1 } } else { if ((a[s] ~= .) & (b[s] ~= .)) { args[j] = (1 - ttail(df, b[s]/sqrt(diag(Sigma[1, 1])))) - (1 - ttail(df, a[s]/sqrt(diag(Sigma[1, 1])))) } else if ((a[s] == .) & (b[s] ~= .)) { args[j] = 1 - ttail(df, b[s]/sqrt(diag(Sigma[1, 1]))) } else if ((b[s] == .) & (a[s] ~= .)) { args[j] = 1 - (1 - ttail(df, a[s]/sqrt(diag(Sigma[1, 1])))) } else if ((a[s] == .) & (b[s] == .)) { args[j] = 1 } } } minindex(args, 1, ii, ww) m = i - 1 + ii[1] // Change elements m and i of a, b, Sigma and C 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){ 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] // Compute next column of C and next value of atilde and btilde 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] } atilde[i] = (a[i] - sum(C[i, 1::(i - 1)]:*y[1::(i - 1)]))/C[i, i] btilde[i] = (b[i] - sum(C[i, 1::(i - 1)]:*y[1::(i - 1)]))/C[i, i] } else { C[i, i] = sqrt(Sigma[i, i]) C[2::k, i] = Sigma[2::k, i]/C[i, i] } // Compute next value of u using integrate arg = (df, i) u[i] = (gamma((df + 1)/2)/(gamma((df + i - 1)/2)*((df + i - 1)*pi())^0.5))* integrate_MVTNORM(atilde[i], btilde[i], arg)/ ((1 - ttail(df + i - 1, btilde[i])) - (1 - ttail(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)) } } // Set the final element of C C[k, k] = sqrt(Sigma[k, k] - sum(C[k, 1::(k - 1)]:^2)) // List of the first 100 primes to use in randomised lattice algorithm p = (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) // Initialise all required variables I = 0 V = 0 sqp = p[1::(k - 1)]:^0.5 d = J(1, k, 0) e = J(1, k, 1) f = J(1, k, 0) // First elements of d, e and f are always the same if (a[1] ~= .){ d[1] = 1 - ttail(df, a[1]/C[1, 1]) } if (b[1] ~= .){ e[1] = 1 - ttail(df, b[1]/C[1, 1]) } f[1] = e[1] - d[1] // Perform M shifts of the Quasi-Monte-Carlo integration algorithm for (i = 1; i <= M; i++) { Ii = 0 // We require k - 1 random uniform numbers Delta = runiform(1, k - 1) // Use N samples in each shift for (j = 1; j <= N; j++) { u[1] = invttail(df, 1 - (d[1] + abs(2*(mod(j*sqp[1] + Delta[1], 1)) - 1)*(e[1] - d[1]))) y[1] = u[1] // Loop to compute other values of d, e and f for (l = 2; l <= k; l++) { if (a[l] ~= .) { d[l] = 1 - ttail(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])) } if (b[l] ~= .){ e[l] = 1 - ttail(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] = invttail(df + l - 1, 1 - (d[l] + abs(2*(mod(j*sqp[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 } // Update the values of the variables delta = (Ii - I)/i I = I + delta V = (i - 2)*V/i + delta^2 E = alpha*sqrt(V) } } // Return result return((I, E)) } real vector invmvnormal_mata(real scalar p, real vector mean, real matrix Sigma, string tail,| real scalar M, real scalar N, real scalar alpha, real scalar itermax, real scalar tol) { // Check if optional arguments have been specified if (args() < 9) { M = 12 N = 1000 alpha = 3 itermax = 1000000 tol = 0.000001 } // Perform checks on input variables if ((p <= 0) | (p >= 1)) { "{error} Probability (p) must be between 0 and 1." exit(198) } if (cols(Sigma) ~= rows(Sigma)) { "{error} Covariance matrix Sigma (Sigma) must be square." exit(198) } if (length(mean) ~= cols(Sigma)) { "{error} Mean vector (mean) must be the same length as the dimension of covariance matrix Sigma (Sigma)." exit(198) } if (length(mean) > 100) { "{error} Only multivariate normal distributions of dimension up to 100 are supported." exit(198) } if ((tail ~= "lower") & (tail ~= "upper") & (tail ~= "both")) { "{error} tail must be set to one of lower, upper or both." exit(198) } if ((M < 1) | (mod(M, 1) ~= 0)) { "{error} Number of shifts of the Quasi-Monte-Carlo integration algorithm to use (M) must be a strictly positive integer." exit(198) } if ((N < 1) | (mod(N, 1) ~= 0)) { "{error} Number of samples to use in each shift of the Quasi-Monte-Carlo integration algorithm (N) must be a strictly positive integer." exit(198) } if (alpha <= 0) { "{error} Chosen Monte-Carlo confidence factor (alpha) must be strictly positive." exit(198) } if ((itermax < 1) | (mod(itermax, 1) ~= 0)) { "{error} Number of allowed iterations in the interval bisection quantile finding algorithm (itermax) must be a strictly positive integer." exit(198) } if (tol < 0) { "{error} The tolerance in the interval bisection quantile finding algorithm (tol) must be strictly positive." exit(198) } // Determine if trivial case is desired if ((p == 0) | (p == 1)) { q = . error = 0 flag = 0 fq = 0 iter = 0 } else { // Initialise all required variables C = cholesky(Sigma) k = rows(Sigma) // If we're dealing with the univariate normal case use the standard // function if (k == 1) { if (tail == "upper") { p = 1 - p } else if (tail == "both") { p = 0.5 + 0.5*p } q = invnormal(p) + mean error = 0 flag = 0 fq = 0 iter = 0 } // If we're dealing with an actual multivariate case then use the // Quasi-Monte-Carlo Randomised-Lattice Separation-Of-Variables method, // along with modified interval bisection else { // List of the first 100 primes to use in randomised lattice algorithm primes = (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) sqp = primes[1::(k - 1)]:^0.5 // Determine sensible starting interval to bisect over b = 10 if (tail == "both") { a = 10^-6 } else { a = -10 } fa = qFinder_mvnormal(a, p, k, tail, M, N, alpha, mean, Sigma, sqp) fb = qFinder_mvnormal(b, p, k, tail, M, N, alpha, mean, Sigma, sqp) if (fa == 0) { q = a error = 0 flag = 0 fq = 0 iter = 0 } else if (fb == 0) { q = b error = 0 flag = 0 fq = 0 iter = 0 } else { if (tail == "both") { while ((fa > 0) & (a >= 10^-20)) { a = a/2 fa = qFinder_mvnormal(a, p, k, tail, M, N, alpha, mean, Sigma, sqp) } while ((fb < 0) & (b <= 10^6)) { b = 2*b fb = qFinder_mvnormal(b, p, k, tail, M, N, alpha, mean, Sigma, sqp) } } else { while (((fa > 0 & fb > 0) | (fa < 0 & fb < 0)) & (a >= -10^6)) { a = 2*a b = 2*b fa = qFinder_mvnormal(a, p, k, tail, M, N, alpha, mean, Sigma, sqp) fb = qFinder_mvnormal(b, p, k, tail, M, N, alpha, mean, Sigma, sqp) } } if (fa == 0) { q = a error = 0 flag = 0 fq = 0 iter = 0 } else if (fb == 0) { q = b error = 0 flag = 0 fq = 0 iter = 0 } // Perform modified interval bisection else if ((fa < 0 & fb > 0) | (fa > 0 & fb < 0)) { iter = 1 while (iter <= itermax) { c = a - ((b - a)/(fb - fa))*fa fc = qFinder_mvnormal(c, p, k, tail, M, N, alpha, mean, Sigma, sqp) if ((fc == 0) | ((b - a)/2 < tol)) { break } else { if (((fa < 0) & (fc < 0)) | ((fa > 0) & (fc > 0))) { a = c fa = fc } else { b = c fb = fc } iter = iter + 1 } } q = c error = (b - a)/2 fq = fc if (iter == itermax + 1) { flag = 1 } else { flag = 0 } } else { q = . error = . flag = 2 fq = . iter = . } } } } // Return result return((q, error, flag, fq, iter)) } function qFinder_mvnormal(q, p, k, tail, M, N, alpha, mean, Sigma, sqp) { // Initialise all required variables 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 (k == 2) { for (i = 1; i <= 2; i++) { if (a[i] == .) { a[i] = -8e+307 } if (b[i] == .) { b[i] = 8e+307 } } I = binormal((b[1] - mean[1])/sqrt(Sigma[1, 1]), (b[2] - mean[2])/sqrt(Sigma[2, 2]), Sigma[1,2]/sqrt(Sigma[1,1]*Sigma[2,2])) + binormal((a[1] - mean[1])/sqrt(Sigma[1, 1]), (a[2] - mean[2])/sqrt(Sigma[2, 2]), Sigma[1,2]/sqrt(Sigma[1,1]*Sigma[2,2])) - binormal((b[1] - mean[1])/sqrt(Sigma[1, 1]), (a[2] - mean[2])/sqrt(Sigma[2, 2]), Sigma[1,2]/sqrt(Sigma[1,1]*Sigma[2,2])) - binormal((a[1] - mean[1])/sqrt(Sigma[1, 1]), (b[2] - mean[2])/sqrt(Sigma[2, 2]), Sigma[1,2]/sqrt(Sigma[1,1]*Sigma[2,2])) } else { // Algorithm is for the case with 0 means, so adjust lower and upper // and then re-order variables for maximum efficiency a = a - mean b = b - mean C = J(k, k, 0) y = J(1, k - 1, 0) atilde = J(1, k - 1, 0) btilde = J(1, k - 1, 0) atilde[1] = a[1] btilde[1] = b[1] // Loop over each column of Sigma for (i = 1; i <= k - 1; i++) { // Determine variate with minimum expectation args = J(1, k - i + 1, 0) for (j = 1; j <= k - i + 1; j++){ s = j + i - 1 if (i > 1) { if ((a[s] ~= .) & (b[s] ~= .)) { args[j] = normal((b[s] - sum(C[s, 1::(i - 1)]:*y[1::(i - 1)]))/ sqrt(Sigma[s, s] - sum(C[s, 1::(i - 1)]:^2))) - normal((a[s] - sum(C[s, 1::(i - 1)]:*y[1::(i - 1)]))/ sqrt(Sigma[s, s] - sum(C[s, 1::(i - 1)]:^2))) } else if ((a[s] == .) & (b[s] ~= .)) { args[j] = normal((b[s] - sum(C[s, 1::(i - 1)]:*y[1::(i - 1)]))/ sqrt(Sigma[s, s] - sum(C[s, 1::(i - 1)]:^2))) } else if ((b[s] == .) & (a[s] ~= .)) { args[j] = 1 - normal((a[s] - sum(C[s, 1::(i - 1)]:*y[1::(i - 1)]))/ sqrt(Sigma[s, s] - sum(C[s, 1::(i - 1)]:^2))) } else if ((a[s] == .) & (b[s] == .)) { args[j] = 1 } } else { if ((a[s] ~= .) & (b[s] ~= .)) { args[j] = normal(b[s]/sqrt(Sigma[1, 1])) - normal(a[s]/sqrt(Sigma[1, 1])) } else if ((a[s] == .) & (b[s] ~= .)) { args[j] = normal(b[s]/sqrt(Sigma[1, 1])) } else if ((b[s] == .) & (a[s] ~= .)) { args[j] = 1 - normal(a[s]/sqrt(Sigma[1, 1])) } else if ((a[s] == .) & (b[s] == .)) { args[j] = 1 } } } minindex(args, 1, ii, ww) m = i - 1 + ii[1] // Change elements m and i of a, b, Sigma and C 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){ 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] // Compute next column of C and next value of atilda and btilda 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] } atilde[i] = (a[i] - sum(C[i, 1::(i - 1)]:*y[1::(i - 1)]))/C[i, i] btilde[i] = (b[i] - sum(C[i, 1::(i - 1)]:*y[1::(i - 1)]))/C[i, i] } else { C[i, i] = sqrt(Sigma[i, i]) C[2::k, i] = Sigma[2::k, i]/C[i, i] } // Compute next value of y using normalden y[i] = (normalden(atilde[i]) - normalden(btilde[i]))/ (normal(btilde[i]) - normal(atilde[i])) } // Set the final element of C C[k, k] = sqrt(Sigma[k, k] - sum(C[k, 1::(k - 1)]:^2)) I = 0 d = J(1, k, 0) e = J(1, k, 1) f = J(1, k, 0) w = J(1, k - 1, 0) // First elements of d, e and f are always the same if (a[1] ~= .) { d[1] = normal(a[1]/C[1, 1]) } if (b[1] ~= .) { e[1] = normal(b[1]/C[1, 1]) } f[1] = e[1] - d[1] // Perform M shifts of the Quasi-Monte-Carlo integration algorithm for (i = 1; i <= M; i++) { Ii = 0 // We require k - 1 random uniform numbers Delta = runiform(1, k - 1) // Use N samples in each shift for (j = 1; j <= N; j++) { // Loop to compute other values of d, e and f for (l = 2; l <= k; l++) { y[l - 1] = invnormal(d[l - 1] + abs(2*(mod(j*sqp[l - 1] + Delta[l - 1], 1)) - 1)* (e[l - 1] - d[l - 1])) if (a[l] ~= .) { d[l] = normal((a[l] - sum(C[l, 1::(l - 1)]:*y[1::(l - 1)]))/C[l, l]) } if (b[l] ~= .) { e[l] = normal((b[l] - sum(C[l, 1::(l - 1)]:*y[1::(l - 1)]))/C[l, l]) } f[l] = (e[l] - d[l])*f[l - 1] } Ii = Ii + (f[k] - Ii)/j } // Update the values of the variables I = I + (Ii - I)/i } } return(I - p) } real vector invmvt_mata(real scalar p, real vector delta, real matrix Sigma, real scalar df, string tail,| real scalar M, real scalar N, real scalar alpha, real scalar itermax, real scalar tol) { // Check if optional arguments have been specified if (args() < 9) { M = 12 N = 1000 alpha = 3 itermax = 1000000 tol = 0.000001 } // Perform checks on input variables if ((p <= 0) | (p >= 1)) { "{error} Probability (p) must be between 0 and 1." exit(198) } if (cols(Sigma) ~= rows(Sigma)) { "{error} Scale matrix Sigma (Sigma) must be square." exit(198) } if (length(delta) ~= cols(Sigma)) { "{error} Vector of non-centrality parameters (delta) must be the same length as the dimension of scale matrix Sigma (sigma)." exit(198) } if (length(delta) > 100) { "{error} Only multivariate t distributions of dimension up to 100 are supported." exit(198) } if ((df < 1) | (mod(df, 1) ~= 0)) { "{error} Degree of freedom (df) must be a positive integer." exit(198) } if ((tail ~= "lower") & (tail ~= "upper") & (tail ~= "both")) { "{error} tail must be set to one of lower, upper or both." exit(198) } if ((M < 1) | (mod(M, 1) ~= 0)) { "{error} Number of shifts of the Quasi-Monte-Carlo integration algorithm to use (M) must be a strictly positive integer." exit(198) } if ((N < 1) | (mod(N, 1) ~= 0)) { "{error} Number of samples to use in each shift of the Quasi-Monte-Carlo integration algorithm (N) must be a strictly positive integer." exit(198) } if (alpha <= 0) { "{error} Chosen Monte-Carlo confidence factor (alpha) must be strictly positive." exit(198) } if ((itermax < 1) | (mod(itermax, 1) ~= 0)) { "{error} Number of allowed iterations in the interval bisection quantile finding algorithm (itermax) must be a strictly positive integer." exit(198) } if (tol < 0) { "{error} The tolerance in the interval bisection quantile finding algorithm (tol) must be strictly positive." exit(198) } // Determine if trivial case is desired if ((p == 0) | (p == 1)) { q = . error = 0 flag = 0 fq = 0 iter = 0 } else { // Initialise all required variables k = rows(Sigma) // If we're dealing with the univariate normal case use the standard // function if (k == 1) { if (tail == "upper") { p = 1 - p } else if (tail == "both") { p = 0.5 + 0.5*p } q = invttail(df, 1 - p) + delta error = 0 flag = 0 fq = 0 iter = 0 } // If we're dealing with an actual multivariate case then use the // Quasi-Monte-Carlo Randomised-Lattice Separation-Of-Variables method, // along with optimize else { // List of the first 100 primes to use in randomised lattice algorithm primes = (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) sqp = primes[1::(k - 1)]:^0.5 // Determine sensible starting interval to bisect over b = 20 if (tail == "both") { a = 10^-6 } else { a = -20 } fa = qFinder_mvt(a, p, k, df, tail, M, N, alpha, delta, Sigma, sqp) fb = qFinder_mvt(b, p, k, df, tail, M, N, alpha, delta, Sigma, sqp) if (fa == 0) { q = a error = 0 flag = 0 fq = 0 iter = 0 } else if (fb == 0) { q = b error = 0 flag = 0 fq = 0 iter = 0 } else { if (tail == "both") { while ((fa > 0) & (a >= 10^-20)) { a = a/2 fa = qFinder_mvt(a, p, k, df, tail, M, N, alpha, delta, Sigma, sqp) } while ((fb < 0) & (b <= 10^6)) { b = 2*b fb = qFinder_mvt(b, p, k, df, tail, M, N, alpha, delta, Sigma, sqp) } } else { while (((fa > 0 & fb > 0) | (fa < 0 & fb < 0)) & (a >= -10^6)) { a = 2*a b = 2*b fa = qFinder_mvt(a, p, k, df, tail, M, N, alpha, delta, Sigma, sqp) fb = qFinder_mvt(b, p, k, df, tail, M, N, alpha, delta, Sigma, sqp) } } if (fa == 0) { q = a error = 0 flag = 0 fq = 0 iter = 0 } else if (fb == 0) { q = b error = 0 flag = 0 fq = 0 iter = 0 } // Perform modified interval bisection else if ((fa < 0 & fb > 0) | (fa > 0 & fb < 0)) { iter = 1 while (iter <= itermax) { c = a - ((b - a)/(fb - fa))*fa fc = qFinder_mvt(c, p, k, df, tail, M, N, alpha, delta, Sigma, sqp) if ((fc == 0) | ((b - a)/2 < tol)) { break } else { if (((fa < 0) & (fc < 0)) | ((fa > 0) & (fc > 0))) { a = c fa = fc } else { b = c fb = fc } iter = iter + 1 } } q = c error = (b - a)/2 fq = fc if (iter == itermax + 1) { flag = 1 } else { flag = 0 } } else { q = . error = . flag = 2 fq = . iter = . } } } } // Return result return((q, error, flag, fq, iter)) } function qFinder_mvt(q, p, k, df, tail, M, N, alpha, delta, Sigma, sqp) { // Algorithm is for the case with 0 non-centrality parameters, so adjust // lower and upper and then re-order variables for maximum efficiency if (tail == "lower") { a = J(1, k, .) b = J(1, k, q) - delta } else if (tail == "upper") { a = J(1, k, q) - delta b = J(1, k, .) } else { a = J(1, k, -q) - delta b = J(1, k, q) - delta } C = J(k, k, 0) u = J(1, k - 1, 0) y = J(1, k - 1, 0) atilde = J(1, k - 1, 0) btilde = J(1, k - 1, 0) atilde[1] = a[1] btilde[1] = b[1] // Loop over each column of Sigma for (i = 1; i <= k - 1; i++) { // Determine variate with minimum expectation args = J(1, k - i + 1, 0) for (j = 1; j <= k - i + 1; j++){ s = j + i - 1 if (i > 1) { if ((a[s] ~= .) & (b[s] ~= .)) { args[j] = (1 - ttail(df + i - 1, sqrt((df + i - 1)/ (df + sum(y[1::(i - 1)]:^2)))* ((b[s] - sum(C[s, 1::(i - 1)]:*y[1::(i - 1)]))/ sqrt(Sigma[s, s] - sum(C[s, 1::(i - 1)]:^2))))) - (1 - ttail(df + i - 1, sqrt((df + i - 1)/ (df + sum(y[1::(i - 1)]:^2)))* ((a[s] - sum(C[s, 1::(i - 1)]:*y[1::(i - 1)]))/ sqrt(Sigma[s, s] - sum(C[s, 1::(i - 1)]:^2))))) } else if ((a[s] == .) & (b[s] ~= .)) { args[j] = 1 - ttail(df + i - 1, sqrt((df + i - 1)/ (df + sum(y[1::(i - 1)]:^2)))* ((b[s] - sum(C[s, 1::(i - 1)]:*y[1::(i - 1)]))/ sqrt(Sigma[s, s] - sum(C[s, 1::(i - 1)]:^2)))) } else if ((b[s] == .) & (a[s] ~= .)) { args[j] = 1 - (1 - ttail(df + i - 1, sqrt((df + i - 1)/ (df + sum(y[1::(i - 1)]:^2)))* ((a[s] - sum(C[s, 1::(i - 1)]:*y[1::(i - 1)]))/ sqrt(Sigma[s, s] - sum(C[s, 1::(i - 1)]:^2))))) } else if ((a[s] == .) & (b[s] == .)) { args[j] = 1 } } else { if ((a[s] ~= .) & (b[s] ~= .)) { args[j] = (1 - ttail(df, b[s]/sqrt(diag(Sigma[1, 1])))) - (1 - ttail(df, a[s]/sqrt(diag(Sigma[1, 1])))) } else if ((a[s] == .) & (b[s] ~= .)) { args[j] = 1 - ttail(df, b[s]/sqrt(diag(Sigma[1, 1]))) } else if ((b[s] == .) & (a[s] ~= .)) { args[j] = 1 - (1 - ttail(df, a[s]/sqrt(diag(Sigma[1, 1])))) } else if ((a[s] == .) & (b[s] == .)) { args[j] = 1 } } } minindex(args, 1, ii, ww) m = i - 1 + ii[1] // Change elements m and i of a, b, Sigma and C 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){ 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] // Compute next column of C and next value of atilde and btilde 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] } atilde[i] = (a[i] - sum(C[i, 1::(i - 1)]:*y[1::(i - 1)]))/C[i, i] btilde[i] = (b[i] - sum(C[i, 1::(i - 1)]:*y[1::(i - 1)]))/C[i, i] } else { C[i, i] = sqrt(Sigma[i, i]) C[2::k, i] = Sigma[2::k, i]/C[i, i] } // Compute next value of u using integrate arg = (df, i) u[i] = (gamma((df + 1)/2)/(gamma((df + i - 1)/2)*((df + i - 1)*pi())^0.5))* integrate_MVTNORM(atilde[i], btilde[i], arg)/ ((1 - ttail(df + i - 1, btilde[i])) - (1 - ttail(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)) } } // Set the final element of C C[k, k] = sqrt(Sigma[k, k] - sum(C[k, 1::(k - 1)]:^2)) // Initialise all required variables I = 0 d = J(1, k, 0) e = J(1, k, 1) f = J(1, k, 0) // First elements of d, e and f are always the same if (a[1] ~= .){ d[1] = 1 - ttail(df, a[1]/C[1, 1]) } if (b[1] ~= .){ e[1] = 1 - ttail(df, b[1]/C[1, 1]) } f[1] = e[1] - d[1] // Perform M shifts of the Quasi-Monte-Carlo integration algorithm for (i = 1; i <= M; i++) { Ii = 0 // We require k - 1 random uniform numbers Delta = runiform(1, k - 1) // Use N samples in each shift for (j = 1; j <= N; j++) { u[1] = invttail(df, 1 - (d[1] + abs(2*(mod(j*sqp[1] + Delta[1], 1)) - 1)*(e[1] - d[1]))) y[1] = u[1] // Loop to compute other values of d, e and f for (l = 2; l <= k; l++) { if (a[l] ~= .) { d[l] = 1 - ttail(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])) } if (b[l] ~= .){ e[l] = 1 - ttail(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] = invttail(df + l - 1, 1 - (d[l] + abs(2*(mod(j*sqp[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 } // Update the values of the variables I = I + (Ii - I)/i } return(I - p) } real scalar integrate_MVTNORM(real scalar lower, real scalar upper, real vector arg) { if ((lower == . & upper == .) | (lower == 0 & upper == .) |(lower ~= . & upper ~= .)) { return(integrate_main_MVTNORM(lower, upper, arg)) } else if (lower == . & upper ~= .) { return(integrate_main_MVTNORM(0, upper, arg) + integrate_main_MVTNORM(0, ., arg)) } else if (lower ~= 0 & upper == .) { return(integrate_main_MVTNORM(lower, 0, arg) + integrate_main_MVTNORM(0, ., arg)) } else { return(integrate_main_MVTNORM(lower, upper, arg)) } } real matrix integrate_main_MVTNORM(real lower, real upper, real vector arg) { if (lower ~= . & upper ~= .) { rw = (.9997137268, .9984919506, .9962951347, .993124937, .9889843952, .9838775407, .9778093585, .9707857758, .9628136543, .9539007829, .9440558701, .933288535, .9216092981, .909029571, .895561645, .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, .292617188, .2625881204, .2323024818, .2017898641, .1710800805, .1402031372, .1091892036, .0780685828, .0468716824, .0156289844, -.0156289844, -.0468716824, -.0780685828, -.1091892036, -.1402031372, -.1710800805, -.2017898641, -.2323024818, -.2625881204, -.292617188, -.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, -.895561645, -.909029571, -.9216092981, -.933288535, -.9440558701, -.9539007829, -.9628136543, -.9707857758, -.9778093585, -.9838775407, -.9889843952, -.993124937, -.9962951347, -.9984919506, -.9997137268 \ .0007346345, .0017093927, .0026839254, .0036559612, .0046244501, .005588428, .0065469485, .0074990733, .0084438715, .0093804197, .0103078026, .011225114, .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, .030798379, .0309504789, .0310723374, .0311638357, .0312248843, .0312554235, .0312554235, .0312248843, .0311638357, .0310723374, .0309504789, .030798379, .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, .011225114, .0103078026, .0093804197, .0084438715, .0074990733, .0065469485, .005588428, .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.9841128, 355.2613119, 339.4351019, 325.6912634, 313.329534, 301.9858553, 291.4401336, 281.5463283, 272.20117, 263.3281685, 254.8686293, 246.776241, 239.0136298, 231.550068, 224.3598948, 217.4213933, 210.7159729, 204.2275596, 197.9421331, 191.8473694, 185.9323602, 180.1873909, 174.6037612, 169.1736398, 163.8899456, 158.7462485, 153.7366875, 148.8559014, 144.09897, 139.4613646, 134.938905, 130.5277232, 126.2242308, 122.0250921, 117.9271991, 113.9276512, 110.0237356, 106.2129115, 102.4927949, 98.86114605, 95.31585735, 91.85494326, 88.47653082, 85.17885121, 81.96023222, 78.81909132, 75.75392947, 72.76332543, 69.84593064, 67.00046452, 64.22571012, 61.52051029, 58.88376398, 56.31442299, 53.81148892, 51.37401033, 49.00108021, 46.69183354, 44.44544511, 42.26112748, 40.13812906, 38.07573237, 36.07325241, 34.13003512, 32.245456, 30.41891877, 28.64985415, 26.93771873, 25.28199387, 23.68218476, 22.13781945, 20.64844797, 19.21364158, 17.83299195, 16.50611047, 15.2326276, 14.01219225, 12.84447118, 11.72914849, 10.66592509, 9.654518244, 8.694661114, 7.786102378, 6.928605829, 6.121950031, 5.365927986, 4.660346836, 4.005027582, 3.399804827, 2.844526543, 2.33905385, 1.883260826, 1.47703433, 1.120273835, .8128912841, .5548109376, .345969181, .1863141021, .075803612, .014386147 \ 7.59641e-96, 3.14629e-94, 9.41175e-96, 7.0997e-106, 3.8299e-102, 2.95362e-97, 2.13367e-97, 2.22822e-97, 1.32643e-95, 5.93228e-99, 1.29153e-99, 1.1580e-100, 1.65794e-95, 4.5382e-100, 1.81558e-90, 1.43778e-93, 2.02036e-91, 1.27976e-88, 6.70480e-86, 2.88488e-83, 1.03790e-80, 3.15247e-78, 8.15375e-76, 1.80986e-73, 3.47185e-71, 5.79263e-69, 8.45496e-67, 1.08537e-64, 1.23138e-62, 1.24024e-60, 1.11357e-58, 8.94734e-57, 6.45627e-55, 4.19775e-53, 2.46681e-51, 1.31398e-49, 6.36127e-48, 2.80604e-46, 1.13048e-44, 4.16886e-43, 1.41014e-41, 4.38379e-40, 1.25483e-38, 3.31306e-37, 8.08163e-36, 1.82420e-34, 3.81586e-33, 7.40743e-32, 1.33621e-30, 2.24265e-29, 3.50627e-28, 5.11237e-27, 6.95920e-26, 8.85323e-25, 1.05359e-23, 1.17402e-22, 1.22601e-21, 1.20085e-20, 1.10409e-19, 9.53625e-19, 7.74309e-18, 5.91443e-17, 4.25261e-16, 2.88012e-15, 1.83835e-14, 1.10650e-13, 6.28352e-13, 3.36821e-12, 1.70506e-11, 8.15480e-11, 3.68633e-10, 1.57560e-09, 6.36971e-09, 2.43643e-08, 8.82006e-08, 3.02264e-07, 9.80834e-07, 3.01427e-06, 8.77431e-06, .0000241958, .0000632109, .0001564521, .0003668548, .0008148716, .0017143197, .00341498, .006438951, .0114854424, .0193678281, .0308463086, .0463401336, .0655510093, .0870966385, .1083141121, .1254070908, .1340433397, .130356613, .1121151033, .0796767462, .0363926059) sum = rw[2, ]:*exp(rw[1, ]):*sf(rw[1, ], arg) return(quadrowsum(sum)) } else if (lower == . & upper == .) { rw = (13.40648734, 12.82379975, 12.34296422, 11.91506194, 11.5214154, 11.15240439, 10.80226075, 10.46718542, 10.14450994, 9.832269808, 9.528965823, 9.23342089, 8.944689217, 8.661996168, 8.38469694, 8.112247311, 7.844182384, 7.580100808, 7.319652822, 7.06253106, 6.808463353, 6.557207032, 6.308544361, 6.062278833, 5.818232135, 5.576241649, 5.33615836, 5.097845105, 4.861175092, 4.626030636, 4.392302079, 4.159886855, 3.928688683, 3.698616859, 3.469585636, 3.24151368, 3.01432358, 2.787941424, 2.562296402, 2.337320464, 2.112947996, 1.889115537, 1.665761509, 1.44282597, 1.220250391, .9979774361, .7759507615, .5541148236, .3324146923, .1107958724, -.1107958724, -.3324146923, -.5541148236, -.7759507615, -.9979774361, -1.220250391, -1.44282597, -1.665761509, -1.889115537, -2.112947996, -2.337320464, -2.562296402, -2.787941424, -3.01432358, -3.24151368, -3.469585636, -3.698616859, -3.928688683, -4.159886855, -4.392302079, -4.626030636, -4.861175092, -5.097845105, -5.33615836, -5.576241649, -5.818232135, -6.062278833, -6.308544361, -6.557207032, -6.808463353, -7.06253106, -7.319652822, -7.580100808, -7.844182384, -8.112247311, -8.38469694, -8.661996168, -8.944689217, -9.23342089, -9.528965823, -9.832269808, -10.14450994, -10.46718542, -10.80226075, -11.15240439, -11.5214154, -11.91506194, -12.34296422, -12.82379975, -13.40648734 \ 5.90807e-79, 1.97286e-72, 3.08303e-67, 9.01922e-63, 8.51888e-59, 3.45948e-55, 7.19153e-52, 8.59756e-49, 6.42073e-46, 3.18522e-43, 1.10047e-40, 2.74878e-38, 5.11623e-36, 7.27457e-34, 8.06743e-32, 7.10181e-30, 5.03779e-28, 2.91735e-26, 1.39484e-24, 5.56103e-23, 1.86500e-21, 5.30232e-20, 1.28683e-18, 2.68249e-17, 4.82984e-16, 7.54890e-15, 1.02887e-13, 1.22788e-12, 1.28790e-11, 1.19130e-10, 9.74792e-10, 7.07586e-09, 4.56813e-08, 2.62910e-07, 1.35180e-06, 6.22152e-06, .0000256762, .0000951716, .000317292, .0009526922, .0025792733, .0063030003, .0139156652, .0277791274, .0501758127, .0820518274, .1215379868, .1631300305, .1984628503, .2188926296, .2188926296, .1984628503, .1631300305, .1215379868, .0820518274, .0501758127, .0277791274, .0139156652, .0063030003, .0025792733, .0009526922, .000317292, .0000951716, .0000256762, 6.22152e-06, 1.35180e-06, 2.62910e-07, 4.56813e-08, 7.07586e-09, 9.74792e-10, 1.19130e-10, 1.28790e-11, 1.22788e-12, 1.02887e-13, 7.54890e-15, 4.82984e-16, 2.68249e-17, 1.28683e-18, 5.30232e-20, 1.86500e-21, 5.56103e-23, 1.39484e-24, 2.91735e-26, 5.03779e-28, 7.10181e-30, 8.06743e-32, 7.27457e-34, 5.11623e-36, 2.74878e-38, 1.10047e-40, 3.18522e-43, 6.42073e-46, 8.59756e-49, 7.19153e-52, 3.45948e-55, 8.51888e-59, 9.01922e-63, 3.08303e-67, 1.97286e-72, 5.90807e-79) sum = rw[2, ]:*exp(rw[1, ]:^2):*sf(rw[1, ], arg) return(quadrowsum(sum)) } } real sf(real rowvector s, real vector arg) { return(s:*(1 :+ (s:^2)/(arg[1] + arg[2] - 1)):^(-(arg[1] + arg[2])/2)) } mata mlib create lmvtnorm, replace mata mlib add lmvtnorm *() mata mlib index end