Code: Risk Aversion, Information Choice, and Price Impact

plot--equilibrium-coefficients--25jun2015.R

## 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 = ''))