Covariates in the Covariance Gibbs

gistfile1.r

## Function to implment the MCMC sampler
polling.gibbs.kal <- function(y, size, dates, states, big, ## data
                              F, X=NULL, XL = NULL, Z=NULL,     ## covariates
                              m.0, C.0, b = NULL, B = NULL,  ## hyper-parameters
                              c.0, d.0, e.0, f.0, g.0, h.0,
                              starts, phi, sigma.a, nu, beta,          ## starting values
                              tune, A, tune.phi,
                              iters = 10000, thin = 1, burnin = 0) { ## mcmc stuff

  aa <- proc.time()
  ## pre-gibbs cleaning (calculate variances, setup times, etc)
  require(MCMCpack)
  require(mvtnorm)
  
  days <- dates- min(dates) + 1  
  T <- max(days)
  num.states <- length(unique(states))
  state.names <- unique(states)

  ## starting values
  ## gammas is the matrix (time by state) of state effects
  gammas <- matrix(starts, nrow = T, ncol = num.states)

  ## Create "projection" 
  S <- exp(-mh(Z, diag(phi, ncol(Z))))
  Delta <- sigma.a * S + diag(rep(nu, ncol(S)))

  ## Calculate number of draws
  draws <- (iters-burnin)/thin

  ## set up acceptance ratios
  acc.a <- acc.nu <- 0
  acc.phi <- rep(0, length(phi))
  
  ## holder matrices
  gamma.hold <- array(NA, dim=c(T,num.states,draws))
  phi.hold <- matrix(NA, nrow = ncol(Z), ncol = draws)
  sigma.a.hold <- rep(NA, times = draws)
  nu.hold <- rep(NA, times = draws)
  
  if (!is.null(X)) {
    beta.hold <- matrix(NA, nrow = ncol(XL), ncol = draws)
  }
  for (i in 1:iters) {
    
    ## update alpha.0, gamma.s0
    a <- matrix(0, nrow = T, ncol = num.states)
    m <- matrix(0, nrow = T, ncol = num.states)
    
    C <- array(0, dim = c(num.states,num.states,T))
    R <- array(0, dim = c(num.states,num.states,T))

    m[1,] <- m.0
    C[,,1] <- C.0
    
    ## Begin Forward Filtering
    for (t in 2:T) {

      ## y[,t] = F[[t]] %*% gammas[,t] + Sigma.t
      ## F relates the states (gammas) to the observations
      ## so we only need the rows that correspond to the
      ## states actually being used today
      a[t,] <- m[t-1,]
      R[,,t] <- C[,,t-1] + Delta
    
      ## slight modification for missing data.
      ## if there are no observations, then the posterior
      ## mean/variance is the prior mean/variance
      if (nrow(big[[t]]) > 0) {
        
        f <- F[[t]] %*% a[t,]
        if (!is.null(X)) {
          e <- big[[t]]$y - X[[t]] %*% beta - f
        } else {
          e <- big[[t]]$y - f
        }
        Sigma.t <- diag(big[[t]]$sigma.sq, nrow(big[[t]]))
                                      
        Qt <- F[[t]] %*% R[, ,t] %*% t(F[[t]]) + Sigma.t
        At <- R[,,t] %*% t(F[[t]]) %*% solve(Qt)

        m[t,] <- a[t,] + At %*% e
        C[,,t] <- R[,,t] - R[,,t] %*% t(F[[t]]) %*% solve(Qt) %*% F[[t]] %*% R[,,t]
      } else {
        m[t,] <- a[t,]
        C[,,t] <- R[,,t]
      }
    }
    
    ## Begin Backward Sampling
    gammas[T,] <- rmvnorm(1, m[T,], C[,,T])
    if (!is.null(X)) {
      beta.resid <- big[[T]]$y - F[[T]] %*% gammas[T,]
    }
      
    for (t in (T-1):1) {
      B.t <- C[,,t] %*% solve(R[,,t+1])
      h.t <- m[t,] + B.t %*% (gammas[t+1,] - a[t+1,])
      H.t <- C[,,t] - B.t %*% R[,,t+1] %*% t(B.t)

      gammas[t,] <- rmvnorm(1, h.t, H.t)
      if (!is.null(X)) {
        beta.resid <- c(beta.resid, big[[t]]$y - F[[t]] %*% gammas[t,])
      }
    }

    if (!is.null(X)) {
      ## update betas (polling house effects)
      Sigma.beta.star <- solve(t(XL) %*% diag(1/sigma.sq) %*% XL + B)
      beta.star <- Sigma.beta.star %*% (t(XL) %*% diag(1/sigma.sq) %*% beta.resid + B %*% b)

      beta <- mvrnorm(1, beta.star, Sigma.beta.star)

      ## normalize the house effects
      beta <- beta - mean(beta)
    }
    
    

    ## MH step for sigma.a
    sigma.a.can <- rlnorm(1, log(sigma.a), sd = tune[1,1])
    Delta.can <- sigma.a.can * S + diag(rep(nu, ncol(S)))
    hold.can <- -dlnorm(sigma.a.can, log(sigma.a), tune[1,1], log = TRUE)
    hold.cur <- -dlnorm(sigma.a, log(sigma.a.can), tune[1,1], log = TRUE)
    hold.can <- hold.can + log(dinvgamma(sigma.a.can, c.0, d.0))
    hold.cur <- hold.cur + log(dinvgamma(sigma.a, c.0, d.0))

    for (t in 2:T) {
      hold.can <- hold.can + dmvnorm(gammas[t,], gammas[t-1,], Delta.can, log = TRUE)
      hold.cur <- hold.cur + dmvnorm(gammas[t,], gammas[t-1,], Delta, log = TRUE)
    }


    r.a <- exp(hold.can - hold.cur)
    if (runif(1) <= r.a) {
      sigma.a <- sigma.a.can
      acc.a <- acc.a + 1
    }

    Delta <- sigma.a * S + diag(rep(nu, ncol(S)))

    ## MH step for nu ##
    nu.can <- rlnorm(1, log(nu), sd = tune[2,2])
    Delta.can <- sigma.a * S + diag(rep(nu.can, ncol(S)))
    hold.can <- -dlnorm(nu.can, log(nu), tune[2, 2], log = TRUE)
    hold.cur <- -dlnorm(nu, log(nu.can), tune[2, 2], log = TRUE)
    hold.can <- hold.can + log(dinvgamma(nu.can, g.0, h.0))
    hold.cur <- hold.cur + log(dinvgamma(nu, g.0, h.0))

    for (t in 2:T) {
      hold.can <- hold.can + dmvnorm(gammas[t,], gammas[t-1,], Delta.can, log = TRUE)
      hold.cur <- hold.cur + dmvnorm(gammas[t,], gammas[t-1,], Delta, log = TRUE)
    }

    r.a <- exp(hold.can - hold.cur)
    if (runif(1) <= r.a) {
      nu <- nu.can
      acc.nu <- acc.nu + 1
    }
    Delta <- sigma.a * S + diag(rep(nu, ncol(S)))

    ## MH for phi
    for (p in 1:length(phi)) {
      newphi <- rlnorm(1, log(phi[p]), sd = tune[p + 2, p + 2])
      phi.can <- phi
      phi.can[p] <- newphi
      S.can <- exp(-mh(Z, diag(c(phi.can))))
      Delta.can <- sigma.a * S.can + diag(rep(nu, ncol(S)))
    
      ## Jumping disttribution discount ##
      hold.can <- -dlnorm(newphi, log(phi[p]), tune[p + 2, p + 2], log = TRUE)
      hold.cur <- -dlnorm(phi[p], log(newphi), tune[p + 2, p + 2], log = TRUE)
    
      ## Priors
      hold.can <- hold.can + dhalfcauchy(newphi, A, log = TRUE)
      hold.cur <- hold.cur + dhalfcauchy(phi[p], A, log = TRUE)
    
      for (t in 2:T) {
        hold.can <- hold.can + dmvnorm(gammas[t,], gammas[t-1,], Delta.can, log = TRUE)
        hold.cur <- hold.cur + dmvnorm(gammas[t,], gammas[t-1,], Delta, log = TRUE)
      }
      
      ## Update phi[p] and Delta
      r.a <- exp(hold.can - hold.cur)
      if (runif(1) <= r.a) {
        phi <- phi.can
        acc.phi[p] <- acc.phi[p] + 1
        S <- S.can
      } 
      Delta <- sigma.a * S + diag(rep(nu, ncol(S)))    
    }

         
    acc <- c(acc.a, acc.nu, acc.phi)/i

    if ( (((i-burnin) %% thin)==0) & (i >burnin)) {
      gamma.hold[,,((i-burnin)/thin)] <- gammas
      phi.hold[,((i-burnin)/thin)] <- phi
      sigma.a.hold[((i-burnin)/thin)] <- sigma.a
      nu.hold[((i-burnin)/thin)] <- nu
      if (!is.null(X)) {
        beta.hold[,((i-burnin)/thin)] <- beta
        save(list = c("gamma.hold", "beta.hold", "phi.hold", "sigma.a.hold", "nu.hold", "acc"),
             file = "polls-out.RData")
      }
    }

    if ((i %% 10)==0) cat("\n")
    if ((i %% 100)==0) cat("Minutes Run: ",((proc.time()-aa)/60)[3],"\n")

    cat(i," ")
    #cat("\nCandidate: ",sqrt(sigma.a.can), sqrt(nu.can), phi.can, "\n")
    cat("New: ", sqrt(sigma.a), sqrt(nu), phi, "\n")
    cat("Acceptance rate: ", acc, "\n\n")
    flush.console()

  }
  if (!is.null(X)) {
    out <- list(gammas=gamma.hold, betas = beta.hold, phi = phi.hold,
            sigma.a = sigma.a.hold, nu = nu.hold,
            acc = acc)
  } else {
    out <- list(gammas=gamma.hold, phi = phi.hold,
            sigma.a = sigma.a.hold, nu = nu.hold,
            acc = acc)
  }
  return(out)
}

