Code for Figure 2. Investment-Horizon Spillovers. Alex Chinco and Mao Ye.

motivating-analysis.R

options(width = 80, digits = 6, digits.secs = 6)
rm(list=ls())

library(reshape)

library(foreach)
library(doMC)
registerDoMC(7)

library(ggplot2)
library(ggthemes)
library(grid)
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)

fig_dir <- "insert_directory_name_here"





################################################################################
## how to solve model

## know:    alpha[t] & delta[t]
## guess:   beta[t] & sigma[t]
## compute: beta[t] & sigma[t] -> lambda[t]

## compute: lambda[t] & alpha[t]            -> alpha[t-dt]
## compute: lambda[t] & alpha[t] & delta[t] -> delta[t-dt]
## compute: beta[t] & lambda[t] & sigma[t]  -> sigma[t-dt]
## compute: alpha[t-dt] & sigma[t-dt]       -> lambda[t-dt]
## compute: lambda[t-dt] & sigma[t-dt]      -> beta[t-dt]

## check: (sigma[0] == 1)
## check: (lambda[t] * (1 - alpha[t] * lambda[t]) > 0) for all t

PreviousAlpha <- function(alpha, lambda) {
    return(1 / (4 * lambda * (1 - alpha * lambda)))
}

PreviousDelta <- function(alpha, delta, lambda, dt) {
    return(delta + alpha * lambda^2 * dt)
}

PreviousSigma <- function(beta, lambda, sigma, dt) {
    return(sigma / sqrt(1 - beta * lambda * dt))
}

LambdaSolve <- function(lambda, alpha, sigma, dt) {
    error <- 1 - 2 * (1 - alpha * lambda) * (1 - lambda^2 * dt / sigma^2)
    return(error^2)
}

Beta <- function(lambda, sigma) {
    return(lambda / sigma^2)
}

ComputeParams <- function(sigma_end, num_periods) {
    dt     <- 1/num_periods    
    alpha  <- rep(NA, num_periods)
    delta  <- rep(NA, num_periods)
    sigma  <- rep(NA, num_periods)
    lambda <- rep(NA, num_periods)
    beta   <- rep(NA, num_periods)
    
    alpha[num_periods]  <- 0
    delta[num_periods]  <- 0    
    sigma[num_periods]  <- sigma_end    
    lambda[num_periods] <- optimize(LambdaSolve, c(0, 1000), tol = 1e-8, alpha = 0, sigma = sigma_end, dt = dt)$minimum    
    beta[num_periods]   <- Beta(lambda[num_periods], sigma[num_periods])
    
    for (t in num_periods:2) {
        alpha[t-1]  <- PreviousAlpha(alpha[t], lambda[t])
        delta[t-1]  <- PreviousDelta(alpha[t], delta[t], lambda[t], dt)
        sigma[t-1]  <- PreviousSigma(beta[t], lambda[t], sigma[t], dt)
        lambda[t-1] <- optimize(LambdaSolve, c(0, 1/alpha[t-1]), tol = 1e-8, alpha = alpha[t-1], sigma = sigma[t-1], dt = dt)$minimum
        beta[t-1]   <- Beta(lambda[t-1], sigma[t-1])
    }
    
    return(list("t"      = (2:num_periods - 0.5) / num_periods,
                "alpha"  = alpha[2:num_periods],
                "delta"  = delta[2:num_periods],
                "sigma"  = sigma[2:num_periods],
                "lambda" = lambda[2:num_periods],
                "beta"   = beta[2:num_periods]
                )
           )
}

EquilibriumSolve <- function(sigma_end, num_periods) {
    params      <- ComputeParams(sigma_end, num_periods)
    sigma_start <- params$sigma[1]
    error       <- (sigma_start - 1)
    return(error^2)
}




################################################################################
## marginal benefit
num_iter <- 1e5
freq     <- 2^(1:11)
num_freq <- length(freq) 

