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