#  ex 9.3  HIPS Micro-simulation and Forward Sampling
#MC estimates for single stratum, allowing for parameter uncertainty in revision hazard, h

# n = 10000 updates (1 per simulated set of parameter values) are required for this model; monitor C, BL, BQ

model {

   # Cost and benefit equations in closed form:
   ####################################

   # Costs
   for(t in 1:N) {
     ct[t] <- inprod( pi[t, ] , c[ ] ) / pow( (1+delta.c) , (t-1) )     # discounted at rate delta.c
   } 
   C <- C0 + sum(ct[])

   # Benefits - life expectancy
   for(t in 1:N) {
     blt[t] <- inprod( pi[t, ] , bl[ ] ) / pow( (1+delta.b) , (t-1) )
   } 
   BL <- sum(blt[])

   # Benefits - QALYs
   for(t in 1:N) {
     bqt[t] <- inprod(pi[t,], bq[])/pow((1+delta.b), (t-1))
   } 
   BQ <- sum(bqt[])

   # Markov model probabilities:
   #######################

   # Transition matrix
   for(t in 1:N) {
     Lambda[t,1,1] <- 1 -  gamma[t] - lambda[t]
     Lambda[t,1,2] <- gamma[t] * lambda.op
     Lambda[t,1,3] <- gamma[t] *(1-lambda.op)
     Lambda[t,1,4] <- 0
     Lambda[t,1,5] <- lambda[t] 

     Lambda[t,2,1] <- 0
     Lambda[t,2,2] <- 0 
     Lambda[t,2,3] <- 0 
     Lambda[t,2,4] <- 0 
     Lambda[t,2,5] <- 1 

     Lambda[t,3,1] <- 0
     Lambda[t,3,2] <- 0 
     Lambda[t,3,3] <- 0
     Lambda[t,3,4] <- 1 -  lambda[t]
     Lambda[t,3,5] <- lambda[t]

     Lambda[t,4,1] <- 0
     Lambda[t,4,2] <- rho * lambda.op
     Lambda[t,4,3] <- rho * (1- lambda.op)
     Lambda[t,4,4] <- 1 - rho - lambda[t]
     Lambda[t,4,5] <- lambda[t]

     Lambda[t,5,1] <- 0
     Lambda[t,5,2] <- 0 
     Lambda[t,5,3] <- 0
     Lambda[t,5,4] <- 0
     Lambda[t,5,5] <- 1

     gamma[t] <- h * (t-1)
  }

 
  # Marginal probability of being in each state at time 1
    pi[1,1] <- 1-lambda.op ;  pi[1,2]<-0  ;    pi[1,3] <- 0 ;  pi[1,4]<-0;  pi[1,5] <- lambda.op

  # Marginal probability of being in each state at time t>1
    for(s in 1:S) {
       for(t in 2:N) {
       pi[t,s] <- inprod(pi[(t-1),], Lambda[t,,s])
     }
   }
# age-sex specific revision hazard

    h      <- exp(logh)
    logh  ~ dnorm(logh0, tau)
    logh0 <- log(h0)  
    
 }