avg_profit <- foreach(f = 1:num_freq, .combine = "rbind") %do% {
    num_periods <- freq[f]
    print(num_periods)
    
    set.seed(num_periods)
    sigma_end <- optimize(EquilibriumSolve, c(0, 0.50), tol = 1e-8, num_periods = num_periods)$minimum
    params    <- ComputeParams(sigma_end, num_periods = num_periods)
    
    profit <- foreach(i = 1:num_iter, .combine = "c", .inorder = FALSE) %dopar% {
        set.seed(i * num_periods)
        dt <- 1/num_periods
        v  <- rnorm(1, mean = 0, sd = 1)
        dz <- rnorm(num_periods - 1, mean = 0, sd = sqrt(dt))
        dx <- rep(NA, num_periods - 1)
        x  <- rep(NA, num_periods)
        x[1] <- 0
        
        p    <- rep(NA, num_periods)
        p[1] <- 0
        
        for (t in 1:(num_periods - 1)) {
            dx[t]  <- params$beta[t] * (v - p[t]) * dt
            x[t+1] <- x[t] + dx[t]
            p[t+1] <- p[t] + params$lambda[t] * (dx[t] + dz[t])
        }
        p <- p[2:num_periods]
        
        return(sum((v - p) * dx))
    }
    
    return(data.frame(f = f, value = mean(profit)))
}





################################################################################
## demand rule
num_periods <- 8
sigma_end   <- optimize(EquilibriumSolve, c(0, 0.50), tol = 1e-8, num_periods = num_periods)$minimum
params      <- ComputeParams(sigma_end, num_periods = num_periods)
demand_rule <- data.frame(t = params$t, value = params$beta)





################################################################################
## trading-volume variance
num_periods <- 8
num_iter    <- 1e5

set.seed(num_periods)
sigma_end <- optimize(EquilibriumSolve, c(0, 0.50), tol = 1e-8, num_periods = num_periods)$minimum
params    <- ComputeParams(sigma_end, num_periods = num_periods)

volume <- foreach(i = 1:num_iter, .combine = "cbind") %dopar% {
    set.seed(i * num_periods)
    
    dt <- 1/num_periods
    v  <- rnorm(1, mean = 0, sd = 1)
    dz <- rnorm(num_periods - 1, mean = 0, sd = sqrt(dt))
    dx <- rep(NA, num_periods - 1)
    x  <- rep(NA, num_periods)
    
    p    <- rep(NA, num_periods)
    p[1] <- 0
    
    for (t in 1:(num_periods - 1)) {
        dx[t]  <- params$beta[t] * (v - p[t]) * dt
        p[t+1] <- p[t] + params$lambda[t] * (dx[t] + dz[t])
    }

    return(abs(dx) + abs(dz))
}

volume_var <- data.frame(t = params$t, value = apply(volume, 1, var) * num_periods)





################################################################################
## create plot
obj.vplayout  <- function(x, y) viewport(layout.pos.row = x, layout.pos.col = y)

theme_set(theme_tufte())


name  <- "motivating-analysis"
files <- list("tex" = paste(name,'.tex',sep=''),
              "pdf" = paste(name,'.pdf',sep=''),
              "aux" = paste(name,'.aux',sep=''),
              "log" = paste(name,'.log',sep='')
              )
tikz(file = files$tex, height = 2, width = 7, standAlone=TRUE)
pushViewport(viewport(layout = grid.layout(1, 3)))


title <- "$\\mathrm{E}(\\text{Profit})$ $[\\mathdollar/\\text{sh}]$"
x <- list("limits" = c(0.35, 11),
          "breaks" = c(1, 3, 5, 7, 9, 11),
          "labels" = c("$f = 1 \\phantom{=f}$", "$3$", "$5$", "$7$", "$9$", "$11$")
          )
