Created
August 28, 2015 15:13
-
-
Save alexchinco/e07ded72159c6f3d312f to your computer and use it in GitHub Desktop.
Solution to multiscale noisy-rational-expectations model.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| ## 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{pxfonts}\n", | |
| "\\usepackage{bm}\n", | |
| "\\usepackage{xfrac}\n" | |
| ) | |
| ) | |
| setTikzDefaults(overwrite = FALSE) | |
| print(options('tikzLatexPackages')) | |
| library(reshape) | |
| library(vars) | |
| library(scales) | |
| library(zoo) | |
| library(splines) | |
| library(MASS) | |
| library(wmtsa) | |
| library(foreach) | |
| library(doMC) | |
| registerDoMC(8) | |
| ## Define directories | |
| scl.str.DAT_DIR <- "~/Dropbox/research/fast_trading_priced_noise/data/" | |
| scl.str.FIG_DIR <- "~/Dropbox/research/fast_trading_priced_noise/figures/" | |
| ## Define global parameters | |
| set.seed(123456) | |
| scl.int.ITER <- 10^4 | |
| scl.int.T <- 10^4 | |
| vec.flt.SIG_E <- sort(runif(10^5, 0, 1.5)) | |
| ## vec.flt.SIG_E <- seq(0, 1.50, by = 0.00001)[-1] | |
| scl.flt.MU_V <- 1 | |
| scl.flt.SIG_V <- 1 | |
| scl.flt.SIG_Z <- 1 | |
| ## Define equlibrium functions | |
| fun.PSI <- function(SIG_E, SIG_V) { | |
| OUT <- SIG_E^2/(SIG_E^2 + SIG_V^2) | |
| return(OUT) | |
| } | |
| fun.OMG <- function(SIG_Z, GAM_V, PSI, SIG_V) { | |
| OUT <- (SIG_Z^2/GAM_V^2)/(PSI * SIG_V^2 + SIG_Z^2/GAM_V^2) | |
| return(OUT) | |
| } | |
| fun.BET_V <- function(LAM_A, GAM_V, PSI, OMG) { | |
| OUT <- 1/(2 * LAM_A) * (1 - LAM_A * GAM_V * (1 - PSI) * OMG) | |
| return(OUT) | |
| } | |
| fun.BET_S <- function(GAM_V, OMG) { | |
| OUT <- - GAM_V/2 * (1 - OMG) | |
| return(OUT) | |
| } | |
| fun.GAM_V <- function(LAM_A, BET_V, BET_S) { | |
| OUT <- 1/(2 * LAM_A) * (1 - LAM_A * (BET_V + BET_S)) | |
| return(OUT) | |
| } | |
| fun.LAM_A <- function(BET_V, BET_S, GAM_V, SIG_V, SIG_E, SIG_Z) { | |
| NUMER <- (BET_V + BET_S + GAM_V) * SIG_V^2 | |
| DENOM <- (BET_V + BET_S + GAM_V)^2 * SIG_V^2 + BET_V^2 * SIG_E^2 + (1 + BET_S^2/GAM_V^2) * SIG_Z^2 | |
| OUT <- NUMER/DENOM | |
| return(OUT) | |
| } | |
| fun.BET_0 <- function(LAM_A, LAM_0, GAM_0, GAM_V, PSI, OMG, MU_V) { | |
| OUT <- - 1/(2 * LAM_A) * (LAM_0 + LAM_A * (GAM_0 + GAM_V * PSI * OMG * MU_V)) | |
| return(OUT) | |
| } | |
| fun.GAM_0 <- function(LAM_A, LAM_0, BET_0) { | |
| OUT <- - 1/(2 * LAM_A) * (LAM_0 + LAM_A * BET_0) | |
| return(OUT) | |
| } | |
| fun.LAM_0 <- function(LAM_A, BET_V, BET_S, GAM_V, MU_V, BET_0, GAM_0) { | |
| OUT <- (1 - LAM_A * (BET_V + BET_S + GAM_V)) * MU_V - LAM_A * (BET_0 + GAM_0) | |
| return(OUT) | |
| } | |
| ## Define error functions | |
| fun.ERR_SLOPE <- function(ARGS, VAL) { | |
| BET_V <- ARGS[1] | |
| BET_S <- ARGS[2] | |
| GAM_V <- ARGS[3] | |
| LAM_A <- ARGS[4] | |
| SIG_V <- VAL[1] | |
| SIG_E <- VAL[2] | |
| SIG_Z <- VAL[3] | |
| PSI <- fun.PSI(SIG_E, SIG_V) | |
| OMG <- fun.OMG(SIG_Z, GAM_V, PSI, SIG_V) | |
| BET_V_HAT <- fun.BET_V(LAM_A, GAM_V, PSI, OMG) | |
| BET_S_HAT <- fun.BET_S(GAM_V, OMG) | |
| GAM_V_HAT <- fun.GAM_V(LAM_A, BET_V, BET_S) | |
| LAM_A_HAT <- fun.LAM_A(BET_V, BET_S, GAM_V, SIG_V, SIG_E, SIG_Z) | |
| ERROR <- sum(c((BET_V_HAT - BET_V)^2, | |
| (BET_S_HAT - BET_S)^2, | |
| (GAM_V_HAT - GAM_V)^2, | |
| (LAM_A_HAT - LAM_A)^2 | |
| ) | |
| ) | |
| return(ERROR) | |
| } | |
| fun.ERR_LEVEL <- function(ARGS, VAL) { | |
| BET_0 <- ARGS[1] | |
| GAM_0 <- ARGS[2] | |
| SIG_V <- VAL[1] | |
| MU_V <- LAM_0 <- VAL[2] | |
| SIG_E <- VAL[3] | |
| SIG_Z <- VAL[4] | |
| BET_V <- VAL[5] | |
| BET_S <- VAL[6] | |
| GAM_V <- VAL[7] | |
| LAM_A <- VAL[8] | |
| PSI <- fun.PSI(SIG_E, SIG_V) | |
| OMG <- fun.OMG(SIG_Z, GAM_V, PSI, SIG_V) | |
| BET_0_HAT <- fun.BET_0(LAM_A, LAM_0, GAM_0, GAM_V, PSI, OMG, MU_V) | |
| GAM_0_HAT <- fun.GAM_0(LAM_A, LAM_0, BET_0) | |
| LAM_0_HAT <- fun.LAM_0(LAM_A, BET_V, BET_S, GAM_V, MU_V, BET_0, GAM_0) | |
| ERROR <- sum(c((BET_0_HAT - BET_0)^2, | |
| (GAM_0_HAT - GAM_0)^2, | |
| (LAM_0_HAT - LAM_0)^2 | |
| ) | |
| ) | |
| return(ERROR) | |
| } | |
| ## Estimate equilibrium parameters | |
| if (TRUE == FALSE) { | |
| scl.int.SIG_E_LEN <- length(vec.flt.SIG_E) | |
| mat.flt.SLOPE <- foreach(e=1:scl.int.SIG_E_LEN, .combine = "rbind", .inorder = FALSE) %dopar% { | |
| if (e %in% c(1,round(seq(0,scl.int.SIG_E_LEN, length.out = 11)))) { | |
| print(round(e/scl.int.SIG_E_LEN * 100)) | |
| } | |
| vec.flt.VAL <- c(scl.flt.SIG_V, | |
| vec.flt.SIG_E[e], | |
| scl.flt.SIG_Z | |
| ) | |
| vec.flt.PAR <- c(1, -1, 1, 1) | |
| obj.opt.OUT <- optim(par = vec.flt.PAR, | |
| fn = fun.ERR_SLOPE, | |
| method = "BFGS", | |
| VAL = vec.flt.VAL | |
| ) | |
| vec.flt.OUT <- c(vec.flt.SIG_E[e], obj.opt.OUT$par, obj.opt.OUT$value, obj.opt.OUT$convergence) | |
| return(as.numeric(vec.flt.OUT)) | |
| } | |
| mat.df.SLOPE <- as.data.frame(mat.flt.SLOPE) | |
| names(mat.df.SLOPE) <- c("sigE", "betV", "betS", "gamV", "lamA", "err", "fit") | |
| mat.df.SLOPE <- mat.df.SLOPE[(mat.df.SLOPE$fit == 0.00) & (mat.df.SLOPE$err < 0.01) & (mat.df.SLOPE$betS < 0.00), c("sigE", "betV", "betS", "gamV", "lamA")] | |
| scl.int.SIG_E_LEN <- dim(mat.df.SLOPE)[1] | |
| mat.flt.LEVEL <- foreach(e=1:scl.int.SIG_E_LEN, .combine = "rbind", .inorder = FALSE) %dopar% { | |
| if (e %in% c(1,round(seq(0,scl.int.SIG_E_LEN, length.out = 11)))) { | |
| print(round(e/scl.int.SIG_E_LEN * 100)) | |
| } | |
| scl.flt.SIG_E <- mat.df.SLOPE[e, 1] | |
| scl.flt.BET_V <- mat.df.SLOPE[e, 2] | |
| scl.flt.BET_S <- mat.df.SLOPE[e, 3] | |
| scl.flt.GAM_V <- mat.df.SLOPE[e, 4] | |
| scl.flt.LAM_A <- mat.df.SLOPE[e, 5] | |
| vec.flt.VAL <- c(scl.flt.SIG_V, | |
| scl.flt.MU_V, | |
| scl.flt.SIG_E, | |
| scl.flt.SIG_Z, | |
| scl.flt.BET_V, | |
| scl.flt.BET_S, | |
| scl.flt.GAM_V, | |
| scl.flt.LAM_A | |
| ) | |
| vec.flt.PAR <- c(-1, -1) | |
| obj.opt.OUT <- optim(par = vec.flt.PAR, | |
| fn = fun.ERR_LEVEL, | |
| method = "BFGS", | |
| VAL = vec.flt.VAL | |
| ) | |
| vec.flt.OUT <- c(scl.flt.SIG_E, obj.opt.OUT$par, obj.opt.OUT$value, obj.opt.OUT$convergence) | |
| return(as.numeric(vec.flt.OUT)) | |
| } | |
| mat.df.LEVEL <- as.data.frame(mat.flt.LEVEL) | |
| names(mat.df.LEVEL) <- c("sigE", "bet0", "gam0", "err", "fit") | |
| mat.df.LEVEL <- mat.df.LEVEL[(mat.df.LEVEL$fit == 0.00) & (mat.df.LEVEL$err < 0.01), c("sigE", "bet0", "gam0")] | |
| mat.df.LEVEL$lam0 <- scl.flt.MU_V | |
| mat.df.PARAM <- merge(mat.df.SLOPE, mat.df.LEVEL, by = "sigE") | |
| save(mat.df.PARAM, file = paste(scl.str.DAT_DIR, "data--equilibrium-parameters--26aug2015.Rdata")) | |
| } | |
| load(file = paste(scl.str.DAT_DIR, "data--equilibrium-parameters--26aug2015.Rdata")) | |
| ## Plot equilibrium parameters | |
| if (TRUE == FALSE) { | |
| mat.df.PLOT <- mat.df.PARAM[, c("sigE", "betV", "betS", "gamV", "lamA")] | |
| names(mat.df.PLOT) <- c("sigE", | |
| "$\\beta_{\\tilde{v}}$", | |
| "$\\beta_{s}$", | |
| "$\\gamma_{v}$", | |
| "$\\lambda_a$" | |
| ) | |
| mat.df.PLOT$sigE <- mat.df.PLOT$sigE/scl.flt.SIG_V | |
| mat.df.PLOT <- melt(mat.df.PLOT, c("sigE")) | |
| theme_set(theme_bw()) | |
| scl.str.RAW_FILE <- "plot--equilibrium-parameters--26aug2015" | |
| 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() | |
| obj.gg2.PLOT <- obj.gg2.PLOT + scale_colour_brewer(palette="Set1") | |
| obj.gg2.PLOT <- obj.gg2.PLOT + geom_smooth(data = mat.df.PLOT, | |
| aes(x = sigE, | |
| y = value | |
| ), | |
| method = "lm", | |
| formula = y ~ ns(x, 25), | |
| se = FALSE, | |
| size = 1.75 | |
| ) | |
| ## obj.gg2.PLOT <- obj.gg2.PLOT + geom_point(data = mat.df.PLOT, | |
| ## aes(x = sigE, | |
| ## y = value | |
| ## ), | |
| ## size = 0.50, | |
| ## alpha = 0.50 | |
| ## ) | |
| obj.gg2.PLOT <- obj.gg2.PLOT + xlab("$\\sigma_{\\epsilon}/\\sigma_v$") | |
| obj.gg2.PLOT <- obj.gg2.PLOT + ylab("") | |
| obj.gg2.PLOT <- obj.gg2.PLOT + scale_x_continuous(limits = c(0, 1.5), breaks = c(0.00, 0.50, 1.00, 1.50)) | |
| obj.gg2.PLOT <- obj.gg2.PLOT + facet_wrap(~ variable, nrow = 1, scales = "free_y") | |
| obj.gg2.PLOT <- obj.gg2.PLOT + theme(plot.margin = unit(c(1,0.15,0.15,-1), "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 = "none" | |
| ) | |
| obj.gg2.PLOT <- obj.gg2.PLOT + ggtitle("Equilibrium Parameters") | |
| 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 = '')) | |
| } | |
| ## Simulate trader profits | |
| if (TRUE == FALSE) { | |
| vec.flt.BREAKS <- seq(0, 1.5, by = 0.01) | |
| scl.int.BKT_LEN <- length(vec.flt.BREAKS) - 1 | |
| mat.df.SMOOTH <- mat.df.PARAM | |
| mat.df.SMOOTH$bkt <- NA | |
| for (b in 1:scl.int.BKT_LEN) { | |
| vec.int.ROWS <- which((mat.df.SMOOTH$sigE >= vec.flt.BREAKS[b]) & (mat.df.SMOOTH$sigE < vec.flt.BREAKS[b+1])) | |
| if (length(vec.int.ROWS) > 0) { | |
| mat.df.SMOOTH[vec.int.ROWS, ]$bkt <- b | |
| } | |
| } | |
| mat.df.SMOOTH <- ddply(mat.df.SMOOTH, | |
| c("bkt"), | |
| function(X)c(mean(X$sigE, na.rm = TRUE), | |
| mean(X$bet0, na.rm = TRUE), | |
| mean(X$betV, na.rm = TRUE), | |
| mean(X$betS, na.rm = TRUE), | |
| mean(X$gam0, na.rm = TRUE), | |
| mean(X$gamV, na.rm = TRUE), | |
| mean(X$lam0, na.rm = TRUE), | |
| mean(X$lamA, na.rm = TRUE) | |
| ) | |
| ) | |
| names(mat.df.SMOOTH) <- c("bkt", "sigE", "bet0", "betV", "betS", "gam0", "gamV", "lam0", "lamA") | |
| scl.int.BKT_LEN <- length(unique(mat.df.SMOOTH$bkt)) | |
| mat.flt.OUTTER <- foreach(e=1:scl.int.BKT_LEN, .combine = "rbind", .inorder = FALSE) %dopar% { | |
| ## if (e %in% c(1,round(seq(0,scl.int.BKT_LEN, length.out = 11)))) { | |
| ## print(round(e/scl.int.BKT_LEN * 100)) | |
| ## } | |
| print(e) | |
| scl.flt.SIG_E <- mat.df.SMOOTH[e, 2] | |
| scl.flt.BET_0 <- mat.df.SMOOTH[e, 3] | |
| scl.flt.BET_V <- mat.df.SMOOTH[e, 4] | |
| scl.flt.BET_S <- mat.df.SMOOTH[e, 5] | |
| scl.flt.GAM_0 <- mat.df.SMOOTH[e, 6] | |
| scl.flt.GAM_V <- mat.df.SMOOTH[e, 7] | |
| scl.flt.LAM_0 <- mat.df.SMOOTH[e, 8] | |
| scl.flt.LAM_A <- mat.df.SMOOTH[e, 9] | |
| mat.flt.INNER <- foreach(i=1:scl.int.ITER, .combine = "rbind") %do% { | |
| scl.flt.V <- rnorm(1, 0, scl.flt.SIG_V) | |
| vec.flt.EPS <- rnorm(scl.int.T, 0, scl.flt.SIG_E) | |
| vec.flt.S <- rnorm(scl.int.T, scl.flt.V, scl.flt.SIG_Z/scl.flt.GAM_V) | |
| vec.flt.V_TIL <- scl.flt.V + vec.flt.EPS | |
| vec.flt.X <- scl.flt.BET_0 + scl.flt.BET_V * vec.flt.V_TIL + scl.flt.BET_S * vec.flt.S | |
| scl.flt.Y <- scl.flt.GAM_0 + scl.flt.GAM_V * scl.flt.V | |
| vec.flt.Z <- rnorm(scl.int.T, 0, scl.flt.SIG_Z) | |
| vec.flt.A <- vec.flt.X + scl.flt.Y + vec.flt.Z | |
| vec.flt.P <- scl.flt.LAM_0 + scl.flt.LAM_A * vec.flt.A | |
| vec.flt.INNER <- c(mean((vec.flt.V_TIL - vec.flt.P) * vec.flt.X), | |
| mean((scl.flt.V - vec.flt.P) * scl.flt.Y), | |
| mean((scl.flt.V - vec.flt.P) * vec.flt.Z), | |
| mean(-(scl.flt.V - vec.flt.P) * vec.flt.A) | |
| ) | |
| return(vec.flt.INNER) | |
| } | |
| vec.flt.OUTTER <- c(scl.flt.SIG_E, | |
| mean(mat.flt.INNER[, 1]), | |
| mean(mat.flt.INNER[, 2]), | |
| mean(mat.flt.INNER[, 3]), | |
| mean(mat.flt.INNER[, 4]) | |
| ) | |
| return(vec.flt.OUTTER) | |
| } | |
| mat.df.PROFITS <- as.data.frame(mat.flt.OUTTER) | |
| names(mat.df.PROFITS) <- c("sigE", "shortRun", "longRun", "noise", "marketMaker") | |
| save(mat.df.PROFITS, file = paste(scl.str.DAT_DIR, "data--trader-profits--26aug2015.Rdata", sep = "")) | |
| } | |
| load(file = paste(scl.str.DAT_DIR, "data--trader-profits--26aug2015.Rdata", sep = "")) | |
| ## Plot trader profits | |
| if (TRUE == TRUE) { | |
| mat.df.PLOT <- mat.df.PROFITS | |
| names(mat.df.PLOT) <- c("sigE", | |
| "Short-Run Trader", | |
| "Long-Run Trader", | |
| "Noise Trader", | |
| "Market Maker" | |
| ) | |
| mat.df.PLOT$sigE <- mat.df.PLOT$sigE/scl.flt.SIG_V | |
| mat.df.PLOT <- melt(mat.df.PLOT, c("sigE")) | |
| theme_set(theme_bw()) | |
| scl.str.RAW_FILE <- "plot--expected-profits--26aug2015" | |
| 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() | |
| obj.gg2.PLOT <- obj.gg2.PLOT + scale_colour_brewer(palette="Set1") | |
| obj.gg2.PLOT <- obj.gg2.PLOT + geom_smooth(data = mat.df.PLOT, | |
| aes(x = sigE, | |
| y = value | |
| ), | |
| method = "lm", | |
| formula = y ~ ns(x, 6), | |
| se = FALSE, | |
| size = 1.75 | |
| ) | |
| ## obj.gg2.PLOT <- obj.gg2.PLOT + geom_point(data = mat.df.PLOT, | |
| ## aes(x = sigE, | |
| ## y = value | |
| ## ), | |
| ## size = 0.50, | |
| ## alpha = 0.50 | |
| ## ) | |
| obj.gg2.PLOT <- obj.gg2.PLOT + xlab("$\\sigma_{\\epsilon}/\\sigma_v$") | |
| obj.gg2.PLOT <- obj.gg2.PLOT + ylab("") | |
| obj.gg2.PLOT <- obj.gg2.PLOT + scale_x_continuous(limits = c(0, 1.5), breaks = c(0.00, 0.50, 1.00, 1.50)) | |
| obj.gg2.PLOT <- obj.gg2.PLOT + facet_wrap(~ variable, nrow = 1, scales = "free_y") | |
| obj.gg2.PLOT <- obj.gg2.PLOT + theme(plot.margin = unit(c(1,0.15,0.15,-1), "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 = "none" | |
| ) | |
| obj.gg2.PLOT <- obj.gg2.PLOT + ggtitle("Expected Profits") | |
| 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 = '')) | |
| } | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment