Skip to content

Instantly share code, notes, and snippets.

@alexchinco

alexchinco/plot--inferring-horizon-from-volume--12jul2016.R

Created July 13, 2016 14:35
Show Gist options
  • Save alexchinco/026529ed8bbaecb824ac12fea7a534a5 to your computer and use it in GitHub Desktop.
Save alexchinco/026529ed8bbaecb824ac12fea7a534a5 to your computer and use it in GitHub Desktop.
Download ZIP
Code for blog post on inferring traders' investment horizons from aggregate trading volume.
Raw
plot--inferring-horizon-from-volume--12jul2016.R
###########################################################################
###########################################################################
## @purpose: Verify that the wavelet-variance estimator picks up traders'
## investment horizons in simulated data.
## ------------------------------------------------------------------------
## @author: Alex Chinco
## @date: 12-JUL-2016
###########################################################################
###########################################################################
###########################################################################
###########################################################################
## @section: Prep workspace
###########################################################################
###########################################################################
options(width=200, digits=6, digits.secs=6)
rm(list=ls())
library(foreign)
library(grid)
library(plyr)
library(ggplot2)
library(reshape)
library(vars)
library(scales)
library(wmtsa)
library(foreach)
library(doMC)
registerDoMC(8)
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{pxfonts}\n",
"\\usepackage{bm}\n",
"\\usepackage{xfrac}\n"
)
)
setTikzDefaults(overwrite = FALSE)
print(options('tikzLatexPackages'))
scl.str.DAT_DIR <- "~/Dropbox/research/fast_trading_priced_noise/data/"
scl.str.FIG_DIR <- "~/Dropbox/research/fast_trading_priced_noise/figures/"
set.seed(12356)
###########################################################################
###########################################################################
## @section: Set global parameters
###########################################################################
###########################################################################
scl.flt.BET <- 0.90
scl.flt.SIG_EP <- 1
scl.flt.GAM <- 0.75
scl.flt.SIG_XI <- 1
scl.flt.F0 <- 0
scl.flt.KAP <- 0.25
scl.flt.RHO <- 1
scl.flt.TET <- 0.00
###########################################################################
###########################################################################
## @section: Define global functions
###########################################################################
###########################################################################
fun.MU_H <- function(BET, GAM, H) {
scl.flt.OUT <- BET * (1 - GAM)^(H - 1)
return(scl.flt.OUT)
}
fun.VAR_H <- function(BET, SIG_EP, GAM, SIG_XI, H) {
scl.flt.OUT <- SIG_EP^2
if (H >= 2) {
for (h in 1:(H-1)) {
scl.flt.OUT <- scl.flt.OUT + BET^2 * (1 - GAM)^(2 * ((H-1) - h)) * SIG_XI^2
}
}
return(scl.flt.OUT)
}
fun.CREATE_DATA <- function(BET, SIG_EP, GAM, SIG_XI, T, F0) {
vec.flt.EP <- rnorm(T, 0, SIG_EP)
vec.flt.XI <- rnorm(T, 0, SIG_XI)
vec.flt.R <- rep(NA, T)
vec.flt.F <- rep(NA, T)
vec.flt.dF <- rep(NA, T)
for (tPlus1 in 1:T) {
if (tPlus1 == 1) {
scl.flt.Ft <- F0
} else {
scl.flt.Ft <- vec.flt.F[tPlus1-1]
}
vec.flt.R[tPlus1] <- BET * scl.flt.Ft + vec.flt.EP[tPlus1]
vec.flt.dF[tPlus1] <- - GAM * scl.flt.Ft + vec.flt.XI[tPlus1]
vec.flt.F[tPlus1] <- scl.flt.Ft + vec.flt.dF[tPlus1]
}
mat.df.DATA <- data.frame(t = seq(1, T),
fLag = c(F0, vec.flt.F[-T]),
f = vec.flt.F,
df = vec.flt.dF,
r = vec.flt.R,
ep = vec.flt.EP,
xi = vec.flt.XI
)
return(mat.df.DATA)
}
###########################################################################
###########################################################################
## @section: Simulate data and compute moments
###########################################################################
###########################################################################
if (TRUE == FALSE) {
scl.int.T <- 2^18
vec.int.H <- 2^c(0,1,2,4,6,8)
scl.int.H_LEN <- length(vec.int.H)
mat.flt.PLOT <- foreach(h=1:scl.int.H_LEN, .combine = "rbind") %do% {
## for (h in 1:scl.int.H_LEN) {
scl.int.H <- vec.int.H[h]
print(paste(scl.int.H, ": ", round(Sys.time()), sep = ""))
mat.flt.OUTTER <- foreach(i=1:10, .combine = "rbind") %dopar% {
## for (i in 1:10) {
mat.df.TEMP <- fun.CREATE_DATA(scl.flt.BET, scl.flt.SIG_EP, scl.flt.GAM, scl.flt.SIG_XI, scl.int.T, scl.flt.F0)
scl.int.NUM_INT <- scl.int.T/scl.int.H
mat.df.TEMP$i <- sort(rep(seq(1, scl.int.NUM_INT), scl.int.H))
vec.int.SAMPLE <- sample(scl.int.NUM_INT, 2^10)
mat.df.TEMP <- mat.df.TEMP[mat.df.TEMP$i %in% vec.int.SAMPLE, ]
mat.df.TEMP <- ddply(mat.df.TEMP, c("i"), function(X)c(head(X$fLag,1), tail(X$r,1)))
names(mat.df.TEMP) <- c("i", "ft", "rtplush")
obj.lm.RESULTS <- lm(rtplush ~ 1 + ft, data = mat.df.TEMP)
scl.flt.MU_HAT <- obj.lm.RESULTS$coef[2]
scl.flt.VAR_HAT <- var(obj.lm.RESULTS$residuals)
vec.flt.INNER <- c(scl.flt.MU_HAT, scl.flt.VAR_HAT)
return(vec.flt.INNER)
}
vec.flt.OUTTER <- apply(mat.flt.OUTTER, 2, mean)
scl.flt.MU <- fun.MU_H(scl.flt.BET, scl.flt.GAM, scl.int.H)
scl.flt.VAR <- fun.VAR_H(scl.flt.BET, scl.flt.SIG_EP, scl.flt.GAM, scl.flt.SIG_XI, scl.int.H)
vec.flt.OUTTER <- c(scl.int.H, vec.flt.OUTTER, scl.flt.MU, scl.flt.VAR)
return(vec.flt.OUTTER)
}
mat.df.PLOT <- as.data.frame(mat.flt.PLOT)
names(mat.df.PLOT) <- c("h", "muHat", "varHat", "mu", "var")
mat.df.PLOTA <- mat.df.PLOT[, c("h", "mu", "var")]
names(mat.df.PLOTA) <- c("h", "$\\mu_{H,t}/f_t$", "$\\sigma_H^2$")
mat.df.PLOTA <- melt(mat.df.PLOTA, c("h"))
mat.df.PLOTA$type <- "Predicted"
mat.df.PLOTB <- mat.df.PLOT[, c("h", "muHat", "varHat")]
names(mat.df.PLOTB) <- c("h", "$\\mu_{H,t}/f_t$", "$\\sigma_H^2$")
mat.df.PLOTB <- melt(mat.df.PLOTB, c("h"))
mat.df.PLOTB$type <- "Realized"
mat.df.PLOT <- rbind(mat.df.PLOTA, mat.df.PLOTB)
theme_set(theme_bw())
scl.str.RAW_FILE <- "plot--predicted-vs-realized-moments--12jul2016"
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(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 = 1.25,
linetype = 4,
alpha = 0.75
)
obj.gg2.PLOT <- obj.gg2.PLOT + geom_path(aes(x = log(h,2),
y = value,
group = type,
colour = type,
linetype = type
),
size = 2.25,
alpha = 0.75
)
obj.gg2.PLOT <- obj.gg2.PLOT + xlab("$\\log_2(H)$")
obj.gg2.PLOT <- obj.gg2.PLOT + ylab("")
obj.gg2.PLOT <- obj.gg2.PLOT + labs(colour = "", linetype = "")
obj.gg2.PLOT <- obj.gg2.PLOT + coord_cartesian(xlim = c(0, 8), ylim = c(0, 2))
obj.gg2.PLOT <- obj.gg2.PLOT + scale_x_continuous(breaks = c(0,1,2,4,6,8))
obj.gg2.PLOT <- obj.gg2.PLOT + scale_y_continuous(breaks = c(0, 0.50, 1.00, 1.50, 2))
obj.gg2.PLOT <- obj.gg2.PLOT + facet_wrap(~variable, ncol = 2)
obj.gg2.PLOT <- obj.gg2.PLOT + theme(plot.margin = unit(c(1,0.15,0.15,-1), "lines"),
axis.text = element_text(size = 10),
axis.title = element_text(size = 10),
plot.title = element_text(vjust = 1.75),
panel.grid.minor = element_blank(),
legend.position = c(0.90, 0.40),
legend.background = element_blank()
)
obj.gg2.PLOT <- obj.gg2.PLOT + ggtitle("Predicted vs. Realized Moments")
print(obj.gg2.PLOT)
dev.off()
system(paste('lualatex', file.path(scl.str.TEX_FILE)), ignore.stdout = TRUE)
system(paste('mv ', scl.str.PDF_FILE, ' ', scl.str.FIG_DIR, sep = ''))
system(paste('rm ', scl.str.TEX_FILE, sep = ''))
system(paste('rm ', scl.str.AUX_FILE, sep = ''))
system(paste('rm ', scl.str.LOG_FILE, sep = ''))
}
###########################################################################
###########################################################################
## @section: Simulate trading-volume time series when all trading occurs at
## a single horizon and compute the horizon-specific variance.
###########################################################################
###########################################################################
if (TRUE == FALSE) {
scl.int.T <- 2^12
vec.int.H <- 2^c(0,2,4,6,8)
if (TRUE == FALSE) {
mat.flt.PLOT <- foreach(h=1:length(vec.int.H), .combine = "rbind") %do% {
## for (h in 1:5) {
scl.int.H <- vec.int.H[h]
print(paste(scl.int.H, ": ", round(Sys.time()), sep = ""))
scl.flt.MU <- fun.MU_H(scl.flt.BET, scl.flt.GAM, scl.int.H)
scl.flt.VAR <- fun.VAR_H(scl.flt.BET, scl.flt.SIG_EP, scl.flt.GAM, scl.flt.SIG_XI, scl.int.H)
mat.flt.OUTTER <- foreach(i=1:100, .combine = "rbind") %dopar% {
## for (i in 1:100) {
mat.df.DATA <- fun.CREATE_DATA(scl.flt.BET, scl.flt.SIG_EP, scl.flt.GAM, scl.flt.SIG_XI, scl.int.T, scl.flt.F0)
scl.int.NUM_INT <- scl.int.T/scl.int.H
mat.df.DATA$i <- sort(rep(seq(1, scl.int.NUM_INT), scl.int.H))
mat.df.TEMP <- ddply(mat.df.DATA, c("i"), function(X)head(X$fLag,1))
names(mat.df.TEMP) <- c("i", "fLag")
mat.df.TEMP$x <- (scl.flt.MU * mat.df.TEMP$fLag)/(scl.flt.RHO^2 * scl.flt.VAR - scl.flt.KAP/scl.int.H)
mat.df.DATA <- mat.df.DATA[, c("t", "i")]
mat.df.TEMP <- mat.df.TEMP[, c("i", "x")]
mat.df.DATA <- merge(mat.df.DATA, mat.df.TEMP, by = c("i"))
vec.flt.X <- as.numeric(mat.df.DATA$x)
obj.wmtsa.WVAR <- wavVar(vec.flt.X, wavelet = "haar")
vec.flt.WVAR <- summary(obj.wmtsa.WVAR)$vmat[2,]
scl.flt.VAR <- sum(vec.flt.WVAR)
vec.flt.INNER <- c(scl.int.H, vec.flt.WVAR[1:9]/scl.flt.VAR)
return(vec.flt.INNER)
}
vec.flt.OUTTER <- apply(mat.flt.OUTTER, 2, mean)
return(vec.flt.OUTTER)
}
mat.df.PLOT <- as.data.frame(mat.flt.PLOT)
save(mat.df.PLOT, file = paste(scl.str.DAT_DIR, "data--all-trading-at-single-horizon--12jul2016.Rdata", sep = ""))
}
load(paste(scl.str.DAT_DIR, "data--all-trading-at-single-horizon--12jul2016.Rdata", sep = ""))
names(mat.df.PLOT) <- c("h", seq(0,8))
mat.df.PLOT <- melt(mat.df.PLOT, c("h"))
mat.df.PLOT$variable <- as.numeric(as.character(mat.df.PLOT$variable))
mat.df.PLOT$h <- factor(mat.df.PLOT$h, labels = log(vec.int.H,2))
theme_set(theme_bw())
scl.str.RAW_FILE <- "plot--all-trading-at-single-horizon--12jul2016"
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 = 3, 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 = 1.25,
linetype = 4,
alpha = 0.75
)
obj.gg2.PLOT <- obj.gg2.PLOT + geom_path(aes(x = variable,
y = value,
group = h,
colour = h,
linetype = h
),
size = 2.25,
alpha = 0.75
)
obj.gg2.PLOT <- obj.gg2.PLOT + xlab("$\\log_2(H)$")
obj.gg2.PLOT <- obj.gg2.PLOT + ylab("$\\mathrm{WVar}_H/\\mathrm{Var}$")
obj.gg2.PLOT <- obj.gg2.PLOT + labs(colour = "$\\log_2(H^\\star)$", linetype = "$\\log_2(H^\\star)$")
obj.gg2.PLOT <- obj.gg2.PLOT + coord_cartesian(xlim = c(0, 8), ylim = c(0, 0.40))
obj.gg2.PLOT <- obj.gg2.PLOT + scale_x_continuous(breaks = c(0,2,4,6,8))
obj.gg2.PLOT <- obj.gg2.PLOT + scale_y_continuous(breaks = c(0, 0.10, 0.20, 0.30, 0.40))
obj.gg2.PLOT <- obj.gg2.PLOT + theme(plot.margin = unit(c(1,0.15,0.15,0.15), "lines"),
axis.text = element_text(size = 10),
axis.title = element_text(size = 10),
plot.title = element_text(vjust = 1.75),
panel.grid.minor = element_blank(),
legend.background = element_blank()
)
obj.gg2.PLOT <- obj.gg2.PLOT + ggtitle("All Trading at a Single Horizon, $H^\\star$")
print(obj.gg2.PLOT)
dev.off()
system(paste('lualatex', file.path(scl.str.TEX_FILE)), ignore.stdout = TRUE)
system(paste('mv ', scl.str.PDF_FILE, ' ', scl.str.FIG_DIR, sep = ''))
system(paste('rm ', scl.str.TEX_FILE, sep = ''))
system(paste('rm ', scl.str.AUX_FILE, sep = ''))
system(paste('rm ', scl.str.LOG_FILE, sep = ''))
}
###########################################################################
###########################################################################
## @section: Simulate trading-volume time series when trading is split between
## 2 different horizons and compute the horizon-specific variance.
###########################################################################
###########################################################################
if (TRUE == TRUE) {
scl.int.T <- 2^12
scl.int.H_FAST <- 2^2
scl.int.H_SLOW <- 2^6
scl.flt.MU_FAST <- fun.MU_H(scl.flt.BET, scl.flt.GAM, scl.int.H_FAST)
scl.flt.VAR_FAST <- fun.VAR_H(scl.flt.BET, scl.flt.SIG_EP, scl.flt.GAM, scl.flt.SIG_XI, scl.int.H_FAST)
scl.flt.MU_SLOW <- fun.MU_H(scl.flt.BET, scl.flt.GAM, scl.int.H_SLOW)
scl.flt.VAR_SLOW <- fun.VAR_H(scl.flt.BET, scl.flt.SIG_EP, scl.flt.GAM, scl.flt.SIG_XI, scl.int.H_SLOW)
mat.flt.PLOT <- foreach(n=1:1000, .combine = "rbind") %dopar% {
## for (n in 1:100) {
mat.df.DATA <- fun.CREATE_DATA(scl.flt.BET, scl.flt.SIG_EP, scl.flt.GAM, scl.flt.SIG_XI, scl.int.T, scl.flt.F0)
scl.int.NUM_INT_FAST <- scl.int.T/scl.int.H_FAST
mat.df.DATA$iFast <- sort(rep(seq(1, scl.int.NUM_INT_FAST), scl.int.H_FAST))
mat.df.FAST <- ddply(mat.df.DATA, c("iFast"), function(X)head(X$fLag,1))
names(mat.df.FAST) <- c("iFast", "fLag")
mat.df.FAST$xFast <- (scl.flt.MU_FAST * mat.df.FAST$fLag)/(scl.flt.RHO^2 * scl.flt.VAR_FAST - scl.flt.KAP/scl.int.H_FAST)
scl.int.NUM_INT_SLOW <- scl.int.T/scl.int.H_SLOW
mat.df.DATA$iSlow <- sort(rep(seq(1, scl.int.NUM_INT_SLOW), scl.int.H_SLOW))
mat.df.SLOW <- ddply(mat.df.DATA, c("iSlow"), function(X)head(X$fLag,1))
names(mat.df.SLOW) <- c("iSlow", "fLag")
mat.df.SLOW$xSlow <- (scl.flt.MU_SLOW * mat.df.SLOW$fLag)/(scl.flt.RHO^2 * scl.flt.VAR_SLOW - scl.flt.KAP/scl.int.H_SLOW)
mat.df.DATA <- mat.df.DATA[, c("t", "iFast", "iSlow")]
mat.df.FAST <- mat.df.FAST[, c("iFast", "xFast")]
mat.df.SLOW <- mat.df.SLOW[, c("iSlow", "xSlow")]
mat.df.DATA <- merge(mat.df.DATA, mat.df.FAST, by = c("iFast"))
mat.df.DATA <- merge(mat.df.DATA, mat.df.SLOW, by = c("iSlow"))
mat.df.DATA <- mat.df.DATA[order(mat.df.DATA$t), ]
mat.df.DATA$xFast <- 0.50 * mat.df.DATA$xFast/sd(mat.df.DATA$xFast)
mat.df.DATA$xSlow <- 0.50 * mat.df.DATA$xSlow/sd(mat.df.DATA$xSlow)
vec.flt.X <- as.numeric(mat.df.DATA$xFast + mat.df.DATA$xSlow)
obj.wmtsa.WVAR <- wavVar(vec.flt.X, wavelet = "haar")
vec.flt.WVAR <- summary(obj.wmtsa.WVAR)$vmat[2,]
scl.flt.VAR <- sum(vec.flt.WVAR)
vec.flt.OUT <- vec.flt.WVAR[1:9]/scl.flt.VAR
return(vec.flt.OUT)
}
vec.flt.PLOT <- apply(mat.flt.PLOT, 2, mean)
mat.df.PLOT <- data.frame(h = seq(0,8),
value = vec.flt.PLOT
)
theme_set(theme_bw())
scl.str.RAW_FILE <- "plot--trading-at-two-different-horizons--12jul2016"
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 = 3, 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 = 1.25,
linetype = 4,
alpha = 0.75
)
obj.gg2.PLOT <- obj.gg2.PLOT + geom_path(aes(x = h,
y = value
),
size = 2.25,
alpha = 0.75
)
obj.gg2.PLOT <- obj.gg2.PLOT + xlab("$\\log_2(H)$")
obj.gg2.PLOT <- obj.gg2.PLOT + ylab("$\\mathrm{WVar}_H/\\mathrm{Var}$")
obj.gg2.PLOT <- obj.gg2.PLOT + coord_cartesian(xlim = c(0, 8), ylim = c(0, 0.20))
obj.gg2.PLOT <- obj.gg2.PLOT + scale_x_continuous(breaks = c(0,2,4,6,8))
obj.gg2.PLOT <- obj.gg2.PLOT + scale_y_continuous(breaks = c(0, 0.05, 0.10, 0.15, 0.20))
obj.gg2.PLOT <- obj.gg2.PLOT + theme(plot.margin = unit(c(1,0.15,0.15,0.15), "lines"),
axis.text = element_text(size = 10),
axis.title = element_text(size = 10),
plot.title = element_text(vjust = 1.75),
panel.grid.minor = element_blank(),
legend.background = element_blank()
)
obj.gg2.PLOT <- obj.gg2.PLOT + ggtitle("Trading at Two Different Horizons, $\\log_2(H^\\star) = \\{ \\, 2, \\, 6 \\, \\}$")
print(obj.gg2.PLOT)
dev.off()
system(paste('lualatex', file.path(scl.str.TEX_FILE)), ignore.stdout = TRUE)
system(paste('mv ', scl.str.PDF_FILE, ' ', scl.str.FIG_DIR, sep = ''))
system(paste('rm ', scl.str.TEX_FILE, sep = ''))
system(paste('rm ', scl.str.AUX_FILE, sep = ''))
system(paste('rm ', scl.str.LOG_FILE, sep = ''))
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment