# Multivariate skew-t distribution

**URL:** <https://discourse.mc-stan.org/t/multivariate-skew-t-distribution/3583>\
**Category:** Modeling\
**Created:** [March 16, 2018, 4:28pm UTC](https://discourse.mc-stan.org/t/multivariate-skew-t-distribution/3583 "2018-03-16T16:28:22Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![KlausHoejlund](https://avatars.discourse-cdn.com/v4/letter/k/43a26b/32.png) [@KlausHoejlund](https://discourse.mc-stan.org/u/KlausHoejlund)\
**Post date:** [March 16, 2018, 4:28pm UTC](https://discourse.mc-stan.org/t/multivariate-skew-t-distribution/3583/1 "2018-03-16T16:28:22Z")

</div>

Hi. I am trying to incorporate non-gaussian returns in asset-allocation models. In regards to that, i need to make draws from a multivariate skew t-distribution. It is based on Sahu, Dey & Branco (2003) and used in Jondau & Rockinger (2012), and I have attached an image. I have the priors of the parameters, but the problem is how to draw from the multivariate skew-t distribution shown in equation 4 on the picture. I don’t think that it’s implemented in Stan.

Can one of you explain to me how this can be implemented?

 ![Udklip](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/d/d68ccbaff36b339cc009323ddd282386ff75a021.PNG)

---

<div class="post-metadata">

**Author:** ![bgoodri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bgoodri/32/4451_2.png) [@bgoodri](https://discourse.mc-stan.org/u/bgoodri)\
**Post date:** [March 16, 2018, 4:49pm UTC](https://discourse.mc-stan.org/t/multivariate-skew-t-distribution/3583/2 "2018-03-16T16:49:14Z")

</div>

There is no multivariate skew-t distribution in the Stan Math library, but if equation (4) is supposed to represent a PDF, then just work out its logarithm analytically and implement that as a user-defined function in the Stan language.

---

<div class="post-metadata">

**Author:** ![KlausHoejlund](https://avatars.discourse-cdn.com/v4/letter/k/43a26b/32.png) [@KlausHoejlund](https://discourse.mc-stan.org/u/KlausHoejlund)\
**Post date:** [March 16, 2018, 5:03pm UTC](https://discourse.mc-stan.org/t/multivariate-skew-t-distribution/3583/3 "2018-03-16T17:03:05Z")

</div>

Thank you very much. Can you elaborate a bit more on what you mean by: “then just work out its logarithm analytically”?

---

<div class="post-metadata">

**Author:** ![bgoodri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bgoodri/32/4451_2.png) [@bgoodri](https://discourse.mc-stan.org/u/bgoodri)\
**Post date:** [March 16, 2018, 5:22pm UTC](https://discourse.mc-stan.org/t/multivariate-skew-t-distribution/3583/4 "2018-03-16T17:22:10Z")

</div>

Stan uses the logarithm of PDFs and PMFs for numerical and algorithmic reasons. Equation (4) appears to be a PDF that is the product of conditionally independent observations at time `t`. The logarithm of a product is the sum of the logarithms, the logarithm of a ratio is the difference between the logarithm of the numerator and the logarithm of the denominator, and the logarithm of a power is equal to the exponent multiplied by the logarithm of the base. Also, the logarithm of the Gamma function is denoted `lgamma` in Stan and lots of other computer languages. With those rules, you can work out what is the logarithm of the right-hand side of (4) and implement it as a user-defined function in the Stan language in a numerically stable way.

---

<div class="post-metadata">

**Author:** ![Rva92](https://avatars.discourse-cdn.com/v4/letter/r/aca169/32.png) [@Rva92](https://discourse.mc-stan.org/u/Rva92)\
**Post date:** [March 18, 2018, 12:48pm UTC](https://discourse.mc-stan.org/t/multivariate-skew-t-distribution/3583/5 "2018-03-18T12:48:06Z")

</div>

Thanks bgoordi, that is a great solution. I would like to check the function which I have written, is there a way to generate a random sample from the specified function for some parameter of nu and xi, which I then can plot in R to check, whether the distribution has the desired shape?

I have tried writting the function as:

functions {  
//Multivariate skewed t-distribution (Sahu, Dey & Branco, 2003)  
real MskT\_lpdf(real z, real nu, real xi){  
real a = (tgamma((nu-1)/2)_sqrt(nu-2))/(sqrt(pi())tgamma(nu/2))(xi-1/xi);  
real b = sqrt(pow(xi,2)+1/pow(xi,2)-1-pow(a,2));  
real k = (-a/b\<= z ? (b_z+a)_xi : (b_z+a)/xi);  
return log(2_2/(xi+1/xi)) + lgamma((nu+1)/2)-(log(sqrt(pi()_(nu-2)))+lgamma(nu/2))-((nu+1)/2)\*log(1+(pow(k,2)/(nu-2)));  
}  
}

data{  
real nu;  
real xi;  
}

model{  
real m~MskT\_lpdf(nu = nu, xi = xi);  
}

I.e. how can i get the stan program when called for R, to return a random sample drawn from the MskT distribution?

---

<div class="post-metadata">

**Author:** ![bgoodri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bgoodri/32/4451_2.png) [@bgoodri](https://discourse.mc-stan.org/u/bgoodri)\
**Post date:** [March 18, 2018, 2:27pm UTC](https://discourse.mc-stan.org/t/multivariate-skew-t-distribution/3583/6 "2018-03-18T14:27:50Z")

</div>

That function seems to be univariate but perhaps it is a valid special case. In any event, you need to call `expose_stan_functions` and verify that the antilog of MskT\_lpdf integrates to 1. See  
[http://mc-stan.org/users/documentation/case-studies/gpareto\_functions.html](http://mc-stan.org/users/documentation/case-studies/gpareto_functions.html)

---

<div class="post-metadata">

**Author:** ![Rva92](https://avatars.discourse-cdn.com/v4/letter/r/aca169/32.png) [@Rva92](https://discourse.mc-stan.org/u/Rva92)\
**Post date:** [March 18, 2018, 3:04pm UTC](https://discourse.mc-stan.org/t/multivariate-skew-t-distribution/3583/7 "2018-03-18T15:04:07Z")

</div>

After 2 korrektion it did indeed integrate to one - that is a nice check - but is there some way, that I in a similar manner can draw e.g. 1000 observations from the distribution and plot them as a histogram to check that it has the sought shape? Either directly in stan or indirectly in R?

---

<div class="post-metadata">

**Author:** ![sid](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/sid/32/5082_2.png) [@sid](https://discourse.mc-stan.org/u/sid)\
**Post date:** [March 18, 2018, 3:16pm UTC](https://discourse.mc-stan.org/t/multivariate-skew-t-distribution/3583/8 "2018-03-18T15:16:58Z")

</div>

Some ideas here as well: [https://stackoverflow.com/questions/43860236/multivariate-skew-normal-in-stan](https://stackoverflow.com/questions/43860236/multivariate-skew-normal-in-stan)

---

<div class="post-metadata">

**Author:** ![bgoodri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bgoodri/32/4451_2.png) [@bgoodri](https://discourse.mc-stan.org/u/bgoodri)\
**Post date:** [March 18, 2018, 4:08pm UTC](https://discourse.mc-stan.org/t/multivariate-skew-t-distribution/3583/9 "2018-03-18T16:08:25Z")

</div>

You would need to write another function in the Stan language (or R) to draw from a multivariate skewed t distribution.

---

<div class="post-metadata">

**Author:** ![Rva92](https://avatars.discourse-cdn.com/v4/letter/r/aca169/32.png) [@Rva92](https://discourse.mc-stan.org/u/Rva92)\
**Post date:** [March 18, 2018, 4:37pm UTC](https://discourse.mc-stan.org/t/multivariate-skew-t-distribution/3583/10 "2018-03-18T16:37:20Z")

</div>

Thank you very much. I am now trying to convert it to the multivariate case - the one stated by Klaus. However, I cannot get it to loop correctly - how should I declare the parameter vectors (nu and xi) and the and the random variabel, z, such that it matches with the number of random variables in the vector z?

functions {  
//Multivariate skewed t-distribution (Sahu, Dey & Branco, 2003)  
real MskT\_lpdf(vector z, vector nu, vector xi){  
vector uni\_pdf;  
for (i in 1:rows(z){  
real a = (tgamma((nu[i]-1)/2)_sqrt(nu[i]-2))/(sqrt(pi())tgamma(nu[i]/2))(xi[i]-1/xi);  
real b = sqrt(pow(xi[i],2)+1/pow(xi[i],2)-1-pow(a,2));  
real k = (-a/b \>= z[i] ? (b_z+a)_xi[i] : (b_z+a)/xi[i]);  
uni\_pdf[i] = log(2_b/(xi[i]+1/xi[i])) + lgamma((nu[i]+1)/2)-(log(sqrt(pi()_(nu[i]-2)))+lgamma(nu[i]/2))-((nu[i]+1)/2)\*log(1+(pow(k,2)/(nu[i]-2)));  
}  
return sum(uni\_pdf)  
}  
}

I apologize for the many questions and hope You can help :)

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [March 19, 2018, 6:05am UTC](https://discourse.mc-stan.org/t/multivariate-skew-t-distribution/3583/11 "2018-03-19T06:05:06Z")

</div>

> [@Rva92](#):
>
> is there some way, that I in a similar manner can draw e.g. 1000 observations from the distribution and plot them as a histogram

Yes, just define a parameter and draw from it:

```no-highlight
functions {
  real foo_lpdf(...) { ... }
}
data {
  int<lower = 0> K;
  ...
}
parameters {
  vector[K] mu;
}
model {
  mu ~ foo(...);
}

```

That’s not as efficient as writing an independent random number generator, but it will test that your function is defined properly.

---

<div class="post-metadata">

**Author:** ![Rva92](https://avatars.discourse-cdn.com/v4/letter/r/aca169/32.png) [@Rva92](https://discourse.mc-stan.org/u/Rva92)\
**Post date:** [March 19, 2018, 5:40pm UTC](https://discourse.mc-stan.org/t/multivariate-skew-t-distribution/3583/12 "2018-03-19T17:40:10Z")

</div>

Thank you so much - it works perfectly! Awesome program and awesome support!

---

<div class="post-metadata">

**Author:** ![Rva92](https://avatars.discourse-cdn.com/v4/letter/r/aca169/32.png) [@Rva92](https://discourse.mc-stan.org/u/Rva92)\
**Post date:** [March 21, 2018, 7:53am UTC](https://discourse.mc-stan.org/t/multivariate-skew-t-distribution/3583/13 "2018-03-21T07:53:22Z")

</div>

Now I am trying to run the multivariate case:  
//This code represents our model written in stan (C++)

functions {  
//Multivariate skewed t-distribution (Sahu, Dey & Branco, 2003)  
real MskT\_lpdf(vector z, vector nu, vector xi){  
int N = rows(nu);  
vector[K] l;  
for(n in 1:K){  
real a = (tgamma((nu[n]-1)/2)_sqrt(nu[n]-2))/(sqrt(pi())tgamma(nu[n]/2))(xi[n]-1/xi[n]); // To ensure zero mean  
real b = sqrt(pow(xi[n],2)+1/pow(xi[n],2)-1-pow(a,2)); // To ensure unit variance  
real k = (-a/b \>= z[n] ? (b_z[n]+a)_xi[n] : (b_z[n]+a)/xi[n]);  
l[n] = log(2_b/(xi[n]+1/xi[n])) + lgamma((nu[n]+1)/2)-(log(sqrt(pi()_(nu[n]-2)))+lgamma(nu[n]/2))-((nu[n]+1)/2)\*log(1+(pow(k,2)/(nu[n]-2)));  
}  
return sum(l);  
}  
}

However, when I try running it from R by checking whether it integrates to 1, using:  
rstan::expose\_stan\_functions(“C:/Users/rva92/Dropbox/Mine ting/Uni/10. Semester (Finance master thesis)/4 - Kodning/Stan functions/Multivariate skewed t distribution.stan”)

skw\_t \<- function(z, nu, xi, K) {  
exp(sapply(z, FUN = MskT\_lpdf, nu = nu, xi = xi))  
}

xi = c(4, 2, 0.5, 3)  
nu = c(4, 5, 3, 4)  
K = as.numeric(length(nu))  
integrate(skw\_t, lower = -Inf, upper = Inf, nu = nu, xi = xi, K = K)

I get the following message:  
Error in FUN(X[[i]], …) :  
[]: accessing element out of range. index 2 out of range; expecting index to be between 1 and 1; index position = 1z

Any idea on, where it went wrong?

---

<div class="post-metadata">

**Author:** ![jmh530](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jmh530/32/4579_2.png) [@jmh530](https://discourse.mc-stan.org/u/jmh530)\
**Post date:** [March 22, 2018, 9:32pm UTC](https://discourse.mc-stan.org/t/multivariate-skew-t-distribution/3583/14 "2018-03-22T21:32:30Z")

</div>

Multivariate t is only appropriate in very limited situations. That is, where you have reason to believe that the degrees of freedom parameter is the same for every asset, both on the marginal distribution and on the joint distribution. In the classical world, a grouped T distribution was proposed as a solution, but it seems like copula modelling has become more important.

You might consider instead a factor structure, such as CAPM or APT. Assuming the factors are uncorrelated, then you just need to fit a univariate t distribution to the factors. The joint distribution may have some non-normal distribution that way, even if the returns, conditional on the factor, are assumed to be normal.

Another issue with multivariate t is that it is a little more tricky to set up a time-varying version. You’re sort of forced to choose whether the multivariate dependence being t is more important than the time-varying volatility. Another reason to prefer a factor structure of some kind. A little more flexibility.