model{
  for (i in 1:p){
    for (j in 1:A){
      for (k in 1:(T+h)){
        dxt[i,j,k]~dnegbin(prob[i,j,k],phi)
        prob[i,j,k]<-phi/(phi+q[i,j,k]*ext[i,j,k])
        q[i,j,k]<-exp(kappa1[i,k]+kappa2[i,k]*(x[j]-xbar)+gamma[i,k-j+A])
      }
    }
  }

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

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

  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-2)){
    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-2)]*gamma_rest_mat[i,]))/(1-C)	
    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_rev[3:C]*gamma_rest_mat[i,]))/(1-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 1:T){
  kappa1[i,j]<-kappa1_rest_mat[i,j]
  }

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

  }

  for (i in 1:p){
    for (j in 2:T){
    prior_mean_kappa1[i,j]<-eta1_mean[j]+rho1*(kappa1_rest_mat[i,j-1]-eta1_mean[j-1])	
    kappa1_rest_mat[i,j]~dnorm(prior_mean_kappa1[i,j],i_sigma2_kappa1)
    }
    prior_mean_kappa1[i,1]<-eta1_mean[1]+0 # the 0 should have been rho1*(kappa1[i,1]-eta1_mean[1]) but it does not work because the model could not be uniquely defined	
    kappa1_rest_mat[i,1]~dnorm(prior_mean_kappa1[i,1],i_sigma2_kappa1)
  }

  for (i in 1:p){
    for (j in 2:h){
    prior_mean_kappa1_forecast[i,j]<-eta1_mean[T+j]+rho1*(kappa1_forecast_rest_mat[i,j-1]-eta1_mean[T+j-1])	
    kappa1_forecast_rest_mat[i,j]~dnorm(prior_mean_kappa1_forecast[i,j],i_sigma2_kappa1)
    }
    prior_mean_kappa1_forecast[i,1]<-eta1_mean[T+1]+rho1*(kappa1_rest_mat[i,T-1]-eta1_mean[T])	
    kappa1_forecast_rest_mat[i,1]~dnorm(prior_mean_kappa1_forecast[i,1],i_sigma2_kappa1)
  }

  eta1_mean<-c(matrix_kappa_X%*%eta1)


  for (i in 1:p){
  for (j in 1:T){
  kappa2[i,j]<-kappa2_rest_mat[i,j]
  }

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

  }

  for (i in 1:p){
    for (j in 2:T){
    prior_mean_kappa2[i,j]<-eta2_mean[j]+rho2*(kappa2_rest_mat[i,j-1]-eta2_mean[j-1])	
    kappa2_rest_mat[i,j]~dnorm(prior_mean_kappa2[i,j],i_sigma2_kappa2)
    }
    prior_mean_kappa2[i,1]<-eta2_mean[1]+0 # the 0 should have been rho2*(kappa2[i,1]-eta2_mean[1]) but it does not work because the model could not be uniquely defined	
    kappa2_rest_mat[i,1]~dnorm(prior_mean_kappa2[i,1],i_sigma2_kappa2)
  }

  for (i in 1:p){
    for (j in 2:h){
    prior_mean_kappa2_forecast[i,j]<-eta2_mean[T+j]+rho2*(kappa2_forecast_rest_mat[i,j-1]-eta2_mean[T+j-1])	
    kappa2_forecast_rest_mat[i,j]~dnorm(prior_mean_kappa2_forecast[i,j],i_sigma2_kappa2)
    }
    prior_mean_kappa2_forecast[i,1]<-eta2_mean[T+1]+rho2*(kappa2_rest_mat[i,T-1]-eta2_mean[T])	
    kappa2_forecast_rest_mat[i,1]~dnorm(prior_mean_kappa2_forecast[i,1],i_sigma2_kappa2)
  }

  eta2_mean<-c(matrix_kappa_X%*%eta2)

  eta1~dmnorm(prior_mean_eta1,prior_prec_eta1) #supply prior_prec_eta1<-solve(matrix(c(400,0,0,2),nrow=2));prior_mean_eta1<-c(0,0)
  rho1~dnorm(0,1) T(-1,1)
  sigma2_kappa1<-1/i_sigma2_kappa1
  i_sigma2_kappa1~dgamma(1,0.0001)

  eta2~dmnorm(prior_mean_eta2,prior_prec_eta2) #supply prior_prec_eta2<-solve(matrix(c(400,0,0,2),nrow=2));prior_mean_eta2<-c(0,0)
  rho2~dnorm(0,1) T(-1,1)
  sigma2_kappa2<-1/i_sigma2_kappa2
  i_sigma2_kappa2~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)

  phi~dgamma(25,0.05)

}
