Code to solve for the equilibrium parameters in the Kyle (1989) model

plot--kyle-1989-model-solution--12nov2015.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{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(7)



## 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.flt.RHO <- 1

scl.int.N <- 2
scl.int.M <- 1

scl.flt.TAU_E <- scl.int.N
scl.flt.TAU_V <- 1
vec.flt.TAU_Z <- seq(0.01, 2, by = 0.01)










## Define equlibrium functions
fun.VARPHI_I <- function(N, GAMMA_I, TAU_Z, TAU_E) {
    OUT <- (N - 1) * GAMMA_I^2 / ((N - 1) * GAMMA_I^2 + TAU_E/TAU_Z)
    return(OUT)
}

fun.VARPHI_U <- function(N, GAMMA_I, TAU_Z, TAU_E) {
    OUT <- N * GAMMA_I^2 / (N * GAMMA_I^2 + TAU_E/TAU_Z)
    return(OUT)
}

fun.TAU_I <- function(TAU_V, TAU_E, VARPHI_I, N) {
    OUT <- TAU_V + TAU_E + VARPHI_I * (N - 1) * TAU_E
    return(OUT)
}

fun.TAU_U <- function(TAU_V, TAU_E, VARPHI_U, N) {
    OUT <- TAU_V + VARPHI_U * N * TAU_E
    return(OUT)
}




## Define error functions 
fun.OBJECTIVE <- function(ARGS, VAL) {
    GAMMA_I <- ARGS[1]
    ZETA    <- ARGS[2]
    XI_I    <- ARGS[3]
    XI_U    <- ARGS[4]    
    
    RHO    <- VAL[1]
    TAU_V  <- VAL[2]
    TAU_Z  <- VAL[3]
    TAU_E  <- VAL[4]
    N      <- VAL[5]
    M      <- VAL[6]

    VARPHI_I <- fun.VARPHI_I(N, GAMMA_I, TAU_Z, TAU_E)
    VARPHI_U <- fun.VARPHI_U(N, GAMMA_I, TAU_Z, TAU_E)

    TAU_I <- fun.TAU_I(TAU_V, TAU_E, VARPHI_I, N)
    TAU_U <- fun.TAU_U(TAU_V, TAU_E, VARPHI_U, N)

    EQN1_ERR <- 1 - (N * XI_I + M * XI_U)
    EQN2_ERR <- (1 - ZETA) - ((1 - VARPHI_I) * (1 - XI_I))
    EQN3_ERR <- (RHO * GAMMA_I/TAU_E) - ((1 - VARPHI_I) * (1 - 2 * ZETA) / (1 - ZETA))
    EQN4_ERR <- ((XI_U * ZETA * TAU_U/(1 - XI_U)) + (XI_U * RHO * GAMMA_I * TAU_I)/TAU_E) - (ZETA * TAU_U - VARPHI_U * TAU_I)

    ERROR <- sum(c(EQN1_ERR^2, EQN2_ERR^2, EQN3_ERR^2, EQN4_ERR^2))

    if (XI_I < VARPHI_U/N) {
        ERROR <- ERROR + 100
    }
    if (XI_U > (1 - VARPHI_U)/M) {
        ERROR <- ERROR + 100
    }
    LAM <- ZETA * TAU_E / (TAU_I * GAMMA_I)
    if ((VARPHI_U * TAU_E) / (GAMMA_I * TAU_U * LAM) > 1) {
        ERROR <- ERROR + 100
    }
    
    
    return(ERROR)
}






## Estimate equilibrium parameters
if (TRUE == TRUE) {
    scl.int.TAU_Z_LEN <- length(vec.flt.TAU_Z)
    mat.flt.COEF      <- foreach(z=1:scl.int.TAU_Z_LEN, .combine = "rbind", .inorder = FALSE) %dopar% {
        print(z)
        
        vec.flt.VAL <- c(scl.flt.RHO,
                         scl.flt.TAU_V,
                         vec.flt.TAU_Z[z],
                         scl.flt.TAU_E,
                         scl.int.N,
                         scl.int.M
                         )
        vec.flt.PAR <- c(1, 1/4, 1/6, 1/3)        
        obj.opt.OUT <- optim(par     = vec.flt.PAR,
                             fn      = fun.OBJECTIVE,
                             method  = "L-BFGS-B",
                             lower   = c(0, 0, 0, 0),
                             upper   = c(Inf, 1/2, 1/scl.int.N, Inf),
                             VAL     = vec.flt.VAL
                             )
        vec.flt.OUT <- c(vec.flt.TAU_Z[z], obj.opt.OUT$par, obj.opt.OUT$value, obj.opt.OUT$convergence)
        return(as.numeric(vec.flt.OUT))
    }
    mat.df.COEF        <- as.data.frame(mat.flt.COEF)
    names(mat.df.COEF) <- c("tauZ", "gamI", "zet", "xiI", "xiU", "err", "fit")        

    mat.df.COEF$varphiI <- fun.VARPHI_I(scl.int.N, mat.df.COEF$gamI, mat.df.COEF$tauZ, scl.flt.TAU_E)
    mat.df.COEF$varphiU <- fun.VARPHI_U(scl.int.N, mat.df.COEF$gamI, mat.df.COEF$tauZ, scl.flt.TAU_E)

    mat.df.COEF$tauI <- fun.TAU_I(scl.flt.TAU_V, scl.flt.TAU_E, mat.df.COEF$varphiI, scl.int.N)
    mat.df.COEF$tauU <- fun.TAU_U(scl.flt.TAU_V, scl.flt.TAU_E, mat.df.COEF$varphiU, scl.int.N)

    mat.df.COEF$lam <- mat.df.COEF$zet * scl.flt.TAU_E / (mat.df.COEF$tauI * mat.df.COEF$gamI)

    mat.df.COEF$betI <- mat.df.COEF$xiI/mat.df.COEF$lam
    mat.df.COEF$betU <- mat.df.COEF$xiU/mat.df.COEF$lam

    mat.df.COEF$lamI <- 1/((scl.int.N - 1) * mat.df.COEF$betI + scl.int.M * mat.df.COEF$betU)
    mat.df.COEF$lamU <- 1/(scl.int.N * mat.df.COEF$betI + (scl.int.M - 1) * mat.df.COEF$betU)

    mat.df.COEF$muI <- 0
    mat.df.COEF$muU <- 0

    mat.df.COEF$tet <- (mat.df.COEF$varphiU * scl.flt.TAU_E) / (mat.df.COEF$gamI * mat.df.COEF$tauU * mat.df.COEF$lam)
    
    save(mat.df.COEF, file = paste(scl.str.DAT_DIR, "data--kyle-1989-coef--12nov2015.Rdata"))
}
load(file = paste(scl.str.DAT_DIR, "data--kyle-1989-coef--12nov2015.Rdata"))




## Plot endogenous parameters
if (TRUE == TRUE) {
    mat.df.PLOT        <- mat.df.COEF[, c("tauZ", "gamI", "zet", "xiI", "xiU")]
    names(mat.df.PLOT) <- c("tauZ",
                            "$\\gamma_I$",
                            "$\\zeta$",
                            "$\\xi_I$",
                            "$\\xi_U$"
                            )
    mat.df.PLOT        <- melt(mat.df.PLOT, c("tauZ"))
    
    theme_set(theme_bw())
    
    scl.str.RAW_FILE <- "plot--kyle-1989-model-solution--endogenous-param--12nov2015"
    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_path(data = mat.df.PLOT,
                                             aes(x = tauZ,
                                                 y = value
                                                 ),
                                             size  = 1.25,
                                             alpha = 0.50
                                             )
    obj.gg2.PLOT <- obj.gg2.PLOT + xlab("$\\tau_z$")
    obj.gg2.PLOT <- obj.gg2.PLOT + ylab("")
    obj.gg2.PLOT <- obj.gg2.PLOT + scale_x_continuous(limits = c(0, 2.0), breaks = c(0.00, 0.50, 1.00, 1.50, 2.00))
    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("Endogenous 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 = ''))
}





## Plot price-impact coefficients
if (TRUE == TRUE) {
    mat.df.PLOT        <- mat.df.COEF[, c("tauZ", "lam", "lamI", "lamU", "tet")]
    names(mat.df.PLOT) <- c("tauZ",
                            "$\\lambda$",
                            "$\\lambda_I$",
                            "$\\lambda_U$",
                            "$\\theta$"
                            )
    mat.df.PLOT        <- melt(mat.df.PLOT, c("tauZ"))
    
    theme_set(theme_bw())
    
    scl.str.RAW_FILE <- "plot--kyle-1989-model-solution--price-impact-coef--12nov2015"
    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_path(data = mat.df.PLOT,
                                             aes(x = tauZ,
                                                 y = value
                                                 ),
                                             size  = 1.25,
                                             alpha = 0.50
                                             )
    obj.gg2.PLOT <- obj.gg2.PLOT + xlab("$\\tau_z$")
    obj.gg2.PLOT <- obj.gg2.PLOT + ylab("")
    obj.gg2.PLOT <- obj.gg2.PLOT + scale_x_continuous(limits = c(0, 2.0), breaks = c(0.00, 0.50, 1.00, 1.50, 2.00))
    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("Price-Impact Coefficients")    
    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 price-impact coefficients
if (TRUE == TRUE) {
    mat.df.PLOT        <- mat.df.COEF[, c("tauZ", "tauI", "tauU", "varphiI", "varphiU")]
    names(mat.df.PLOT) <- c("tauZ",
                            "$\\tau_I$",
                            "$\\tau_U$",
                            "$\\varphi_I$",
                            "$\\varphi_U$"
                            )
    mat.df.PLOT        <- melt(mat.df.PLOT, c("tauZ"))
    
    theme_set(theme_bw())
    
    scl.str.RAW_FILE <- "plot--kyle-1989-model-solution--information-efficiency-parameters--12nov2015"
    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_path(data = mat.df.PLOT,
                                             aes(x = tauZ,
                                                 y = value
                                                 ),
                                             size  = 1.25,
                                             alpha = 0.50
                                             )
    obj.gg2.PLOT <- obj.gg2.PLOT + xlab("$\\tau_z$")
    obj.gg2.PLOT <- obj.gg2.PLOT + ylab("")
    obj.gg2.PLOT <- obj.gg2.PLOT + scale_x_continuous(limits = c(0, 2.0), breaks = c(0.00, 0.50, 1.00, 1.50, 2.00))
    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("Information-Efficiency 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 = ''))
}