# varying-intercept, varying-slope radon model for earnings by ethnicity and age
# modeling the correlation

model {
  for (i in 1:n){
    y[i] ~ dnorm (y.hat[i], tau.y)
    y.hat[i] <- a[eth[i],age[i]] + b[eth[i],age[i]]*x[i]
  }
  tau.y <- pow(sigma.y, -2)
  sigma.y ~ dunif (0, 100)

  for (j in 1:n.eth){
    for (k in 1:n.age){
      a[j,k] <- B[j,k,1]
      b[j,k] <- B[j,k,2]
      B[j,k,1:2] ~ dmnorm (B.hat[j,k,], Tau.B[,])
      for (l in 1:2){
        B.hat[j,k,l] <- mu[l] + G.eth[j,l] + G.age[k,l]
      }
    }
  }
  mu[1:2] ~ dmnorm (zeros[], Tau.0[,])
  for (j in 1:n.eth){
    G.eth[j,1:2] ~ dmnorm (zeros[], Tau.eth[,])
  }
  for (k in 1:n.age){
    G.age[k,1:2] ~ dmnorm (zeros[], Tau.age[,])
  }

  Tau.B[1:2,1:2] <- inverse(Sigma.B[,])
  for (l in 1:2){
    Sigma.B[l,l] <- pow(sigma.B[l], 2)
    sigma.B[l] ~ dunif (0, 100)
  }
  Sigma.B[1,2] <- rho.B*sigma.B[1]*sigma.B[2]
  Sigma.B[2,1] <- Sigma.B[1,2]
  rho.B ~ dunif (-1, 1)

  Tau.eth[1:2,1:2] <- inverse(Sigma.eth[,])
  for (l in 1:2){
    Sigma.eth[l,l] <- pow(sigma.eth[l], 2)
    sigma.eth[l] ~ dunif (0, 100)
  }
  Sigma.eth[1,2] <- rho.eth*sigma.eth[1]*sigma.eth[2]
  Sigma.eth[2,1] <- Sigma.eth[1,2]
  rho.eth ~ dunif (-1, 1)

  Tau.age[1:2,1:2] <- inverse(Sigma.age[,])
  for (l in 1:2){
    Sigma.age[l,l] <- pow(sigma.age[l], 2)
    sigma.age[l] ~ dunif (0, 100)
  }
  Sigma.age[1,2] <- rho.age*sigma.age[1]*sigma.age[2]
  Sigma.age[2,1] <- Sigma.age[1,2]
  rho.age ~ dunif (-1, 1)
}