# How much Nickell bias given infinite N, no X, # and various T and phi? # # Christopher Adolph faculty.washington.edu/cadolph # 26 June 2015 rm(list=ls()) library(tile) # compute Nickell bias given the true AR(1) coefficient # and the number of time periods T # (from Nickell (1981), equation 18) nickell <- function(phi, T) { A <- 1 - 1/T * (1-phi^T)/(1-phi) bias <- -(1+phi)/(T-1) * A * (1 - 2*phi/((1-phi)*(T-1)) * A)^(-1) bias } # Scenarios to calculate (all combinations) T <- 2:200 phi <- c(0.2, 0.5, 0.8, 0.95) # Storage array arrayOut <- array(NA, c(length(T),4,length(phi))) # Compute array of all combinations of T and phi for (i in 1:length(phi)) { arrayOut[,,i] <- cbind(rep(phi[i], length(T)), T, nickell(phi[i], T), 100*abs(nickell(phi[i], T)/phi[i])) } # Plot results using tile graphics collectedTraces <- vector("list", length(phi)) for (i in 1:length(phi)) collectedTraces[[i]] <- linesTile(x=T, y=arrayOut[,3,i], plot=1) phiLabs <- textTile(x=rep(1.5,length(phi)), y=arrayOut[1,3,], labels=paste0("phi=",phi), plot=1) legend <- textTile(x=rep(85,3), y=c(-0.645,-0.70,-0.755), labels=c("Bias in fixed effects estimator", "with a lagged outcome", "and no other covariates" ), col=c("black", "black", "black"), cex=0.8, plot=1) xat <- c(2, 3, 4, 5, 7, 10, 15, 20, 30, 50, 100, 200) tile(collectedTraces, phiLabs, legend, limits=c(1,220,-1.0,0.05), xaxis=list(log=TRUE, at=xat), yaxis=list(major=FALSE), xaxistitle=list(labels="Periods modeled (T)", y=0.1), yaxistitle=list(labels="Bias of lagged DV (phi)", x=0), gridlines=list(type="y"), height=list(plot="golden"), output=list(file="asympPhiBias",width=8) )