*
* %LOGMVSKEWT(z,xi,nu)
* Log multivariate skew-t density function (with common degrees of
* freedom parameter).
*
* Parameters:
*   Z    a K vector, which is assumed to represent a standardized (mean 0,
*        covariance matrix I) process.
*   xi   a K vector of (positive) skewness parameters (xi(i)=1.0 means no
*        skewness, > 1 means skewed positive, < 1 means skewed negative).
*   nu   degrees of freedom for underlying Student t.
*
* If your model generates mean zero data with covariance matrix
* <<sigma>>, you need to standardize the input vector by pre-multiplying
* by inv(%decomp(sigma)) and adjust the likelihood by subtracting
* .5 * log |sigma|.
*
* Reference: Bauwens & Laurent(2005), "A New Class of Multivariate Skew
* Densities, With Application to Generalized Autoregressive Conditional
* Heteroscedasticity Models," JBES, vol 23, pp 346-354.
*
* Revision Schedule:
*  03/2012 Corrected formula and documented more extensively by Tom Doan,
*          Estima.
*
function %logmvskewt z xi nu
type real      %logmvskewt nu
type vector    z xi
local integer 	i k
local real    	lic zistar zistarsq mi si
*
compute k=%rows(z)
*
* Compute invariant parts of the integrating constant
*
compute lic=k*log(2.0/sqrt(%pi))+%lngamma(.5*(nu+k))-%lngamma(.5*nu)-.5*k*log(nu-2)
compute zistarsq=0.0
do i=1,k
   compute mi=%gamma(.5*(nu-1))*sqrt(nu-2)/(sqrt(%pi)*%gamma(.5*nu))*(xi(i)-1.0/xi(i))
   compute si=sqrt((xi(i)^2+1.0/xi(i)^2-1)-mi^2)
   compute zistar=(si*z(i)+mi)*%if(z(i)>=-mi/si,1.0/xi(i),xi(i))
   compute lic=lic+log(xi(i)*si/(1+xi(i)^2))
   compute zistarsq+=zistar^2
end do i
compute %logmvskewt=lic-.5*(k+nu)*log(1+zistarsq/(nu-2))
end
