# Dynamic panel data models with Stan?

**URL:** <https://discourse.mc-stan.org/t/dynamic-panel-data-models-with-stan/5136>\
**Category:** Modeling\
**Created:** [August 9, 2018, 2:09pm UTC](https://discourse.mc-stan.org/t/dynamic-panel-data-models-with-stan/5136 "2018-08-09T14:09:06Z")\
**Posts on this page:** 20\
**Page:** 2

<div class="post-metadata">

**Author:** ![c5v](https://avatars.discourse-cdn.com/v4/letter/c/e19adc/32.png) [@c5v](https://discourse.mc-stan.org/u/c5v)\
**Post date:** [January 31, 2019, 7:27pm UTC](https://discourse.mc-stan.org/t/dynamic-panel-data-models-with-stan/5136/22 "2019-01-31T19:27:23Z")

</div>

Thank you so much for your explanation! I am going to try this approach.

---

<div class="post-metadata">

**Author:** ![rtrangucci](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/rtrangucci/32/8375_2.png) [@rtrangucci](https://discourse.mc-stan.org/u/rtrangucci)\
**Post date:** [January 31, 2019, 7:28pm UTC](https://discourse.mc-stan.org/t/dynamic-panel-data-models-with-stan/5136/23 "2019-01-31T19:28:53Z")

</div>

Great! Happy to help. Let me know how it goes.

---

<div class="post-metadata">

**Author:** ![c5v](https://avatars.discourse-cdn.com/v4/letter/c/e19adc/32.png) [@c5v](https://discourse.mc-stan.org/u/c5v)\
**Post date:** [February 1, 2019, 12:24pm UTC](https://discourse.mc-stan.org/t/dynamic-panel-data-models-with-stan/5136/24 "2019-02-01T12:24:31Z")

</div>

@rtrangucci Sorry for disturbing you again. But I have been trying to get your approach working in my multilevel context.

My stan code now looks like the following:

```stan
data {
  int<lower=1> I;
  int<lower=1> T;
  matrix[T,I] y;
  matrix[T,I] y_lag;
  
  // Define variables in data test
  // Number of level-1 observations (an integer)
  int<lower=0> Ni;
  // Number of level-2 clusters
  int<lower=0> Nj;
  // Number of level-3 clusters
  int<lower=0> Nk;
  // Numer of models x countries
  int<lower = 0> Njk;
  // Number of fixed effect parameters
  int<lower=0> p;
  // Number of random effect parameters
  //int<lower=0> q;
  int<lower=0> Npars;
  
  // Variables with fixed coefficient
  matrix[Ni,p] fixedEffectVars;
  
  // Cluster IDs
  int<lower=1> modelid[Ni];
  int<lower=1> countryid[Ni];
  
  // Level 3 look up vector for level 2
  int<lower=1> countryLookup[Npars];
  int<lower=1> countryMonthLookup[Njk];
}
parameters {
  real<lower=0,upper=1> delta_raw;
  //real<lower=0> sigma_u;
  real<lower=0> sigma_e;
  //vector[I] u_raw;
  
    // Define parameters to estimate
  // Population intercept (a real number)
  real beta_0;
  
  // Fixed effects
  vector[p] beta;
  
  
  // Level-2 random effect
  real u_0jk[Npars];
  real<lower=0> sigma_u0jk;
  
  // Level-3 random effect
  real u_0k[Nk];
  real<lower=0> sigma_u0k;
}
transformed parameters {
  //vector[I] u = sigma_u * u_raw;
  real delta = delta_raw * 2 - 1;
  
    // Varying intercepts
  real beta_0jk[Npars];
  real beta_0k[Nk];
  
  // Individual mean
  real mu[Ni];

  // Varying intercepts definition
  // Level-3 (10 level-3 random intercepts)
  for (k in 1:Nk) {
    beta_0k[k] = beta_0 + u_0k[k];
  }
  // Level-2 (100 level-2 random intercepts)
  for (j in 1:Npars) {
    beta_0jk[j] = beta_0k[countryLookup[j]] + u_0jk[j];
  }
  
    //Individual mean
  for (i in 1:Ni) {
    mu[i] = beta_0jk[countryMonthLookup[i]] + fixedEffectVars[i,]*beta;
  }
}
model {
  real one_minus_delta_sq = (1 - delta) * (1 + delta);
  int count = 1;
  //vector[I] u_delta = u * (1 - delta);
  for (i in 1:I) {
    // target += normal_lpdf(y[1,i] | u[i],
    // sigma_e / sqrt(one_minus_delta_sq));
    // target += normal_lpdf(y[2:T,i] | u_delta[i]
    // + delta * y_lag[2:T,i], sigma_e);
    for(j in 1:T){
      if(j == 1){
        y[j,i] ~ normal(mu[count], sigma_e / sqrt(one_minus_delta_sq));
      }
      else{
        y[j,i] ~ normal(mu[count] + delta * y_lag[j,i] - delta*mu[count], sigma_e);
      }
      count += 1;
    }
  }
  //u_raw ~ normal(0, 1);
  //sigma_u ~ normal(0, 2);
  sigma_e ~ normal(0, 1);
  delta_raw ~ beta(4, 4);
  
  // Random effects distribution
  u_0k ~ normal(0, sigma_u0k);
  u_0jk ~ normal(0, sigma_u0jk);
  beta[1] ~ normal(-0.25, 1);
}
generated quantities{
    matrix[T,I] y_rep;
    int count_rep = 1;
    real one_minus_delta_sq = (1 - delta) * (1 + delta);
    //vector[I] u_delta = u * (1 - delta);

    for (i in 1:I) {
      for(j in 1:T){
        if(j == 1){
          y_rep[j,i] = normal_rng(mu[count_rep], sigma_e / sqrt(one_minus_delta_sq));
        }
        else{
          y_rep[j,i] = normal_rng(mu[count_rep] + delta * y_lag[j,i] - delta*mu[count_rep], sigma_e);
        }
        count_rep += 1;
      }
    }
}

```

Here `mu[count]` denotes the \mu\_{ij} from my previous post.

When I estimate this model I get the warning that all my transitions after warm-up are divergent. However, all my parameters seem to have a reasonable value and a Rhat of 1, except from the `delta` parameter which is 0.00. This is very strange to me and I don’t really know what is going wrong here. Do you have any idea?

---

<div class="post-metadata">

**Author:** ![rtrangucci](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/rtrangucci/32/8375_2.png) [@rtrangucci](https://discourse.mc-stan.org/u/rtrangucci)\
**Post date:** [February 1, 2019, 4:44pm UTC](https://discourse.mc-stan.org/t/dynamic-panel-data-models-with-stan/5136/25 "2019-02-01T16:44:04Z")

</div>

This is a pretty complex model, so it’s hard for me to see what the problem is. Does the model fit OK without the lagged term? Have you done fake data checking to make sure you can pull out known parameters in a model that doesn’t include the lagged term? If yes, can you fit fake data that includes the AR1 term?

---

<div class="post-metadata">

**Author:** ![c5v](https://avatars.discourse-cdn.com/v4/letter/c/e19adc/32.png) [@c5v](https://discourse.mc-stan.org/u/c5v)\
**Post date:** [February 1, 2019, 7:19pm UTC](https://discourse.mc-stan.org/t/dynamic-panel-data-models-with-stan/5136/26 "2019-02-01T19:19:31Z")

</div>

Yes, without the lagged term the model works perfectly fine and predictive accuracy of the model is already pretty good. So the error only occurs when I add the lagged term in:

```stan
if(j == 1){
          y_rep[j,i] = normal_rng(mu[count_rep], sigma_e / sqrt(one_minus_delta_sq));
        }
        else{
          y_rep[j,i] = normal_rng(mu[count_rep] + delta * y_lag[j,i] - delta*mu[count_rep], sigma_e);
}

```

The parameter estimates are still fine after this model and predictive accuracy is basically the same as without the lagged term. The weird thing is that the `delta` parameter equals 0 and I get all those divergent transitions.

---

<div class="post-metadata">

**Author:** ![rtrangucci](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/rtrangucci/32/8375_2.png) [@rtrangucci](https://discourse.mc-stan.org/u/rtrangucci)\
**Post date:** [February 1, 2019, 7:47pm UTC](https://discourse.mc-stan.org/t/dynamic-panel-data-models-with-stan/5136/27 "2019-02-01T19:47:51Z")

</div>

Have you tried non-centering your random effects `u_0k` and `u_0jk`?

Edit: noncentering would involve declaring a new parameter, `u_0k_raw`, e.g. and declaring `u_0k` as a transformed parameter like so:

```
parameters {
  ...
  vector[Nk] u_0k = sigma_u0k * u_0k_raw;
  ...
}
model {
...
u_0k_raw ~ normal(0, 1);
...
}

```

---

<div class="post-metadata">

**Author:** ![c5v](https://avatars.discourse-cdn.com/v4/letter/c/e19adc/32.png) [@c5v](https://discourse.mc-stan.org/u/c5v)\
**Post date:** [February 1, 2019, 10:04pm UTC](https://discourse.mc-stan.org/t/dynamic-panel-data-models-with-stan/5136/28 "2019-02-01T22:04:08Z")

</div>

Unfortunately, this gives me again only divergent transitions after warm-up. Next to that, the `delta` parameter stays 0

---

<div class="post-metadata">

**Author:** ![rtrangucci](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/rtrangucci/32/8375_2.png) [@rtrangucci](https://discourse.mc-stan.org/u/rtrangucci)\
**Post date:** [February 1, 2019, 11:02pm UTC](https://discourse.mc-stan.org/t/dynamic-panel-data-models-with-stan/5136/29 "2019-02-01T23:02:08Z")

</div>

Did you generate fake data from the AR1 model and fit the fake data, or is this fitted to real data? If not, I’d say you should generate fake data from the AR1 model and make sure you can pull out the right parameter values from the fake data before fitting it to the real data.

---

<div class="post-metadata">

**Author:** ![c5v](https://avatars.discourse-cdn.com/v4/letter/c/e19adc/32.png) [@c5v](https://discourse.mc-stan.org/u/c5v)\
**Post date:** [February 3, 2019, 4:47pm UTC](https://discourse.mc-stan.org/t/dynamic-panel-data-models-with-stan/5136/30 "2019-02-03T16:47:29Z")

</div>

@rtrangucci sorry for my late reply, but I have been trying to get the model working for a fake simulated dataset.

I simulated data according to y\_{it} = (\alpha + u\_{i}) + x\_{i}\beta +\delta y\_{i,t-1} +\epsilon\_{it} using this code:

```
library(dplyr)
set.seed(123)
delta <- 0.8
sigma_u <- 0.5
sigma_e <- 1
I <- 100
T <- 30
u = rnorm(I, 0, sigma_u)
beta = 10
X_init = rnorm(I, 0, 1)#matrix(rnorm(I*T), nrow = T, ncol = I)
Xbeta_init = as.matrix(X_init)%*%beta
initial_values <- data_frame(
  individual = 1:I,
  u = u,
  Xbeta = as.vector(Xbeta_init),
  initial_y = u + as.vector(Xbeta_init) + sigma_e / 
    sqrt(1 - delta^2) * rnorm(I) # p(y_1 | u, sigma_e, delta)
)

simulated_panel <- initial_values %>%
  group_by(individual) %>% 
  do({
    simulate <- function(x) {
      out <- data.frame(y = rep(NA, T),
                        time = rep(NA, T),
                        u = x$u)
      for(t in 1:T) {
        out$time[t] <- t
        if(t == 1) {
          out$y[t] <- x$initial_y
        } else {
          out$y[t] <- rnorm(1, delta * (out$y[t-1] - x$u - x$Xbeta) + x$u + x$Xbeta, sigma_e)
        }
      }
      return(out)
    }
    as.data.frame(simulate(.))
  }) %>% 
  group_by(individual) %>% 
  mutate(lagged_y = lag(y)) %>% ungroup()

y <- matrix(0, T, I)
y_lag <- matrix(0, T, I)
for (i in 1:I) {
  y[,i] <- filter(simulated_panel, individual == i)$y
  y_lag[2:T,i] <- filter(simulated_panel, individual == i)$lagged_y[2:T]
}

```

Afterwards, I estimated my model which looks like this in stan:

```stan
data {
  int<lower=1> I;
  int<lower=1> T;
  matrix[T,I] y;
  matrix[T,I] y_lag;
  
  int<lower=0> Ni;
  
  // Variables with fixed coefficient
  matrix[Ni,p] fixedEffectVars;
}
parameters {
  real<lower=0,upper=1> delta_raw;
  real<lower=0> sigma_u;
  real<lower=0> sigma_e;
  vector[I] u_raw;
  
  // Fixed effects
  vector[p] beta;
}
transformed parameters {
  vector[I] u = sigma_u * u_raw;
  real delta = delta_raw * 2 - 1;
 
  //Individual mean
  for (i in 1:Ni) {
    mu[i] = fixedEffectVars[i,]*beta;// + beta_0jk[countryMonthLookup[i]];
  }
}
model {
  real one_minus_delta_sq = (1 - delta) * (1 + delta);
 
  for (i in 1:I) {
    for(j in 1:T){
      if(j == 1){
        y[j,i] ~ normal(u[i] + mu[i], sigma_e / sqrt(one_minus_delta_sq));
      }
      else{
        y[j,i] ~ normal(u[i] + mu[i] + delta * y_lag[j,i] - delta*(u[i]+mu[i]), sigma_e);
      }
    }
  }
  
  u_raw ~ normal(0, 1);
  sigma_u ~ normal(0, 2);
  sigma_e ~ normal(0, 1);
  delta_raw ~ beta(4, 4);
  
}
generated quantities{
  matrix[T,I] y_rep;
  real one_minus_delta_sq = (1 - delta) * (1 + delta);
  
  for (i in 1:I) {
    for(j in 1:T){
      if(j == 1){
        y_rep[j,i] = normal_rng(u[i] + mu[i], sigma_e / sqrt(one_minus_delta_sq));
      }
      else{
        y_rep[j,i] = normal_rng(u[i] + mu[i] + delta * y_lag[j,i] - delta*(u[i]+mu[i]), sigma_e);
      }
    }
  }
}

```

Using

```
stan_dat <- list(I = I, T = T, y = y, y_lag = y_lag, Ni = 100, p = 1, fixedEffectVars = Xbeta_init)
library(rstan)

fit <- stan(file = "dynamic_test.stan",
            data = stan_dat, iter = 5000, warmup = 2500, chains = 4)

```

In this case the `delta`, `sigma_u` and `sigma_e` parameters are estimated correctly. However, the estimated `beta`parameter equals 1, while it is equal to 10. Also, when changing the `beta` parameter to a different value in the DGP, the estimated value still equals about 1. Next to that I get the warning message that there were 2500 transitions after warmup that exceeded the maximum treedepth. Do you have any idea why this is the case and how to obtain the correct `beta` estimates?

---

<div class="post-metadata">

**Author:** ![rtrangucci](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/rtrangucci/32/8375_2.png) [@rtrangucci](https://discourse.mc-stan.org/u/rtrangucci)\
**Post date:** [February 4, 2019, 6:21pm UTC](https://discourse.mc-stan.org/t/dynamic-panel-data-models-with-stan/5136/31 "2019-02-04T18:21:40Z")

</div>

Sorry, I’m having trouble following your logic in the Stan program with `Ni`, but I put together an example with fixed effects.

```
library(dplyr)
library(rstan)
set.seed(123)
delta <- 0.8
sigma_u <- 0.5
sigma_e <- 1
I <- 100
T <- 30
u = rnorm(I, 0, sigma_u)
p <- 2
beta = rnorm(p)
X <- lapply(1:I, function(x) matrix(rnorm(p * T), nrow = T, ncol = p))
fixed_effects_by_T <- lapply(X, function(x) x %*% beta)

values <- data_frame(
  u = u,
  initial_y = u + sapply(fixed_effects_by_T,function(x) x[1]) + sigma_e / 
    sqrt(1 - delta^2) * rnorm(I), # p(y_1 | u, sigma_e, delta),
  mu = fixed_effects_by_T
)

dfs <- list()
for (i in 1:I) {
  out <- data.frame(y = rep(NA, T),
                    time = rep(NA, T),
                    u = values$u[i],
                    individual = i)
  for(t in 1:T) {
    out$time[t] <- t
    if(t == 1) {
      out$y[t] <- values$initial_y[i]
    } else {
      out$y[t] <- rnorm(1, delta * (out$y[t-1] - values$u[i] - values$mu[[i]][t]) + values$u[i] + values$mu[[i]][t], sigma_e)
    }
  }
  out <- out %>% mutate(
    lagged_y = lag(y)
  )
  dfs[[i]] <- out
}
simulated_panel <- do.call(rbind,dfs)

y <- matrix(0, T, I)
y_lag <- matrix(0, T, I)
for (i in 1:I) {
  y[,i] <- filter(simulated_panel, individual == i)$y
  y_lag[2:T,i] <- filter(simulated_panel, individual == i)$lagged_y[2:T]
}

mod <- stan_model('dynamic-panel_w_predictors.stan')
stan_data <- list(I = I, T = T, p = p, y = y, y_lag = y_lag, X = X)

fit <- sampling(mod, data = stan_data, iter = 2000, cores = 4, chains = 4, control = list(adapt_delta = 0.9))
library(bayesplot)
mcmc_recover_intervals(as.matrix(fit, pars = 'beta'), true = beta)
mcmc_recover_intervals(as.matrix(fit, pars = 'delta'), true = delta)
mcmc_recover_intervals(as.matrix(fit, pars = 'sigma_u'), true = sigma_u)
mcmc_recover_intervals(as.matrix(fit, pars = 'sigma_e'), true = sigma_e)

```

```
data {
  int<lower=1> I;
  int<lower=1> T;
  int<lower=1> p;
  matrix[T,I] y;
  matrix[T,I] y_lag;
  matrix[T,p] X[I];
}
parameters {
  real<lower=0,upper=1> delta_raw;
  real<lower=0> sigma_u;
  real<lower=0> sigma_e;
  vector[I] u_raw;
  
  // Fixed effects
  vector[p] beta;
}
transformed parameters {
  vector[I] u = sigma_u * u_raw;
  real delta = delta_raw * 2 - 1;
  matrix[T, I] mu;
  for (i in 1:I)
    mu[,i] = u[i] + X[i] * beta;
}
model {
  real one_minus_delta_sq = (1 - delta) * (1 + delta);
 
  for (i in 1:I) {
    for(j in 1:T){
      if(j == 1){
        y[j,i] ~ normal(mu[1,i], sigma_e / sqrt(one_minus_delta_sq));
      }
      else{
        y[j,i] ~ normal(mu[j,i] * (1 - delta) + delta * y_lag[j,i], sigma_e);
      }
    }
  }
  
  u_raw ~ normal(0, 1);
  sigma_u ~ normal(0, 2);
  sigma_e ~ normal(0, 1);
  delta_raw ~ beta(4, 4);
  beta ~ normal(0, 5);
}
generated quantities{
  matrix[T,I] y_rep;
  real one_minus_delta_sq = (1 - delta) * (1 + delta);
  
  for (i in 1:I) {
    for(j in 1:T){
      if(j == 1){
        y_rep[j,i] = normal_rng(mu[1,i], sigma_e / sqrt(one_minus_delta_sq));
      }
      else{
        y_rep[j,i] = normal_rng(mu[j,i] * (1 - delta) + delta * y_lag[j,i], sigma_e);
      }
    }
  }
}

```

The model seems to do OK at pulling out the parameters, though there are divergent transitions if you don’t dial up `adapt_delta` to 0.9. This example should be extensible to other models with group-level predictors.

---

<div class="post-metadata">

**Author:** ![c5v](https://avatars.discourse-cdn.com/v4/letter/c/e19adc/32.png) [@c5v](https://discourse.mc-stan.org/u/c5v)\
**Post date:** [February 4, 2019, 7:23pm UTC](https://discourse.mc-stan.org/t/dynamic-panel-data-models-with-stan/5136/32 "2019-02-04T19:23:31Z")

</div>

@rtrangucci Thanks you so much! I am going to make this work in my multilevel context. I will let you know if it works!

sidenote: `Ni` is just the total number of observations

---

<div class="post-metadata">

**Author:** ![c5v](https://avatars.discourse-cdn.com/v4/letter/c/e19adc/32.png) [@c5v](https://discourse.mc-stan.org/u/c5v)\
**Post date:** [February 11, 2019, 6:46pm UTC](https://discourse.mc-stan.org/t/dynamic-panel-data-models-with-stan/5136/33 "2019-02-11T18:46:42Z")

</div>

@rtrangucci I have one final question about your model. I want to make the lagged variable parameter \delta also dependent on i. So the model then becomes y\_{it} = (\alpha+u\_{i}) + x\_{i}\beta + \delta\_{i} + \varepsilon\_{it}. In my multilevel setting I could again estimate \delta\_i using non-centered parameterization and use:

```stan
parameters {
  ...
  real delta_raw;
  real u_delta_raw[Nk];
  real<lower=0> sigma_u_delta
  ...
}
transformed parameters{
  real u_delta[N];
  for(i in 1:Nk){//number of individuals i
  u_delta[i] = sigma_u_delta*u_delta_raw[i];
  }

  for(k in 1:Nk){//number of individuals i
   delta[k] = delta_raw + u_delta[k];
  }
}
model {
...
u_delta_raw ~ normal(0, 1);
sigma_u_delta~normal(0,1);
...
}

```

However, this does not necessarily mean that my \delta\_i is in the interval [-1,1]. Do you maybe know how to incorporate this?

Many thanks in advance!

---

<div class="post-metadata">

**Author:** ![rtrangucci](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/rtrangucci/32/8375_2.png) [@rtrangucci](https://discourse.mc-stan.org/u/rtrangucci)\
**Post date:** [February 14, 2019, 7:19pm UTC](https://discourse.mc-stan.org/t/dynamic-panel-data-models-with-stan/5136/34 "2019-02-14T19:19:23Z")

</div>

Hi @c5v,

This is almost all the way there. You can change this line:

```
  for(k in 1:Nk){//number of individuals i
   delta[k] = delta_raw + u_delta[k];
  }

```

to

```
  for(k in 1:Nk){//number of individuals i
   delta[k] = inv_logit(delta_raw + u_delta[k]) * 2 - 1;
  }

```

which will get you back to the [-1,1] scale.

---

<div class="post-metadata">

**Author:** ![c5v](https://avatars.discourse-cdn.com/v4/letter/c/e19adc/32.png) [@c5v](https://discourse.mc-stan.org/u/c5v)\
**Post date:** [February 26, 2019, 7:50pm UTC](https://discourse.mc-stan.org/t/dynamic-panel-data-models-with-stan/5136/35 "2019-02-26T19:50:27Z")

</div>

@rtrangucci first of all thanks for all the help!

As final (really) question I would like to ask whether you know a paper or some formal description of the model:

> [@rtrangucci](#):
>
> Y\_{i,t}=μ\_{i}+δ(Y\_{i,t−1}−μ\_{i})+ϵ\_{i,t}, ϵ\_{i,t}∼Normal(0,1)

I would like to use it in a paper for university but I cannot find any paper which uses this type of model which I can cite.

This one is not really what you are doing

> [@Saram\_Han](#):
>
> Just in case if someone is interested in this approach who has no econometric background, I found the following paper useful to understand rtrangucci’s code.  
> [https://www.researchgate.net/profile/Hashem\_Pesaran/publication/226673422\_Random\_Coefficient\_Models/links/0fcfd5089366b8067f000000/Random-Coefficient-Models.pdf](https://www.researchgate.net/profile/Hashem_Pesaran/publication/226673422_Random_Coefficient_Models/links/0fcfd5089366b8067f000000/Random-Coefficient-Models.pdf)

---

<div class="post-metadata">

**Author:** ![aornugent](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/aornugent/32/1142_2.png) [@aornugent](https://discourse.mc-stan.org/u/aornugent)\
**Post date:** [March 3, 2019, 9:10pm UTC](https://discourse.mc-stan.org/t/dynamic-panel-data-models-with-stan/5136/36 "2019-03-03T21:10:51Z")

</div>

I’m digging this approach:

> [@rtrangucci](#):
>
> \begin{align} Y\_{i,t} &=μ\_i+δ(Y\_{i,t−1}−μ\_i)+ϵ\_{i,t} \\ ϵ\_{i,t} & \sim Normal(0,1) \end{align}

@rtrangucci - do you know of any smart ways to deal with irregularly spaced observations?

Imputing missing data gets me most of the way there, but my estimates of \sigma\_u and \sigma\_e become non-identifiable as the proportion of missing values increases.

Appreciated,  
Andrew

---

<div class="post-metadata">

**Author:** ![c5v](https://avatars.discourse-cdn.com/v4/letter/c/e19adc/32.png) [@c5v](https://discourse.mc-stan.org/u/c5v)\
**Post date:** [March 26, 2019, 1:37pm UTC](https://discourse.mc-stan.org/t/dynamic-panel-data-models-with-stan/5136/37 "2019-03-26T13:37:27Z")

</div>

Just to make sure, this approach works because now the correlation between \mu\_{i} and \delta(y\_{i,t-1}-\mu\_i) equals 0, right? Since E \left[\mu\_i (\delta(y\_{i,t-1} -\mu\_{i})) \right] = \mu\_{i}\delta (E\left[y\_{i,t-1} \right] - \mu\_i) = 0, using E\left[y\_{i,t-1} \right] = \mu\_i

---

<div class="post-metadata">

**Author:** ![joeHoover](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/joehoover/32/5477_2.png) [@joeHoover](https://discourse.mc-stan.org/u/joeHoover)\
**Post date:** [August 17, 2019, 8:33am UTC](https://discourse.mc-stan.org/t/dynamic-panel-data-models-with-stan/5136/39 "2019-08-17T08:33:27Z")

</div>

Hi @rtrangucci,

Like others, I’ve found this code really useful. So, thank you!

However, I’m having trouble estimating group-specific \delta\_i on some simulated data.

I’ve reduced my model to pretty much exactly what you suggested for the random intercept model with group-specific \delta\_i. However, the estimated \delta\_i seem to be stuck around .08-.09, even when the true mean is .4.

Interestingly, when I remove the random intercept from the model, the range restriction on \delta\_i goes away and the estimation is decent.

So, I’m wondering whether I might have implemented something wrong or whether you might have any insight into this behavior.

@c5v, it looks like you worked quite a bit with these models. Were you able to estimate random \delta\_i?

If anyone has any insight into this, I’d love to hear it!

I’ve attached the simulated [data](https://discourse.mc-stan.org/uploads/short-url/muSDHwb4zKUyD6Z1MGDY5Qtc4CM.csv) (231.7 KB) I’m working with as well as [the code](https://discourse.mc-stan.org/uploads/short-url/1xumxM4XEJXwGJMwbp8WYoqLi7r.R) (798 Bytes) for running and checking the model.

The stan model I’m using is:

```
data {
  int<lower=1> S; // Number of states
  int<lower=1> T; // Number of time points
  int<lower=1> P; // Number of predictors
  matrix[T,S] y; // Outcome variable
  matrix[T,S] y_lag; // Lagged outcome
  matrix[T,P] X[S]; // List of [T, P] matrices containing intercept and predictors
}
parameters {

  real<lower=0,upper=1> delta_raw; // Raw AR(1) coefficient
  real u_delta_raw[S];
  real<lower=0> sigma_u_delta;
  
  
  vector[S] u_raw;
  real<lower=0> sigma_u;
  
  real<lower=0> sigma_e; // Residual error

  vector[P] beta;
}
transformed parameters {
  
  real delta[S];
  real u_delta[S];
  vector[S] u;
  matrix[T, S] mu;
  
  u = sigma_u * u_raw;

  for(s in 1:S){//number of states s
    u_delta[s] = sigma_u_delta*u_delta_raw[s];
  }
  
  for(s in 1:S){//number of states s
   delta[s] = inv_logit(delta_raw + u_delta[s]) * 2 - 1;
  }
  
  for (s in 1:S)
    //mu[,s] = u[s] + X[s] * beta;
    mu[,s] = u[s] + X[s] * beta;
    
}
model {
  

  for (s in 1:S) {
    for(j in 1:T){
      if(j == 1){
        y[j,s] ~ normal(mu[1,s], sigma_e / sqrt(1 - delta[s]) * (1 + delta[s]));
      }
      else{
        y[j,s] ~ normal(mu[j,s] * (1 - delta[s]) + delta[s] * y_lag[j,s], sigma_e);
      }
    }
  }
  
  sigma_u ~ normal(0,2);
  u_raw ~ normal(0,1);
  sigma_e ~ normal(0, 1);
  beta[1] ~ normal(0, 2);
  
  delta_raw ~ beta(6, 6);
  u_delta_raw ~ normal(0, 1);
  sigma_u_delta ~ normal(0, 2);

}
generated quantities{
  matrix[T,S] y_rep;

  for (s in 1:S) {
    for(j in 1:T){
      if(j == 1){
        y_rep[j,s] = normal_rng(mu[1,s], sigma_e / sqrt(1 - delta[s]) * (1 + delta[s]));
      }
      else{
        y_rep[j,s] = normal_rng(mu[j,s] * (1 - delta[s]) + delta[s] * y_lag[j,s], sigma_e);
      }
    }
  }
}

```

---

<div class="post-metadata">

**Author:** ![ignacio](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ignacio/32/1015_2.png) [@ignacio](https://discourse.mc-stan.org/u/ignacio)\
**Post date:** [November 15, 2019, 7:25pm UTC](https://discourse.mc-stan.org/t/dynamic-panel-data-models-with-stan/5136/40 "2019-11-15T19:25:56Z")

</div>

@rtrangucci thanks a lot for your post. I’m having troubles understanding

> [@rtrangucci](#):
>
> ```stan
> 
> for (i in 1:I) { 
> target += normal_lpdf(y[1,i] | u[i], sigma_e / sqrt(one_minus_delta_sq)); 
> }
> 
> ```

How do you get that target for t=1 from this equation?

Y\_{i,t} = \mu\_i + \delta (Y\_{i,t-1} - \mu\_i) + \epsilon\_{i,t}

I understand the target for `2:T` but the first one is confusing me.

---

<div class="post-metadata">

**Author:** ![Max\_Mantei](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/max_mantei/32/7912_2.png) [@Max\_Mantei](https://discourse.mc-stan.org/u/Max_Mantei)\
**Post date:** [November 15, 2019, 9:57pm UTC](https://discourse.mc-stan.org/t/dynamic-panel-data-models-with-stan/5136/41 "2019-11-15T21:57:12Z")

</div>

Hi @ignacio!

I’m not the great @rtrangucci, but let me give it a try.

So, unfortunately the notation changed a bit over the course of this thread—but, basically what’s `u[i]` in the code is \mu\_i in equation, right? So far so good.

In the answers above you can find that \mathbb{E}[Y\_{i,t}]=\mu\_i, so in the code that is the `u[i]`—the mean of `y[1,i]`.

For the variance we have

\begin{align} \text{Var}[Y\_{i,t}]&=\text{Var}[\delta Y\_{i,t-1}]+\text{Var}[\epsilon\_{i,t}] \\ \text{Var}[Y\_{i,t}]&=\delta^2\text{Var}[Y\_{i,t-1}] +\sigma^2\_\epsilon\\ \text{Var}[Y\_{i,t}] - \delta^2\text{Var}[Y\_{i,t-1}]&=\sigma^2\_\epsilon. \end{align}

Now, I think we need to assume \text{Var}[Y\_{i,t}] = \text{Var}[Y\_{i,t-1}], which is reasonable (assume iid residuals / a stationary process). Then,

\begin{align} \text{Var}[Y\_{i,t}] - \delta^2\text{Var}[Y\_{i,t}]&=\sigma^2\_\epsilon\\ \text{Var}[Y\_{i,t}](1 - \delta^2)&=\sigma^2\_\epsilon\\ \text{Var}[Y\_{i,t}]&=\frac{\sigma^2\_\epsilon}{(1 - \delta^2)}\\ \text{Sd}[Y\_{i,t}]&=\frac{\sigma\_\epsilon}{\sqrt{1 - \delta^2}}, \end{align}

which is the `sigma_e / sqrt(one_minus_delta_sq)` part of the code.

I learned this at the [Helsinki StanCon Tutorial with Jonah](https://github.com/jgabry/stancon2018helsinki_intro/blob/master/Pest_Control_Example.Rmd). Have a look—he also discusses a GP formulation of this AR(1) process.

Cheers! :)

---

<div class="post-metadata">

**Author:** ![ignacio](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ignacio/32/1015_2.png) [@ignacio](https://discourse.mc-stan.org/u/ignacio)\
**Post date:** [November 15, 2019, 10:23pm UTC](https://discourse.mc-stan.org/t/dynamic-panel-data-models-with-stan/5136/42 "2019-11-15T22:23:25Z")

</div>

Thanks a lot @Max_Mantei! Do you know by any chance if there is a video for @jonah’s tutorial?

[Previous page](https://discourse.mc-stan.org/t/dynamic-panel-data-models-with-stan/5136.md?page=1)

[Next page](https://discourse.mc-stan.org/t/dynamic-panel-data-models-with-stan/5136.md?page=3)