# Helper functions for Monte Carlo of FE estimator # # Christopher Adolph faculty.washington.edu/cadolph # 27 June 2015 # Simulate panel data and run through the within estimator simulatePanel <- function(sims, N, T, phi, beta, sigma, sigmaX, phiX, sigmaAlpha, burnin=100) { require(plm) phihat <- betahat <- rep(NA, sims) res <- vector("list", sims) for (isim in 1:sims) { # Make data # Loop over units y <- x <- units <- periods <- NULL for (iobs in 1:N) { alpha <- rnorm(1, mean=0, sd=sigmaAlpha) currX <- makeTS(n=T+burnin, ar=phiX, varY=sigmaX, burnin=0) x <- c(x, currX[(burnin+1):(burnin+T)]) y <- c(y, alpha + makeTS(n=T, ar=phi, x=currX, beta=beta, varY=sigma, burnin=burnin)) units <- c(units, rep(iobs, T)) periods <- c(periods, 1:T) } curdata <- as.data.frame(cbind(units, periods, y, x)) curdata <- pdata.frame(curdata, c("units","periods")) # Estimate beta and phi with FE estimator res[[isim]] <- plm(y~lag(y)+x, curdata, model="within") phihat[isim] <- coef(res[[isim]])[1] betahat[isim] <- coef(res[[isim]])[2] } # Results bias <- list(beta=mean(betahat - beta), phi=mean(phihat - phi), longrun=mean(betahat/(1-phihat) - beta/(1-phi))) rmse <- list(beta=sqrt(mean( (betahat - beta)^2 )), phi=sqrt(mean( (phihat - phi)^2 )), longrun=sqrt(mean( (betahat/(1-phihat) - beta/(1-phi))^2 ))) list(bias=bias, rmse=rmse) } # Simple univariate time series simulator: # Allows deterministic trends, seasonality, and ARMA(p,q) # process, including potentially nonstationary and/or # noninvertible processes, as well as covariates # # Christopher Adolph (faculty.washington.edu/cadolph) # 30 May 2015 # # Usage: # # y <- makeTS(n = 1000, # Simulate 1000 periods # ar = c(0.67), # AR(1); phi = 0.67 # ma = NULL, # MA(0) # trend = 0.01, # Determ. trend (+0.01 units/period) # seasons = c(0,0,0, 0,0,0, 0,0,0, 0,0,0), # 12 period cycle # adjust = c("level"), # Additive cycles; alternative is "factor" # # which multiplies last round # x = NULL, # n x k matrix of exogenous covariates # beta = NULL, # k-vector of parameters relating x to y # varY = 1, # Error term variance # initY = 0, # Prior level of time series, for lags # initE = 0, # Prior error, for lags # burnin = 10 # Discard this many iterations before saving # ) # makeTS <- function(n=1, ar=NULL, ma=NULL, trend=NULL, seasons=NULL, adjust=c("factor", "level"), x=NULL, beta=NULL, varY=1, initY=0, initE=0, burnin=0 ) { nseasons <- length(seasons) doS <- length(seasons)>0 doT <- length(trend)>0 doA <- length(ar)>0 doM <- length(ma)>0 if (!is.null(x)) {doX <- TRUE; x <- as.matrix(x)} else {doX <- FALSE} y <- rep(NA, n + burnin) while (length(initY)