model{
  for (i in 1:p){
    for (j in 1:A){
      for (k in 1:T){
        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(alpha[j]+kappa[i,k]*beta[j]+gamma[k-j+A])
      }
    }
  }


  for (j in 2:(C-2)){
  gamma[j]<-gamma_rest[j-1]
  }

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

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

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

    kappa[i,1]<--sum(kappa_rest_mat[i,])

  }

  for (i in 1:p){
    for (j in 2:(T-1)){
    prior_mean_kappa[i,j]<-eta_mean[j+1]+rho*(kappa_rest_mat[i,j-1]-eta_mean[j])	
    kappa_rest_mat[i,j]~dmnorm(prior_mean_kappa[i,j],i_sigma2_kappa)
    }
    prior_mean_kappa[i,1]<-eta_mean[2]+0 # the 0 should have been rho*(kappa[i,1]-eta_mean[1]) 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)
  }

  eta_mean<-c(matrix_kappa_X%*%eta)

  for (j in 2:A){
  beta[j]<-beta_rest[j-1]
  }

  beta[1]<-1-sum(beta_rest)


    beta_rest[1:(A-1)]~dmnorm(prior_mean_beta,prior_prec_beta)

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

  eta~dmnorm(prior_mean_eta,prior_prec_eta) #supply prior_prec_eta<-solve(matrix(c(400,0,0,2),nrow=2));prior_mean_eta<-c(0,0)
  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)

  phi~dgamma(25,0.05)

}