y <- list("limits" = c(0.5, 1),
          "breaks" = c(0.5, 1),
          "labels" = c("$\\frac{\\sigma}{2}$", "$\\sigma$")
          )
plot <- ggplot(data = avg_profit)
plot <- plot + geom_path(aes(x = f,
                             y = value
                             ),
                         size  = 2,
                         alpha = 0.75
                         )
plot <- plot + ggtitle(title)
plot <- plot + coord_cartesian(xlim = x$limits, ylim = y$limits)
plot <- plot + scale_x_continuous(breaks = x$breaks, labels = x$labels)
plot <- plot + scale_y_continuous(breaks = y$breaks, labels = y$labels)
plot <- plot + theme(plot.margin = unit(c(0.15, 0.50, 0.15, 0.50), "lines"),
                     axis.text   = element_text(size = 12),
                     axis.ticks  = element_blank(),
                     axis.title  = element_blank(),
                     plot.title  = element_text(size = 16)
                     )
print(plot, vp = obj.vplayout(1, 1))


title <- "Demand $[\\text{sh}/\\mathdollar]$"
x <- list("limits" = c(1/8, 1),
          "breaks" = c(0.2, 0.4, 0.6, 0.8, 1),
          "labels" = c("$t = 0.2 \\phantom{=t}$", "$0.4$", "$0.6$", "0.8", "$1.0$")
          )
y <- list("limits" = c(1, 5),
          "breaks" = c(1, 2, 3, 5)
          )
plot <- ggplot(data = demand_rule)
plot <- plot + geom_point(aes(x = t,
                              y = value
                              ),
                          size  = 2,
                          alpha = 0.75
                          )
plot <- plot + ggtitle(title)
plot <- plot + coord_cartesian(xlim = x$limits, ylim = y$limits)
plot <- plot + scale_x_continuous(breaks = x$breaks, labels = x$labels)
plot <- plot + scale_y_continuous(breaks = y$breaks)
plot <- plot + theme(plot.margin = unit(c(0.15, 0.50, 0.15, 0.50), "lines"),
                     axis.text   = element_text(size = 12),
                     axis.ticks  = element_blank(),
                     axis.title  = element_blank(),
                     plot.title  = element_text(size = 16)
                     )
print(plot, vp = obj.vplayout(1, 2))


title <- "$\\mathrm{Var}(\\text{Volume})$ $[\\text{sh}^{\\scriptscriptstyle 2}]$"
x <- list("limits" = c(1/8, 1),
          "breaks" = c(0.2, 0.4, 0.6, 0.8, 1),
          "labels" = c("$t = 0.2 \\phantom{=t}$", "$0.4$", "$0.6$", "0.8", "$1.0$")
          )
y <- list("limits" = c(0.40, 0.71),
          "breaks" = c(0.4, 0.5, 0.7)
          )
plot <- ggplot(data = volume_var)
plot <- plot + geom_point(aes(x = t,
                              y = value
                              ),
                          size  = 2,
                          alpha = 0.75
                          )
plot <- plot + ggtitle(title)
plot <- plot + coord_cartesian(xlim = x$limits, ylim = y$limits)
plot <- plot + scale_x_continuous(breaks = x$breaks, labels = x$labels)
plot <- plot + scale_y_continuous(breaks = y$breaks)
plot <- plot + theme(plot.margin = unit(c(0.15, 0.50, 0.15, 0.50), "lines"),
                     axis.text   = element_text(size = 12),
                     axis.ticks  = element_blank(),
                     axis.title  = element_blank(),
                     plot.title  = element_text(size = 16)
                     )
print(plot, vp = obj.vplayout(1, 3))



dev.off()


system(paste('lualatex', file.path(files$tex)), ignore.stdout = TRUE)
system(paste('mv ', files$pdf, ' ', fig_dir, sep = ''))
system(paste('rm ', files$tex, sep = ''))
system(paste('rm ', files$aux, sep = ''))
system(paste('rm ', files$log, sep = ''))