*! Date : 03 Jan 2019 *! Version : 1.3 *! Authors : Michael J Grayling & Adrian P Mander /* 16/04/15 v1.0 Basic version complete. 08/09/17 v1.1 Added method option. 09/09/17 v1.2 Changed so df = 0 or df = . corresponds to mvnormalden. Added support for non-integer df. 03/01/19 v1.3 Converted delta and sigma to be optional with internal defaults. Removed method option for simplicity. */ program define rmvt, rclass version 15.0 syntax , [n(integer 1) DELta(numlist) Sigma(string) df(real 1)] ///// Check input variables //////////////////////////////////////////////////// local lendelta:list sizeof delta if (`n' < 1) { di "{error}Number of random vectors to generate (n) must be a strictly positive integer." exit(198) } if ("`sigma'" != "") { if (colsof(`sigma') != rowsof(`sigma')) { di "{error}Scale matrix Sigma (sigma) must be square." exit(198) } if ((`lendelta' != 0) & (`lendelta' != colsof(`sigma'))) { di "{error}Vector of non-centrality parameters (delta) 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) } } else { if (`lendelta' != 0) { mat chol = I(`lendelta') } else { mat chol = I(2) } } ///// Perform main computations //////////////////////////////////////////////// if (`lendelta' != 0) { local matadelta "" foreach l of local delta { if "`matadelta'" == "" local matadelta "`l'" else local matadelta "`matadelta',`l'" } mat delta = (`matadelta') } else if ("`sigma'" != "") { mat delta = J(1, colsof(`sigma'), 0) } else { mat delta = J(1, 2, 0) } mata: rmvt_void(`n', `df') ///// Output /////////////////////////////////////////////////////////////////// return mat rmvt = rmatrix end ///// Mata ///////////////////////////////////////////////////////////////////// mata: void rmvt_void(n, df) { delta = st_matrix("delta") C = st_matrix("chol") st_matrix("rmatrix", rmvt_mata(n, delta, ., df, C)) } 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)) } } end