'sim' method for 'plm' objects

[ruettenauer]

Dear Yu-Sung Su,

I recently wrote some lines of code to use the sim function with plm objects. I didn't test the code in different scenarios, but it might help to implement the sim function for class plm objects (if this is an option).

sim.plm<-function(object, n.sims=100)
{
  object.class <- class(object)[[1]]
  summ <- summary (object)
  coef <- summ$coef[,1:2,drop=FALSE]
  dimnames(coef)[[2]] <- c("coef.est","coef.sd")
  # sigma.hat <- summ$sigma 
  # TR: define sigma by hand
  NN <- nrow(object$model)
  PP <- nrow(coef)
  sigma.hat <- sqrt(deviance(object) / (NN-PP))
  # TR: end              
  beta.hat <- coef[,1,drop = FALSE]
  # V.beta <- summ$cov.unscaled
  V.beta <- vcov(summ)/sigma.hat^2 # TR: unscale scaled vcov
  # n <- summ$df[1] + summ$df[2]
  # k <- summ$df[1]
  n <- nrow(summ$model) # TR: define n
  k <- nrow(summ$coefficients) # TR: define k
  sigma <- rep (NA, n.sims)
  beta <- array (NA, c(n.sims,k))
  dimnames(beta) <- list (NULL, rownames(beta.hat))
  for (s in 1:n.sims){
    sigma[s] <- sigma.hat*sqrt((n-k)/rchisq(1,n-k))
    beta[s,] <- MASS::mvrnorm(1, beta.hat, V.beta * sigma[s]^2)
  }

  ans <- new("sim",
             coef = beta,
             sigma = sigma)
  return (ans)
}

#### Example ####

library(arm)
library(plm)

data(Cigar)

plm.mod<-plm(sales~ price + pop + pimin + price*pimin,
             model="within", effect="individual",
             index=c("state", "year"), data=Cigar)

summary(plm.mod)

plm.mod.sim<-sim.plm(plm.mod, n.sims=500)

Best, Tobias