Hi,
First, let me say that I am a beginner with Stan and I find it amazing both from an ease of use and documentation perspective.
I am trying to fit a 4 state regime switching model to some market data, and I am getting divergent transitions and very low n_eff values. The diagnostic charts in ShinyStan show that the chains are not mising well. I read the Workflow guide below and I understand that I need to change the parametrisation in my model, but I don’t really understand why the new parametrisation fixes the issue in the example.
https://betanalpha.github.io/assets/case_studies/rstan_workflow.html
I would welcome any comment on my implementation of the model, and any suggestion on how to fix my divergence issues.
To briefly describe my model, the value of r_t is distributed according to \mathcal{N}( \mu_i, \sigma_i), with i \in \{1,2,3,4\}. The rows of the transition matrix P are distributed according to a dirichlet distribution. Finally, there are constraints on the signs of the \mu_i. \Pi is the vector of unconditional probability of each state, and I try to impose constraints on average returns in states 1,2 and 3,4.
data {
int<lower=1> T; // number of observations (length) - time series of returns
real r[T]; // weekly return
// 1-2 bear regime
// 3-4 bull regime
real m[4]; // mu prior mean m = c(-0.7, 0.2, -0.2,0.3)
real n[4]; // mu prior sd n = c(1, 1, 1, 1)
real v[4]; // sigma2 prior v = c(1,1,1,1)
real s[4]; // sigma2 prior s = c(1,1,1,1)
vector<lower=0>[3] alpha[4]; // transit prior rows
// {p11, p12, p14} ~ Dir(8, 1.5, 0.5),
// {p21, p22, p24} ~ Dir(1.5, 8, 0.5),
// {p31, p33, p34} ~ Dir(0.5, 8, 1.5),
// {p41, p43, p44} ~ Dir(0.5, 1.5, 8).
simplex[4] s0; // inital state
}
parameters {
simplex[3] reppp[4]; // transit probs row i
real<upper=0> mu1;// bear market
real<lower=0> mu2;// bear rally
real<upper=0> mu3;// bull correction
real<lower=0> mu4;// bull market
real<lower=0> sigma2[4]; // sd
// simplex[4] s0; // inital state
}
transformed parameters {
vector[4] eta[T]; // First, let’s define the likelihood contribution of an observation r[t] under each state
simplex[4] si[T]; // probability of being in each state
vector[T] f; // likelihood
vector[4] PI; // unconditional state probs
matrix[4,4] P; // transit probs
real mu_bull;
real mu_bear;
// put together the transition matrix
for (j in 1:4){
P[j, 1] = reppp[j,1];
if ( j == 1 || j==2 ){
P[j, 2] = reppp[j,2];
P[j, 3] = 0;
} else {
P[j, 2] = 0;
P[j, 3] = reppp[j,2];
}
P[j, 4] = reppp[j,3];
}
// fill in etas (likelihood contribution of an observation r[t] under each state)
{
real mu[4]; // expected returns
mu[1] = mu1;
mu[2] = mu2;
mu[3] = mu3;
mu[4] = mu4;
for(t in 1:T) {
for(j in 1:4 ){
eta[t,j] = exp(normal_lpdf(r[t] | mu[j], sigma2[j]));
}
}
}
// work out likelihood contributions
f[1] = dot_product( P*s0, eta[1]);
si[1] = ((P*s0) .* eta[1]) / f[1];
// and for the rest
for(t in 2:T) {
f[t] = dot_product( P*si[t-1], eta[t]);
si[t] = ((P*si[t-1]) .* eta[t]) / f[t];
}
// unconditional probabilities PI
{
matrix[4,5] A;
vector[5] ee;
matrix[4,4] AAprime;
vector[4] AE;
A = append_col(P - diag_matrix(rep_vector(1., 4)),rep_vector(1., 4));
ee = rep_vector(0., 5);
ee[5] = 1.0;
// PI = (A' * A) \ (A' * ee);
AAprime = tcrossprod(A);
AE = A * ee;
PI = AAprime \ AE;
mu_bear = PI[1]/(PI[1]+PI[2])*mu1 + PI[2]/(PI[1]+PI[2])*mu2;
mu_bull = PI[3]/(PI[3]+PI[4])*mu3 + PI[4]/(PI[3]+PI[4])*mu4;
}
}
// T observations (length)
model {
for (kk in 1:4){
reppp[kk] ~ dirichlet(alpha[kk]); // Sample the transition probability
sigma2[kk] ~ cauchy( v[kk], s[kk] );
}
mu1 ~ normal( m[1], n[1] );
mu2 ~ normal( m[2], n[2] );
mu3 ~ normal( m[3], n[3] );
mu4 ~ normal( m[4], n[4] );
// likelihood is really easy here!
target += sum(log(f));
if( mu_bull < 0){ target += negative_infinity(); }
if( mu_bear > 0){ target += negative_infinity(); }
}