Code simulating the accuracy of the binomial approximation to the normal distribution via Berry-Esseen.

binom-approx-to-normal.R

## Prep workspace
rm(list=ls())
library(foreign)
library(grid)
library(plyr)
library(ggplot2)
library(tikzDevice)
library(reshape)
library(vars)
library(scales)   
library(zoo)   


scl.str.DAT_DIR  <- "~/Dropbox/research/trading_on_coincidences/data/"
scl.str.FIG_DIR  <- "~/Dropbox/research/trading_on_coincidences/figures/"








## Define simulation parameters
vec.int.H <- c(10^2, 10^3, 10^4)

scl.flt.PI    <- 1/100
scl.flt.DELTA <- 1

scl.int.ITER <- 1000






## Simulate sum of binomial shocks
mat.flt.DATA <- matrix(rep(NA, scl.int.ITER * 3), nrow = scl.int.ITER)

for (h in 1:3) {
    for (i in 1:scl.int.ITER) {
        scl.int.H     <- vec.int.H[h]
        scl.flt.X_POS <- sum(scl.flt.DELTA/sqrt(scl.int.H) * as.numeric(runif(scl.int.H,0,1) < scl.flt.PI))
        scl.flt.X_NEG <- - sum(scl.flt.DELTA/sqrt(scl.int.H) * as.numeric(runif(scl.int.H,0,1) < scl.flt.PI))    

        mat.flt.DATA[i, h] <- (scl.flt.X_POS + scl.flt.X_NEG)/(scl.flt.DELTA * sqrt(2 * scl.flt.PI * (1 - scl.flt.PI)))
    }
}





## Compute empirical CDF
vec.flt.S     <- seq(-4, 4, by = 0.01)
scl.int.S_LEN <- length(vec.flt.S)
mat.dfm.PLOT  <- data.frame(s = rep(vec.flt.S, 3),
                            H = sort(rep(vec.int.H, scl.int.S_LEN)),
                            F = NA
                            )
for (h in 1:3) {
    for (s in 1:scl.int.S_LEN) {
        mat.dfm.PLOT[(mat.dfm.PLOT$H == vec.int.H[h]) & (mat.dfm.PLOT$s == vec.flt.S[s]), ]$F <- sum(mat.flt.DATA[,h] <= vec.flt.S[s])/scl.int.ITER
    }
}
mat.dfm.PLOT$H2  <- factor(mat.dfm.PLOT$H, labels = c("$H = 10^2$", "$H = 10^3$", "$H = 10^4$"))
mat.dfm.PLOT$PHI <- pnorm(mat.dfm.PLOT$s, 0, 1)

mat.dfm.LBL <- data.frame(s   = c(-0.55, 2),
                          H   = c(10^2,10^2),
                          F   = c(0.55,0.80),
                          lbl = c("$\\Phi(x)$","$F_H(x)$"),
                          col = c(brewer_pal(palette = "Set1")(2)[2], brewer_pal(palette = "Set1")(2)[1])
                          )
mat.dfm.LBL$H2  <- factor(mat.dfm.LBL$H, labels = c("$H = 10^2$"))

mat.dfm.MAX_ERR <- ddply(mat.dfm.PLOT,
                         c("H"),
                         function(X)c(max(abs(X$F - X$PHI)))
                         )
vec.flt.MAX_ERR <- mat.dfm.MAX_ERR[,2]
mat.dfm.MAX_ERR <- data.frame(s   = rep(1.5, 3),
                              H   = c(10^2,10^3,10^4),
                              F   = rep(0.05, 3),
                              lbl =
                              c(paste("$\\max_x | F_H(x) - \\Phi(x) | = ", round(vec.flt.MAX_ERR[1],2), "$", sep=""),
                                paste("$\\max_x | F_H(x) - \\Phi(x) | = ", round(vec.flt.MAX_ERR[2],2), "$", sep=""),
                                paste("$\\max_x | F_H(x) - \\Phi(x) | = ", round(vec.flt.MAX_ERR[3],2), "$", sep="")
                                )
                              )
mat.dfm.MAX_ERR$H2  <- factor(mat.dfm.MAX_ERR$H, labels = c("$H = 10^2$", "$H = 10^3$", "$H = 10^4$"))


