Run simulation to show AR(1) coefficient bias in small samples.

AR(1) Coefficient Bias

## Prep workspace
options(width=120, digits=6, digits.secs=6)
rm(list=ls())

library(foreign)
library(grid)
library(plyr)
library(ggplot2)
library(tikzDevice)
print(options('tikzLatexPackages'))
options(tikzLatexPackages =
        c("\\usepackage{tikz}\n",
          "\\usepackage[active,tightpage,psfixbb]{preview}\n",
          "\\PreviewEnvironment{pgfpicture}\n",
          "\\setlength\\PreviewBorder{0pt}\n",
          "\\usepackage{amsmath}\n",
          "\\usepackage{xfrac}\n"
          )
        )
setTikzDefaults(overwrite = FALSE)
print(options('tikzLatexPackages'))
library(reshape)
library(vars)
library(scales)   
library(zoo)   

library(foreach)
library(doMC)
registerDoMC(8)



## Define directories
scl.str.DAT_DIR  <- "~/Dropbox/research/distant_speculators_and_mispricing/data/"
scl.str.FIG_DIR  <- "~/Dropbox/research/distant_speculators_and_mispricing/figures/"



## Define global parameters
set.seed(12345)

vec.int.T     <- c(21, 60, 100)
scl.int.T_LEN <- length(vec.int.T)

scl.int.BETA_LEN <- 100
vec.flt.BETA     <- seq(-1, 1, length.out = scl.int.BETA_LEN)

scl.int.ITER <- 10^4







## Define equlibrium functions
if (TRUE == TRUE) {
    mat.flt.SIM <- foreach(b = 1:scl.int.BETA_LEN, .combine='rbind') %:% 
        foreach(t = 1:scl.int.T_LEN, .combine='rbind') %dopar% {
            print(b)

            vec.flt.BETA_HAT  <- rep(NA, scl.int.ITER)
            vec.flt.PRED_BIAS <- rep(NA, scl.int.ITER)

            vec.flt.COV_TRUE <- rep(NA, scl.int.ITER)
            vec.flt.COV_EST  <- rep(NA, scl.int.ITER)
            
            vec.flt.VAR_TRUE <- rep(NA, scl.int.ITER)
            vec.flt.VAR_EST  <- rep(NA, scl.int.ITER)
            
            for (i in 1:scl.int.ITER) {
                vec.flt.X    <- rep(NA, vec.int.T[t])
                vec.flt.X[1] <- 0
                vec.flt.ERR  <- rnorm(vec.int.T[t], 0, 1)
                for (s in 2:vec.int.T[t]) {
                    vec.flt.X[s] <- vec.flt.BETA[b] * vec.flt.X[s-1] + vec.flt.ERR[s]
                }
                obj.lm.RESULT        <- lm(vec.flt.X[2:vec.int.T[t]] ~ vec.flt.X[1:(vec.int.T[t]-1)])
                vec.flt.BETA_HAT[i]  <- summary(obj.lm.RESULT)$coef[2,1]

                scl.flt.MU_X_HAT     <- mean(vec.flt.X[1:(vec.int.T[t]-1)])
                scl.flt.MU_Y_HAT     <- mean(vec.flt.X[2:vec.int.T[t]])
                scl.flt.NUMER        <- 1/(vec.int.T[t]-1) * sum( (vec.flt.ERR[2:vec.int.T[t]] - scl.flt.MU_Y_HAT) * (scl.flt.MU_X_HAT - vec.flt.X[1:(vec.int.T[t]-1)]) )
                scl.flt.DENOM        <- 1/(vec.int.T[t]-1) * sum( (scl.flt.MU_X_HAT - vec.flt.X[1:(vec.int.T[t]-1)]) * (scl.flt.MU_X_HAT - vec.flt.X[1:(vec.int.T[t]-1)]) )                
                vec.flt.PRED_BIAS[i] <- - scl.flt.NUMER/scl.flt.DENOM

                vec.flt.COV_TRUE[i] <- 1/(vec.int.T[t]-1) * sum(vec.flt.ERR[2:vec.int.T[t]] * (0 - vec.flt.X[1:(vec.int.T[t]-1)]))
                vec.flt.COV_EST[i]  <- scl.flt.NUMER

                vec.flt.VAR_TRUE[i] <- 1/(vec.int.T[t]-1) * sum( vec.flt.X[1:(vec.int.T[t]-1)] * vec.flt.X[1:(vec.int.T[t]-1)] )                
                vec.flt.VAR_EST[i]  <- scl.flt.DENOM
                
                
            }
            return(c(vec.flt.BETA[b], vec.int.T[t],
                     mean(vec.flt.BETA_HAT) - vec.flt.BETA[b], mean(vec.flt.PRED_BIAS),
                     mean(vec.flt.COV_EST),
                     mean(vec.flt.VAR_EST) - mean(vec.flt.VAR_TRUE)
                     )
                   )
        }
    mat.df.SIM          <- as.data.frame(mat.flt.SIM)
    names(mat.df.SIM)   <- c("beta", "T", "Actual", "Predicted", "covBias", "varBias")
    save(mat.df.SIM, file = paste(scl.str.DAT_DIR, "data--ar1-coef-bias-prediction.Rdata", sep = ""))
}
load(file = paste(scl.str.DAT_DIR, "data--ar1-coef-bias-prediction.Rdata", sep = ""))











