Help needed please with State Space based regression model for balanced Panel Data

hello all, I am trying to fit a stochastic frontier model (not much different from linear regression) for a balanced panel data. I am allowing the inefficiency (u) to be correlated with the previous lag. I have tried the following code but I am getting very few neff samples and divergence.

Could someone please help?

data {
int<lower=1> N;// Number of observations
int<lower=1> J; // Number of groups
int<lower=1> P; // Number of predictors i
int<lower=1> T; // Number of time period
real Y[N]; // dependent variables
matrix[N,P] X; // matrix of independent variables☺
int<lower=1,upper=J> dhb_id[N];
int<lower=1,upper=T> TIME[N];
}

parameters {
real alpha;
real delta;
vector[P] beta;
real<lower=0,upper=1> rho;
real<lower=0> sigma;
real<lower=0> phi;
vector [N]u;
}

transformed parameters {
}

model {
alpha ~ normal(0,1);
beta ~ normal (0,1);
sigma ~ gamma(1,1);
delta ~ normal (0,1);
phi ~ gamma(1,1);
rho ~ beta(3,4);

for(n in 1:N) {
if(TIME[n]==1) {
u[n] ~ normal( delta/(1-rho), phi/(1-rho * rho)); /// impose stationary and intilize the model
}
else
u[n] ~ normal(delta + rho * u[n-1], phi);
}

Y ~ normal(alpha+ X*beta- exp(u) , sigma);

}

Please if someone could help me, I will be very gratefull.

This looks like the type of thing that could be non-centered: https://mc-stan.org/users/documentation/case-studies/divergences_and_bias.html

I think in models like these that amounts to making u a transformed parameter and writing the time series in terms of increments.

Check the reparameterization of h here: https://mc-stan.org/docs/2_25/stan-users-guide/stochastic-volatility-models.html, it is probably similar to what you’d end up doing for u in your model.

Hi Ben,
Thanks of the link. I have look through it. It is slightly different :)

Would you know or any one else on how to convert the following code into non centred parameterisation using transformed parameters?

for(n in 1:N) {
if(TIME[n]==1) {
target += normal_lpdf(u[n] | delta/(1-rho), phi /sqrt(1 - rho * rho));

}
else
target += normal_lpdf(u[n] | delta + rho*u[n-1], phi);

}

Well here’s a basic centered parameterization of theta:

parameters {
  real theta;
  real mu;
  real<lower = 0.0> phi;
}

model {
  mu ~ normal(0, 1);
  phi ~ normal(0, 1);
  theta ~ normal(mu, phi);
}

And here is the non-centered parameterization:

parameters {
  real z;
  real mu;
  real<lower = 0.0> phi;
}

transformed parameters {
  real theta = mu + phi * z;
}

model {
  mu ~ normal(0, 1);
  phi ~ normal(0, 1);
  z ~ normal(0, 1);
}

The point is in either case, theta takes on the prior distribution normal(mu, phi), just in the first it is sampled directly and in the second it is a transformed parameter.

In your case you have:

target += normal_lpdf(u[n] | delta + rho*u[n-1], phi);

which is like:

u[n] ~ normal(delta + rho*u[n-1], phi);

which if u[n] is a transformed parameter and z[n] is a normal(0, 1) then you could probably do something like:

u[n] = delta + rho * u[n + 1] + phi * z[n];

But check out the divergences case study for the details on the whys/hows of this: https://mc-stan.org/users/documentation/case-studies/divergences_and_bias.html

It takes a little time but is worth understanding cause it comes up a lot.