# Spatial modeling via joint specification

# This function implements a two-stage hierarchical model
# using a second-stage ML weighted least squares
# betas are the first-stage estimates
# se are the first-stage standard errors
# V is the first-stage covariance matrix
# Z is the second-stage design matrix
# SIGMA is the second-stage covariance matrix
# pi is vector of coefficients for the second-stage

hld.joint.Z.SIGMA <- function(Z, SIGMA, tau, betas, se) {
  n <- length(betas)
  V <- diag(se^2)
  T <- (tau^2)*SIGMA
  W <- solve(V + T)
  pi <-  solve(t(Z) %*% W %*% Z) %*% t(Z) %*% W %*% betas
  second.betas <- Z %*% pi
  I <- diag(1, nrow=(n))
  B <- W %*% V
  hld.betas <- B %*% second.betas + (I - B) %*% betas
  # Interval Estimation
  H <- Z %*% solve( t(Z) %*% W %*% Z ) %*% t(Z) %*% W
  C <- V %*% (I - (t(I - H) %*% B) )
}