Hi, I am trying to estimate the following model:

Where g(.|υt,ηt) is a skewed t distribution (see code) and υ is the degrees of freedom, η is the asymmetry parameter. Σt^(1/2) is the Cholesky decomposition of the variance-covariance matrix – but this is not important in the one asset case. The variance of the asset is estimated via a GJR GARCH model, while the degrees of freedom and asymmetry parameter are also estimated using GARCH like models – they are stated as:

I have tried by implementing the following code - it runs, but it cannot start the iteration correctly, i.e. it cannot find appropriate starting values - I have tried also incorporating initial value as suggest by Stan, but it does not fix it. Any suggestions on, how to optimize the code?

```
functions {
//Multivariate skewed t-distribution (Sahu, Dey & Branco, 2003)
real MskT_lpdf(real z, real nu, real xi){
real a = exp((lgamma((nu-1)/2)+log(sqrt(nu-2))) - (log(sqrt(pi()))+lgamma(nu/2))+log((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*b/(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 {
int<lower=0> T; //number of observations
real r[T]; //return data - notice, cannot med -100%, as this will results in negative prices
real<lower=0> sigma1; //First sigma
//GARCH prior parameters
real omg_p;
real omg_q;
real psi_p;
real psi_q;
real alpha_q;
real alpha_p;
real lam_p;
real lam_q;
real mu_mu;
real mu_sigma;
//evelution of nu and xi parameters
real c0_p;
real c0_q;
real c1_plus_p;
real c1_plus_q;
real c1_minus_p;
real c1_minus_q;
real c2_p;
real c2_q;
real d0_p;
real d0_q;
real d1_plus_p;
real d1_plus_q;
real d1_minus_p;
real d1_minus_q;
real d2_p;
real d2_q;
}
parameters {
real mu;
real<lower=0> omg; //Cannot be negative to ensure positive variance
real<lower=0,upper=1> alpha; //Cannot be negative to ensure positive variance
real<lower=0> psi; //Cannot be negative to ensure positive variance
real<lower=0> lam; //Cannot be negative to ensure positive variance
real c0;
real c1_minus;
real c1_plus;
real c2;
real d0;
real d1_minus;
real d1_plus;
real d2;
}
transformed parameters{
real<lower=0> sigma[T];
real<lower=0, upper=(1-alpha-psi/2)> bet; //Upper limit to ensure stationarity
real nu[T];
real xi[T];
real inn[T];
real eps[T];
sigma[1] = sigma1;
bet = lam - alpha - psi/2;
inn[1] = 0.2;
nu[1] = 4;
xi[1] = 0.5;
for (t in 2:T){
sigma[t] = sqrt(omg + (alpha + psi*step(mu-r[t - 1]))*pow(r[t - 1] - mu,2) + bet*pow(sigma[t-1],2));
nu[t] = exp(log(c0)+c1_minus*fabs(inn[t-1])*step(-inn[t-1])+c1_plus*fabs(inn[t-1])*(1-step(-inn[t-1]))+c2*log(nu[t-1]-4))-4;
xi[t] = exp(log(d0)+d1_minus*inn[t-1]*step(-inn[t-1])+d1_plus*inn[t-1]*(1-step(-inn[t-1]))+d2*log(xi[t-1]));
inn[t] = (r[t]-mu) / sigma[t];
}
}
model{
//Data distribution - model
for(i in 2:T){
inn[i] ~ MskT(nu[i], xi[i]);
}
//prior distributions
omg~beta(omg_p, omg_q);
alpha~beta(psi_p, psi_q);
psi~beta(alpha_q, alpha_q);
lam~beta(lam_p, lam_q);
mu~normal(mu_mu, mu_sigma);
c0~lognormal(c0_p, c0_q);
c1_minus~normal(c1_minus_p, c1_minus_q);
c1_plus~normal(c1_plus_p, c1_plus_q);
c2~normal(c2_p, c2_q)T[-1,1];
d0~lognormal(d0_p, d0_q);
d1_minus~normal(d1_minus_p, d1_minus_q);
d1_plus~normal(d1_plus_p, d1_plus_q);
d2~normal(d2_p, d2_q)T[-1,1];
}
```

[edit: escaped code]