vec.flt.PRED_ERR <- 0.50/sqrt(2 * vec.int.H * scl.flt.PI * (1 - scl.flt.PI))
mat.dfm.PRED_ERR <- data.frame(s   = rep(1.5, 3),
                               H   = c(10^2,10^3,10^4),
                               F   = rep(0.20, 3),
                               lbl =
                               c(paste("$\\frac{0.50}{\\sqrt{2 \\cdot H \\cdot \\pi \\cdot (1 - \\pi)}} = ", round(vec.flt.PRED_ERR[1],3), "$", sep=""),
                                 paste("$\\frac{0.50}{\\sqrt{2 \\cdot H \\cdot \\pi \\cdot (1 - \\pi)}} = ", round(vec.flt.PRED_ERR[2],3), "$", sep=""),
                                 paste("$\\frac{0.50}{\\sqrt{2 \\cdot H \\cdot \\pi \\cdot (1 - \\pi)}} = ", round(vec.flt.PRED_ERR[3],3), "$", sep="")
                                 )
                               )
mat.dfm.PRED_ERR$H2  <- factor(mat.dfm.PRED_ERR$H, labels = c("$H = 10^2$", "$H = 10^3$", "$H = 10^4$"))



## Create realized dividend plot
theme_set(theme_bw())

scl.str.RAW_FILE <- 'normal-approx-to-binom'
scl.str.TEX_FILE <- paste(scl.str.RAW_FILE,'.tex',sep='')
scl.str.PDF_FILE <- paste(scl.str.RAW_FILE,'.pdf',sep='')
scl.str.PNG_FILE <- paste(scl.str.RAW_FILE,'.png',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, width = 7, standAlone=TRUE)

obj.gg2.PLOT <- ggplot()
obj.gg2.PLOT <- obj.gg2.PLOT + scale_colour_brewer(palette="Set1")
obj.gg2.PLOT <- obj.gg2.PLOT + geom_path(data = mat.dfm.PLOT,
                                         aes(x     = s,
                                             y     = PHI,
                                             group = H2
                                             ),
                                         size   = 0.25,
                                         colour = brewer_pal(palette = "Set1")(2)[1]
                                         )
obj.gg2.PLOT <- obj.gg2.PLOT + geom_point(data = mat.dfm.PLOT,
                                          aes(x     = s,
                                              y     = F,
                                              group = H2
                                              ),
                                          size   = 0.50,
                                          colour = brewer_pal(palette = "Set1")(2)[2]
                                          )
obj.gg2.PLOT <- obj.gg2.PLOT + geom_text(data = mat.dfm.LBL,
                                         aes(x      = s,
                                             y      = F,
                                             label  = lbl,
                                             group  = H2,
                                             colour = col
                                             ),
                                         size = 3
                                         )
obj.gg2.PLOT <- obj.gg2.PLOT + geom_text(data = mat.dfm.MAX_ERR,
                                         aes(x      = s,
                                             y      = F,
                                             label  = lbl,
                                             group  = H2
                                             ),
                                         size = 2
                                         )
obj.gg2.PLOT <- obj.gg2.PLOT + geom_text(data = mat.dfm.PRED_ERR,
                                         aes(x      = s,
                                             y      = F,
                                             label  = lbl,
                                             group  = H2
                                             ),
                                         size = 2
                                         )
obj.gg2.PLOT <- obj.gg2.PLOT + xlab("")
obj.gg2.PLOT <- obj.gg2.PLOT + ylab("$F_H(x)$")
obj.gg2.PLOT <- obj.gg2.PLOT + facet_wrap(~ H2, nrow = 1)
obj.gg2.PLOT <- obj.gg2.PLOT + theme(plot.margin     = unit(c(1,1,-1,0), "lines"),
                                     plot.title      = element_text(vjust = 1.75),
                                     legend.position = "none"
                                     )

obj.gg2.PLOT <- obj.gg2.PLOT + ggtitle("Normal Approximation to Binomial Distribution")

print(obj.gg2.PLOT)
dev.off()

system(paste('pdflatex', file.path(scl.str.TEX_FILE)), ignore.stdout = TRUE)
system(paste('convert -density 600', file.path(scl.str.PDF_FILE), ' ', file.path(scl.str.PNG_FILE)))
system(paste('mv ', scl.str.PNG_FILE, ' ', scl.str.FIG_DIR, sep = ''))
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 = ''))