## Create beta bias plot
mat.df.PLOT        <- mat.df.SIM[, c("beta", "T", "Actual")]
names(mat.df.PLOT) <- c("rho", "T", "value")

mat.df.PLOT$T    <- factor(mat.df.PLOT$T,
                            levels = c(21, 60, 100),
                            labels = c("$T = 21$", "$T = 60$", "$T = 100$")
                            )


theme_set(theme_bw())

scl.str.RAW_FILE <- "plot--ar1-coef-bias"
scl.str.TEX_FILE <- paste(scl.str.RAW_FILE,'.tex',sep='')
scl.str.PDF_FILE <- paste(scl.str.RAW_FILE,'.pdf',sep='')
scl.str.AUX_FILE <- paste(scl.str.RAW_FILE,'.aux',sep='')
scl.str.LOG_FILE <- paste(scl.str.RAW_FILE,'.log',sep='')

tikz(file = scl.str.TEX_FILE, height = 2.25, width = 7, standAlone=TRUE)    
obj.gg2.PLOT <- ggplot(data = mat.df.PLOT)                       
obj.gg2.PLOT <- obj.gg2.PLOT + scale_colour_brewer(palette="Set1")
obj.gg2.PLOT <- obj.gg2.PLOT + geom_hline(yintercept = 0,
                                          size       = 0.50,
                                          linetype   = 4
                                          )
obj.gg2.PLOT <- obj.gg2.PLOT + geom_path(aes(x = rho,
                                             y = value
                                             ),
                                         size = 1.25
                                         )
obj.gg2.PLOT <- obj.gg2.PLOT + xlab("$\\rho$")
obj.gg2.PLOT <- obj.gg2.PLOT + ylab("$\\mathrm{E}(\\widehat{\\rho}) - \\rho$")
obj.gg2.PLOT <- obj.gg2.PLOT + scale_x_continuous(limits = c(-1, 1),
                                                  breaks = c(-2/3, -1/3, 0, 1/3, 2/3),
                                                  labels = c("$-\\sfrac{2}{3}$", "$-\\sfrac{1}{3}$", "$0$", "$\\sfrac{1}{3}$", "$\\sfrac{2}{3}$")
                                                  )
obj.gg2.PLOT <- obj.gg2.PLOT + facet_wrap(~ T, nrow = 1)
obj.gg2.PLOT <- obj.gg2.PLOT + theme(plot.margin       = unit(c(1,0.15,0.15,0.15), "lines"),
                                     axis.text         = element_text(size = 8),
                                     axis.title        = element_text(size = 10),
                                     plot.title        = element_text(vjust = 1.75),                                     
                                     panel.grid.minor  = element_blank(),
                                     legend.position   = c(0.90, 0.275),
                                     legend.background = element_blank()
                                     )
obj.gg2.PLOT <- obj.gg2.PLOT + ggtitle("Bias in Estimated $\\mathrm{AR}(1)$ Slope Coefficient")    
print(obj.gg2.PLOT)
dev.off()

system(paste('lualatex', file.path(scl.str.TEX_FILE)), ignore.stdout = TRUE)
system(paste('rm ', scl.str.TEX_FILE, sep = ''))
system(paste('mv ', scl.str.PDF_FILE, ' ', scl.str.FIG_DIR, sep = ''))
system(paste('rm ', scl.str.AUX_FILE, sep = ''))
system(paste('rm ', scl.str.LOG_FILE, sep = ''))