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

  for (j in 1:A){
      for (k in 1:(T+h)){
        dxt[p,j,k]~dpois(q[p,j,k]*ext[p,j,k])
        q[p,j,k]<-exp(alpha[p,j]+Beta[j]*Kappa[k])
      }
    }


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

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

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

  }

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

  for (i in 1:(p-1)){
    for (j in 2:h){
    prior_mean_kappa_forecast[i,j]<-eta_mean_kappa[T+j]+rho_kappa*(kappa_forecast_rest_mat[i,j-1]-eta_mean_kappa[T+j-1])	
    kappa_forecast_rest_mat[i,j]~dmnorm(prior_mean_kappa_forecast[i,j],i_sigma2_kappa)
    }
    prior_mean_kappa_forecast[i,1]<-eta_mean_kappa[T+1]+rho_kappa*(kappa_rest_mat[i,T-1]-eta_mean_kappa[T])	
    kappa_forecast_rest_mat[i,1]~dmnorm(prior_mean_kappa_forecast[i,1],i_sigma2_kappa)
  }

  eta_mean_kappa<-c(matrix_kappa_X%*%eta_kappa)

  #Kappa<-c(-sum(Kappa_rest),Kappa_rest)
  #Kappa_rest~dmnorm(prior_mean_Kappa,prior_prec_Kappa)

  for (j in 2:T){
  Kappa[j]<-Kappa_rest[j-1]
  }

  Kappa[1]<--sum(Kappa_rest)

  for (k in 1:h){
    Kappa[k+T]<-Kappa_forecast_rest[k]
  }

    for (j in 2:(T-1)){
    prior_mean_Kappa[j]<-eta_mean_Kappa[j+1]+rho_Kappa*(Kappa_rest[j-1]-eta_mean_Kappa[j])	
    Kappa_rest[j]~dmnorm(prior_mean_Kappa[j],i_sigma2_Kappa)
    }
    prior_mean_Kappa[1]<-eta_mean_Kappa[2]+0 # the 0 should have been rho*(Kappa[1]-eta_mean_Kappa[1]) but it does not work because the model could not be uniquely defined	
    Kappa_rest[1]~dmnorm(prior_mean_Kappa[1],i_sigma2_Kappa)

    for (j in 2:h){
    prior_mean_Kappa_forecast[j]<-eta_mean_Kappa[T+j]+rho_Kappa*(Kappa_forecast_rest[j-1]-eta_mean_Kappa[T+j-1])	
    Kappa_forecast_rest[j]~dmnorm(prior_mean_Kappa_forecast[j],i_sigma2_Kappa)
    }
    prior_mean_Kappa_forecast[1]<-eta_mean_Kappa[T+1]+rho_Kappa*(Kappa_rest[T-1]-eta_mean_Kappa[T])	
    Kappa_forecast_rest[1]~dmnorm(prior_mean_Kappa_forecast[1],i_sigma2_Kappa)

  eta_mean_Kappa<-c(matrix_kappa_X%*%eta_Kappa)

  for (i in 1:(p-1)){
  for (j in 2:A){
  beta[i,j]<-beta_rest_mat[i,j-1]
  }

  beta[i,1]<-1-sum(beta_rest_mat[i,])

  }

  for (i in 1:(p-1)){

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

  }

  Beta<-c(1-sum(Beta_rest),Beta_rest)
  Beta_rest~dmnorm(prior_mean_Beta,prior_prec_Beta)
  
  for (i in 1:p){
  for (j in 1:A){
  alpha[i,j]~dnorm(0,0.001)
  }
  }

  eta_kappa~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_kappa~dnorm(0,1) T(-1,1)
  sigma2_kappa<-1/i_sigma2_kappa
  i_sigma2_kappa~dgamma(1,0.0001)

  eta_Kappa~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_Kappa~dnorm(0,1) T(-1,1)
  sigma2_Kappa<-1/i_sigma2_Kappa
  i_sigma2_Kappa~dgamma(1,0.0001)

  phi~dgamma(25,0.05)

}