## Mahalanobis distance matrix
mh <- function(x, cov) {
  S <- matrix(NA, nrow(x), nrow(x))
  for (i in 1:nrow(S)) {
    for (j in 1:i) {
      S[i,j] <- sqrt(mahalanobis(x[i,], x[j,], cov=cov))
    }
  }
  S[upper.tri(S)] <- t(S)[upper.tri(S)]
  return(S)
}

dhalfcauchy <- function(x, A, log = FALSE) {
  out <- -A - log(1 + x/(A^2))
  if (!log) out <- exp(out)
  return(out)
}

## Fake Covariate Data
N <- 9
Z <- as.matrix(expand.grid(1:3, 1:3))/10
s.a <- 0.01^2
nu <- 0.005^2
phi <- 0.1
D <- s.a * exp(-mh(Z, diag(phi, ncol(Z))))


## Create data with gammas as the latent variable, y as the observations
## I assume the observational variance is known
T <- 10
F.t <- list()
big <- list()
gs <- matrix(NA, ncol = T, nrow = N)
y <- matrix(NA, ncol = T, nrow = N)
gs[,1] <- runif(N, 0.4, 0.6)
y[,1] <- gs[,1] +  rnorm(N, 0, sqrt((gs[,1] *(1-gs[,1]))/500))
for (t in 2:T) {
  gs[,t] <- rmvnorm(1, gs[,t-1], D)
  y[,t] <- gs[,t] +  rnorm(N, 0, sqrt((gs[,t] *(1-gs[,t]))/500))
  F.t[[t]] <- diag(N)
  big[[t]] <- data.frame(y=y[,t], state=1:N, sigma.sq= sqrt(y[,t] * (1 - y[,t]))/sqrt(500))
}
  

