Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
Code: Risk Aversion, Information Choice, and Price Impact
## 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(foreach)
scl.str.DAT_DIR <- "~/Dropbox/research/index_trading_and_noise/code/"
scl.str.FIG_DIR <- "~/Dropbox/research/index_trading_and_noise/figures/"
## Global params
scl.flt.SIGV <- 1
scl.flt.MUS <- 10
scl.flt.SIGS <- 1
scl.int.N <- 10
scl.flt.PHI <- 1
scl.flt.LAM <- 1
## Coefficient values
fun.ALP0 <- function(GAM0, GAM1, I, PHI, MUS, SIGS) {
NUMER <- GAM0 - GAM1 * MUS
DENOM <- (I + 1) + PHI * GAM1 * SIGS^2
OUT <- - (NUMER / DENOM) * (1/GAM1)
return(OUT)
}
fun.ALP1 <- function(GAM1, I, PHI, SIGS) {
NUMER <- 1
DENOM <- (I + 1) + PHI * GAM1 * SIGS^2
OUT <- (NUMER/DENOM) * (1/GAM1)
return(OUT)
}
fun.BET0 <- function(ALP0, ALP1, I, E, PHI, SIGV, MUS, SIGS) {
NUMER <- (I * ALP0) * (PHI/E) * SIGS^2 + (I * ALP1) * MUS
DENOM <- SIGS^2 + I^2 * ALP1^2 * SIGV^2
OUT <- (NUMER/DENOM) * SIGV^2
return(OUT)
}
fun.BET1 <- function(ALP1, I, E, PHI, SIGV, SIGS) {
NUMER <- (PHI/E) * SIGS^2 + (I * ALP1)
DENOM <- SIGS^2 + I^2 * ALP1^2 * SIGV^2
OUT <- (I * ALP1) * (NUMER/DENOM) * SIGV^2
return(OUT)
}
fun.BET2 <- function(ALP1, BET1, I) {
OUT <- BET1/(I * ALP1)
return(OUT)
}
fun.GAM0 <- function(BET0, BET1, BET2, SIGV, MUS, SIGS) {
KAP <- (BET1^2 * SIGV^2)/(BET2^2 * SIGS^2 + BET1^2 * SIGV^2)
OUT <- (KAP/BET1) * (BET0 - BET2 * MUS) / (KAP/BET1 - 1)
return(OUT)
}
fun.GAM1 <- function(BET1, BET2, E, PHI, SIGV, SIGS) {
KAP <- (BET1^2 * SIGV^2)/(BET2^2 * SIGS^2 + BET1^2 * SIGV^2)
OUT <- (PHI/E) * (KAP - 1) * SIGV^2 / (KAP/BET1 - 1)
return(OUT)
}
## Utility constants
fun.Ai <- function(ALP1, GAM1, I, PHI) {
A11 <- 1 - GAM1 * I * ALP1
A12 <- (1/2) * GAM1
A21 <- A12
A22 <- 0
OUT <- - PHI * ALP1 * matrix(c(A11, A12, A21, A22), ncol = 2)
return(OUT)
}
fun.Bi <- function(ALP0, ALP1, GAM0, GAM1, I, PHI) {
B1 <- ALP0 - GAM0 * ALP1 - 2 * GAM1 * I * ALP0 * ALP1
B2 <- GAM1 * ALP0
OUT <- - PHI * matrix(c(B1, B2), ncol = 1)
return(OUT)
}
fun.Ci <- function(ALP0, GAM0, GAM1, I, PHI, LAM) {
OUT <- PHI * (LAM + (GAM0 + GAM1 * I * ALP0) * ALP0)
return(OUT)
}
fun.Au <- function(ALP1, GAM1, I, U, PHI) {
A11 <- - (1 - GAM1 * I * ALP1) * (I * ALP1)
A12 <- 1/2 * (1 - 2 * GAM1 * I * ALP1)
A21 <- A12
A22 <- GAM1
OUT <- - PHI/U * matrix(c(A11, A12, A21, A22), ncol = 2)
return(OUT)
}
fun.Bu <- function(ALP0, ALP1, GAM0, GAM1, I, U, PHI) {
B1 <- GAM0 * I * ALP1 - I * ALP0 + 2 * GAM1 * I^2 * ALP1 * ALP0
B2 <- - GAM0 - 2 * GAM1 * I * ALP0
OUT <- - PHI/U * matrix(c(B1, B2), ncol = 1)
return(OUT)
}
fun.Cu <- function(ALP0, GAM0, GAM1, I, U, PHI) {
OUT <- - (PHI/U) * (I * ALP0) * (GAM0 + GAM1 * I * ALP0)
return(OUT)
}
## Expected utility
fun.UTIL <- function(SIGV, SIGS, A, B, C) {
SIG <- diag(c(SIGV^2, SIGS^2))
DENOM <- det(diag(c(1,1)) - 2 * SIG %*% A)^(1/2)
NUMER <- exp(1/2 * t(B) %*% solve(diag(c(1,1)) - 2 * SIG %*% A) %*% SIG %*% B + C)
OUT <- - NUMER/DENOM
return(OUT)
}
## Optimization functions
fun.OPT_ALP1 <- function(x, PARAM) {
I <- PARAM[1]
U <- PARAM[2]
PHI <- PARAM[3]
SIGV <- PARAM[4]
MUS <- PARAM[5]
SIGS <- PARAM[6]
BET1 <- fun.BET1(x, I, U, PHI, SIGV, SIGS)
BET2 <- fun.BET2(x, BET1, I)
GAM1 <- fun.GAM1(BET1, BET2, U, PHI, SIGV, SIGS)
ALP1 <- fun.ALP1(GAM1, I, PHI, SIGS)
OUT <- (x - ALP1)^2
return(OUT)
}
fun.OPT_ALP0 <- function(x, PARAM) {
I <- PARAM[1]
U <- PARAM[2]
PHI <- PARAM[3]
SIGV <- PARAM[4]
MUS <- PARAM[5]
SIGS <- PARAM[6]
ALP1 <- PARAM[7]
BET1 <- PARAM[8]
BET2 <- PARAM[9]
GAM1 <- PARAM[10]
BET0 <- fun.BET0(x, ALP1, I, U, PHI, SIGV, MUS, SIGS)
GAM0 <- fun.GAM0(BET0, BET1, BET2, SIGV, MUS, SIGS)
ALP0 <- fun.ALP0(GAM0, GAM1, I, PHI, MUS, SIGS)
OUT <- (x - ALP0)^2
return(OUT)
}
fun.COEF <- function(PARAM) {
I <- PARAM[1]
U <- PARAM[2]
PHI <- PARAM[3]
SIGV <- PARAM[4]
MUS <- PARAM[5]
SIGS <- PARAM[6]
ALP1 <- as.numeric(optimize(fun.OPT_ALP1,
c(0, 1),
PARAM = c(I, U, PHI, SIGV, MUS, SIGS)
)[1]
)
BET1 <- fun.BET1(ALP1, I, U, PHI, SIGV, SIGS)
BET2 <- fun.BET2(ALP1, BET1, I)
GAM1 <- fun.GAM1(BET1, BET2, U, PHI, SIGV, SIGS)
ALP0 <- as.numeric(optimize(fun.OPT_ALP0,
c(0, 1),
PARAM = c(I, U, PHI, SIGV, MUS, SIGS, ALP1, BET1, BET2, GAM1)
)[1]
)
BET0 <- fun.BET0(ALP0, ALP1, I, U, PHI, SIGV, MUS, SIGS)
GAM0 <- fun.GAM0(BET0, BET1, BET2, SIGV, MUS, SIGS)
return(c(ALP0, BET0, GAM0, ALP1, BET1, GAM1))
}
fun.OPT_INFO <- function(x, PARAM) {
N <- PARAM[1]
U <- N - x
PHI <- PARAM[2]
SIGV <- PARAM[3]
MUS <- PARAM[4]
SIGS <- PARAM[5]
LAM <- PARAM[6]
COEF <- fun.COEF(c(x, U, PHI, SIGV, MUS, SIGS))
ALP0 <- COEF[1]
BET0 <- COEF[2]
GAM0 <- COEF[3]
ALP1 <- COEF[4]
BET1 <- COEF[5]
GAM1 <- COEF[6]
Ai <- fun.Ai(ALP1, GAM1, x, PHI)
Bi <- fun.Bi(ALP0, ALP1, GAM0, GAM1, x, PHI)
Ci <- fun.Ci(ALP0, GAM0, GAM1, x, PHI, LAM)
Au <- fun.Au(ALP1, GAM1, x, U, PHI)
Bu <- fun.Bu(ALP0, ALP1, GAM0, GAM1, x, U, PHI)
Cu <- fun.Cu(ALP0, GAM0, GAM1, x, U, PHI)
Ui <- fun.UTIL(SIGV, SIGS, Ai, Bi, Ci)
Uu <- fun.UTIL(SIGV, SIGS, Au, Bu, Cu)
OUT <- (Ui - Uu)^2
return(OUT)
}
## Create plot data
fun.WRAPPER <- function(PARAM) {
N <- PARAM[1]
PHI <- PARAM[2]
SIGV <- PARAM[3]
MUS <- PARAM[4]
SIGS <- PARAM[5]
LAM <- PARAM[6]
I <- as.numeric(optimize(fun.OPT_INFO,
c(0, N),
PARAM = c(N, PHI, SIGV, MUS, SIGS, LAM)
)[1]
)
U <- N - I
COEF <- fun.COEF(c(I, U, PHI, SIGV, MUS, SIGS))
ALP0 <- COEF[1]
BET0 <- COEF[2]
GAM0 <- COEF[3]
ALP1 <- COEF[4]
BET1 <- COEF[5]
GAM1 <- COEF[6]
Ai <- fun.Ai(ALP1, GAM1, I, PHI)
Bi <- fun.Bi(ALP0, ALP1, GAM0, GAM1, I, PHI)
Ci <- fun.Ci(ALP0, GAM0, GAM1, I, PHI, LAM)
Ui <- fun.UTIL(SIGV, SIGS, Ai, Bi, Ci)
KAP <- (BET1^2 * SIGV^2)/(GAM1^2 * SIGS^2 + BET1^2 * SIGV^2)
EP <- BET0 - GAM1 * MUS
return(c(COEF, KAP, EP, I/N, Ui))
}
vec.flt.SIGS <- seq(0.50, 2, by = 0.01)
scl.flt.SIGS_LEN <- length(vec.flt.SIGS)
mat.flt.DATA <- foreach(s=1:scl.flt.SIGS_LEN, .combine = "rbind") %do% {
vec.flt.PARAM <- c(scl.int.N,
scl.flt.PHI,
scl.flt.SIGV,
scl.flt.MUS,
vec.flt.SIGS[s],
scl.flt.LAM
)
vec.flt.VALUES <- fun.WRAPPER(vec.flt.PARAM)
return(c(vec.flt.SIGS[s], vec.flt.VALUES))
}
mat.df.DATA <- data.frame(mat.flt.DATA)
names(mat.df.DATA) <- c("sigs",
"$\\alpha_0$ $[\\sfrac{\\text{sh}}{\\text{tr}}]$",
"$\\beta_0$ $[\\sfrac{\\mathdollar}{\\text{sh}}]$",
"$\\gamma_0$ $[\\sfrac{\\mathdollar}{\\text{sh}}]$",
"$\\alpha_1$ $[\\sfrac{\\text{sh}^2}{(\\mathdollar \\cdot \\text{tr})}]$",
"$\\beta_1$ $[1]$",
"$\\beta_2 = \\gamma_1$ $[\\sfrac{\\mathdollar}{\\text{sh}^2}]$",
"$\\kappa$ $[1]$",
"$\\mathrm{E}(p)$ $[\\sfrac{\\mathdollar}{\\text{sh}}]$",
"$\\sfrac{\\mathcal{I}}{\\mathcal{N}}$ $[1]$",
"$\\mathrm{E}(\\text{Utility})$ $[1]$"
)
## Plot results
mat.df.PLOT <- melt(mat.df.DATA[, 1:7], c("sigs"))
theme_set(theme_bw())
scl.str.RAW_FILE <- 'plot--parameter-values'
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.50, 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_point(aes(x = sigs,
y = value
),
size = 0.75
)
obj.gg2.PLOT <- obj.gg2.PLOT + xlab("$\\sigma_s$ $[\\text{sh}]$")
obj.gg2.PLOT <- obj.gg2.PLOT + ylab("")
obj.gg2.PLOT <- obj.gg2.PLOT + facet_wrap(~ variable, nrow = 2, scales = "free_y")
obj.gg2.PLOT <- obj.gg2.PLOT + theme(plot.margin = unit(c(1,0.15,0.15,-0.85), "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 = "none"
)
obj.gg2.PLOT <- obj.gg2.PLOT + ggtitle("Parameter Values")
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 = ''))
## Plot results
mat.df.PLOT <- melt(mat.df.DATA[, c(1,8,9,10,11)], c("sigs"))
theme_set(theme_bw())
scl.str.RAW_FILE <- 'plot--summary-statistics'
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_point(aes(x = sigs,
y = value
),
size = 0.75
)
obj.gg2.PLOT <- obj.gg2.PLOT + xlab("$\\sigma_s$ $[\\text{sh}]$")
obj.gg2.PLOT <- obj.gg2.PLOT + ylab("")
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,-0.85), "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 = "none"
)
obj.gg2.PLOT <- obj.gg2.PLOT + ggtitle("Summary Statistics")
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