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] ? (bz[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?