model{
  for (i in 1:p){
    for (j in 1:A){
      for (k in 1:(T+h)){
        dxt[i,j,k]~dpois(q[i,j,k]*ext[i,j,k])
        q[i,j,k]<-exp(alpha[i,j]+kappa[i,k]+beta[i,j]*(t_forecast[k]-tbar)+gamma[i,k-j+A])
      }
    }
  }

  for (i in 1:p){
  for (j in 2:(C-2)){
  gamma[i,j]<-gamma_rest_mat[i,j-1]
  }

  gamma[i,C]<-(sum(cohorts[1:(C-3)]*cohorts_rev[4:C]*gamma_rest_mat[i,]))/(C-1)
  gamma[i,C-1]<-(sum(cohorts[1:(C-3)]*cohorts_rev[3:(C-1)]*gamma_rest_mat[i,]))/(2-C)
  gamma[i,1]<-(sum(cohorts_rev[3:(C-1)]*cohorts_rev[4:C]*gamma_rest_mat[i,]))/((1-C)*(C-2))


  for (k in 1:h){
    gamma[i,k+C]<-gamma_forecast_rest_mat[i,k]
  }

  }

  for (i in 1:p){
    for (j in 3:(C-3)){
    prior_mean_gamma[i,j]<-gamma_rest_mat[i,j-1]+rho_gamma*(gamma_rest_mat[i,j-1]-gamma_rest_mat[i,j-2])	
    gamma_rest_mat[i,j]~dmnorm(prior_mean_gamma[i,j],i_sigma2_gamma)
    }
    prior_mean_gamma[i,2]<-gamma_rest_mat[i,1]	
    gamma_rest_mat[i,2]~dmnorm(prior_mean_gamma[i,2],i_sigma2_gamma*(1-rho_gamma^2))
    prior_mean_gamma[i,1]<-0	
    gamma_rest_mat[i,1]~dmnorm(prior_mean_gamma[i,1],0.01*i_sigma2_gamma)
  }

  for (i in 1:p){
    for (j in 3:h){
    prior_mean_gamma_forecast[i,j]<-gamma_forecast_rest_mat[i,j-1]+rho_gamma*(gamma_forecast_rest_mat[i,j-1]-gamma_forecast_rest_mat[i,j-2])	
    gamma_forecast_rest_mat[i,j]~dmnorm(prior_mean_gamma_forecast[i,j],i_sigma2_gamma)
    }
    prior_mean_gamma_forecast[i,2]<-gamma_forecast_rest_mat[i,1]+rho_gamma*(gamma_forecast_rest_mat[i,1]-gamma[i,C])  #or swap gamma[i,C] with (sum(cohorts[1:(C-3)]*cohorts_rev[4:C]*gamma_rest_mat[i,]))/(C-1)	
    gamma_forecast_rest_mat[i,2]~dmnorm(prior_mean_gamma_forecast[i,2],i_sigma2_gamma)
    prior_mean_gamma_forecast[i,1]<-gamma[i,C]+rho_gamma*(gamma[i,C]-gamma[i,C-1])  #or swap gamma[i,C-1] with (sum(cohorts[1:(C-3)]*cohorts_rev[3:(C-1)]*gamma_rest_mat[i,]))/(2-C)	
    gamma_forecast_rest_mat[i,1]~dmnorm(prior_mean_gamma_forecast[i,1],i_sigma2_gamma)
  }

  for (i in 1:p){
  for (j in 2:(T-1)){
  kappa[i,j]<-kappa_rest_mat[i,j-1]
  }

  kappa[i,1]<-(sum(t_rev[3:T]*kappa_rest_mat[i,]))/(1-T)
  kappa[i,T]<-(sum((t[1:(T-2)])*kappa_rest_mat[i,]))/(1-T)

  for (k in 1:h){
    kappa[i,k+T]<-kappa_forecast_rest_mat[i,k]
  }

  }

  for (i in 1:p){
    for (j in 2:(T-2)){
    prior_mean_kappa[i,j]<-rho*(kappa_rest_mat[i,j-1])	
    kappa_rest_mat[i,j]~dmnorm(prior_mean_kappa[i,j],i_sigma2_kappa)
    }
    prior_mean_kappa[i,1]<-rho*(0) # the 0 should have been (sum(t_rev[3:T]*kappa_rest_mat[i,]))/(1-T) but it does not work because the model could not be uniquely defined	
    kappa_rest_mat[i,1]~dmnorm(prior_mean_kappa[i,1],i_sigma2_kappa)
  }

  for (i in 1:p){
    for (j in 2:h){
    prior_mean_kappa_forecast[i,j]<-rho*(kappa_forecast_rest_mat[i,j-1])	
    kappa_forecast_rest_mat[i,j]~dmnorm(prior_mean_kappa_forecast[i,j],i_sigma2_kappa)
    }
    prior_mean_kappa_forecast[i,1]<-rho*(kappa[i,T])	 #or swap kappa[i,T] with (sum((t[1:(T-2)])*kappa_rest_mat[i,]))/(1-T)
    kappa_forecast_rest_mat[i,1]~dmnorm(prior_mean_kappa_forecast[i,1],i_sigma2_kappa)
  }

  for (i in 1:p){

    beta[i,1:A]~dmnorm(prior_mean_beta,prior_prec_beta)

  }

  for (i in 1:p){
  for (j in 1:A){
  alpha[i,j]~dnorm(0,0.001)
  }
  }

  rho~dnorm(0,1) T(-1,1)
  sigma2_kappa<-1/i_sigma2_kappa
  i_sigma2_kappa~dgamma(1,0.0001)
  rho_gamma~dnorm(0,1) T(-1,1)
  i_sigma2_gamma<-1/sigma2_gamma
  sigma2_gamma<-sigma_gamma^2
  sigma_gamma~dunif(0,1)

}