y <- c(y)
dts <- rep(1:T, each = N)
sts <- rep(1:N, times = T)

## first day priors
m.0 <- c(rep(.5, times = N))
C.0 <- diag(c(rep(.1^2, times = N)))

## house effects
b <- rep(0, times = ncol(pollster.contrast))
B <- diag(rep(1/(.1^2), times = ncol(pollster.contrast)))

## common covariance
sigma.a <- 0.01^2
c.0 <- 1
d.0 <- 0.0005
tune.a <- 0.2

phi <- rep(.1, times = ncol(Z))
e.0 <- 1
f.0 <- 0.01
tune.phi <- c(0.5, 0.8, 0.8)##rep(0.5, ncol(Z))

## nugget 
nu <- 0.005^2
g.0 <- 1
h.0 <- 0.0005
tune.nu <- .25

tune <- diag(c(tune.a, tune.nu, tune.phi))

out <- polling.gibbs.kal(y = y, dates = dts, big = big,states = sts, F = F.t, Z = Z,
                  m.0 =m.0, c.0 = c.0, d.0 = d.0, e.0 = e.0, f.0 = f.0, 
                  C.0 = C.0, g.0 = g.0, h.0 = h.0,
                  starts = 0, phi = phi, sigma.a = sigma.a,
                         beta = NULL, nu = nu,
                         tune = tune, tune.phi = tune.phi, A = 5,
                  iters = 50, burnin=0,thin=1)