DM Monte Carlo.Rmd

---
title: "Monte Carlo"
author: "Nathan Hart"
date: "6/8/2021"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

(.packages())

#install required packages
requiredPackages = c("rmarkdown", "knitr", "FinCal", "ggplot2", "dplyr", "reshape2")
for(p in requiredPackages){
  if(!require(p,character.only = TRUE)) install.packages(p, repos = "http://cran.us.r-project.org")
  library(p,character.only = TRUE)
}
rm(requiredPackages)

```

## R Markdown

### Summary:
Three options for dealing with SO2 pollution:
1. no scrubber
2. wet scrubber
3. dry scrubber

dry scrubber can be turned on and off dynamically at no cost. A threshold must be determined for switching on/off.

### Parameters

Writing down all of the constants or estimates in the case here so they are easy to reference and adjust as needed.

```{r parameters}
#financial params
quarter.discount.rate <- 0.02411
wet.initial.investment.percent <- 0.75

#pricing params
ea.price.per.ton <- 1150 #to play around with for direct calculation
ea.initial.ppt <- 800 #initial value. will simulate with random walk, expected to increase
drift <- 0.0125
sd.error <- 0.06
num.quarters <- 60 #15 years of quarterly price changes

#emissions data
emissions.per.quarter <- 4000
EA.allocation <- 2000
discount.rate <- 0.1

#wet scrubber params
#estimated value, but up to 99% possible
wet.scrub.efficiency <- 0.98 
#cost per ton of SO2 processed
wet.cost.per.ton <- 10 
#120M to install system
wet.investment <- 120000000 
#unit is fiscal quarters, 3-5 but 10% chance of less than one year. 20% chance of more than a year
wet.install.time.quarters <- 4 

#dry scrubber params
dry.investment <- 5000000
dry.install.time.quarters <- 1 #estimate
dry.var.cost.per.ton <- 600
dry.scrub.efficiency <- 0.55 #40%-70% is range, choosing 55 as described by prof
dry.on.threshold <- 600

#import csv coal data
coal.data <- read.csv(file.choose(), header=TRUE, sep=",", stringsAsFactors = TRUE)

```

### Direct Calculations

Not doing any simulation here, just calculating using average or estimated values.

From case: "I can compare the scrubbers at $1150 EA cost"


```{r direct.calculation}
#no scrubs
total.cost.none <- (emissions.per.quarter - EA.allocation) * ea.price.per.ton
print(total.cost.none)

#wet scrubber
wet.emissions <- emissions.per.quarter * (1-wet.scrub.efficiency)
total.cost.wet <- (wet.emissions - EA.allocation) * ea.price.per.ton #negative numbers are earnings
print(total.cost.wet)

#dry scrubber
dry.emissions <- emissions.per.quarter * (1 - dry.scrub.efficiency)
operating.cost.dry <- emissions.per.quarter * dry.var.cost.per.ton
total.cost.dry <- (dry.emissions - EA.allocation) * ea.price.per.ton + operating.cost.dry
print(total.cost.dry)

```

### Simulation setup

Monte Carlo simulations for uncertain variables

Random walk simulation for prices
https://nwfsc-timeseries.github.io/atsa-labs/sec-tslab-random-walks-rw.html

Vectors of payoffs, npv function from FinCal library to calculate NPV

functions for calculating emissions and npv of each solution. 
factoring in on/off decisions, build time, build costs and pay timings

```{r simulation setup}
#parameter number of times to simulate
num.sims <- 1000

#random walk function for EA prices
random.walk.ea <- function(seed){
  ## set random number seed
  set.seed(seed)
  ## initialize arrays. unclear if SD error needs to be multiplied by initial value or not
  prices.sim <- noise <- rnorm(n = num.quarters, mean = drift * ea.initial.ppt, sd = sd.error * ea.initial.ppt)
  ## compute values 2 thru TT
  prices.sim[1] <- ea.initial.ppt
  for (t in 2:num.quarters) {
    #look at last value, add a deviation (gaussian noise)
    #could also add drift here multiplied by most recent value if we are looking for more an exponential drift?
      prices.sim[t] <- (prices.sim[t - 1] + noise[t])
  }
  return(prices.sim)
}

#no scrubber npv calculation function
no.scrub.npv <- function(simulated.prices){
  no.scrub.payoffs <- (emissions.per.quarter - EA.allocation) * simulated.prices
  no.scrub.cost <- npv(quarter.discount.rate, no.scrub.payoffs)
  return(no.scrub.cost)
}

#wet scrubber npv calc function
wet.scrub.npv <- function(simulated.prices){
  wet.emissions <- emissions.per.quarter * (1-wet.scrub.efficiency)
  wet.payoffs <- (wet.emissions - EA.allocation) * simulated.prices
  #don't get emissions benefits until operational, no scrub emissions until then
  wet.payoffs[1:wet.install.time.quarters] <- (emissions.per.quarter - EA.allocation) * simulated.prices[1:wet.install.time.quarters]
  #TODO may be off by one below, might need a initial cost before first quarter. will be close either way
  wet.payoffs[1] <- wet.payoffs[1] + wet.initial.investment.percent * wet.investment
  wet.payoffs[wet.install.time.quarters] <- wet.payoffs[wet.install.time.quarters] + (1-wet.initial.investment.percent) * wet.investment
  wet.total.cost <- npv(quarter.discount.rate, wet.payoffs) #negative numbers are earnings
  return(wet.total.cost)
}

wet.scrub.npv(random.walk.ea(123))

#dry scrubber npv calculation function
dry.scrub.npv <- function(simulated.prices, threshold.ea.price){
  #excess emissions when turned on
  dry.emissions <- (emissions.per.quarter - EA.allocation) * (1 - dry.scrub.efficiency)
  off.emissions <- emissions.per.quarter - EA.allocation
  #vector with emissions based on on/off threshold
  real.excess.emissions <- ifelse(simulated.prices > threshold.ea.price, dry.emissions, off.emissions)
  #print(real.excess.emissions)
  #vector with costs per quarter factoring in on/off decision
  dry.payoffs <- real.excess.emissions * simulated.prices
  #don't get emissions benefits until operational, no scrub emissions until then
  dry.payoffs[1:dry.install.time.quarters] <- off.emissions * simulated.prices[1:dry.install.time.quarters]
  #initial investment to build system
  dry.payoffs[1] <- dry.payoffs[1] + dry.investment
  #variable cost, not counting time under construction
  dry.payoffs[dry.install.time.quarters:num.quarters] <- ifelse(simulated.prices[dry.install.time.quarters:num.quarters] > threshold.ea.price, dry.payoffs[dry.install.time.quarters:num.quarters] + dry.var.cost.per.ton * real.excess.emissions[dry.install.time.quarters:num.quarters], dry.payoffs[dry.install.time.quarters:num.quarters])
  #npv of dry scrubber project
  dry.total.cost <- npv(quarter.discount.rate, dry.payoffs)
  #print(dry.payoffs)
  return(dry.total.cost)
  #TODO: factor in different coal types, terminal value for npv?
}

dry.scrub.npv(random.walk.ea(123), dry.on.threshold)

#calculate total emissions from dry scrub solution
dry.scrub.emissions <- function(simulated.prices){
  
}


sim.data <- data.frame(
  id = c(1:num.sims),
  no.scrubber.npv = c(0),
  wet.scrubber.npv = c(0),
  dry.scrubber.npv = c(0),
  dry.scrubber.emissions = c(0)
  
  )

#show data table for next section, plus output from random walk algorithm
head(sim.data)
sample.walk <- random.walk.ea(22)
plot(sample.walk)
#wet.scrub.npv(sample.walk)
```

### Simulation

Simulation iterations set above, will create many different EA price profiles to get distribution of outcomes

Calculates the NPV and emissions of each option for each potential EA curve.

Emissions for no scrubber and wet scrubber are constant

For dry scrubber: type of coal will have impact on efficiency. Emissions will vary based on price and threshold chosen

Note: positive numbers are costs, as that is more or less the "default" scenario

Format output to make profit/loss clear

```{r simulation}
for (i in 1:num.sims) {
  #use index as new seed to get many possible EA price curves
  prices.sim <- random.walk.ea(i)
  
  #no scrubber
  sim.data$no.scrubber.npv[i] <- no.scrub.npv(prices.sim)
  
  #wet scrubber
  sim.data$wet.scrubber.npv[i] <- wet.scrub.npv(prices.sim)
  
  #dry scrubber
  sim.data$dry.scrubber.npv[i] <- dry.scrub.npv(prices.sim, dry.on.threshold)
  
  #emissions will vary now in addition to costs/savings

}

#no scrub emissions, assuming emissions are constant, implied in case
no.scrub.emissions <- (emissions.per.quarter - EA.allocation) * num.quarters
#print(no.scrub.emissions)

#wet emissions assuming constant emissions
total.emissions.wet <- (wet.emissions - EA.allocation) * num.quarters
#print(total.emissions.wet)

head(sim.data)

plot.data <- melt(sim.data[,2:4])
#View(plot.data)
ggplot(plot.data, aes(x=value, fill = variable), alpha = 0.5) + geom_histogram(position = 'dodge', color = "black")
  
```

### Optimization

Loop through different values of threshold for the dry scrubber to determine optimal threshold.

Ideas: make a 95% confidence interval for EA prices and use those curves



```{r optimization}
thresholds <- seq(500, 1000, 10)
efficiencies <- seq(40,60,5)
opt.sim.count <- 50

optimization.data <- data.frame(threshold.values = thresholds, npv = c(0), efficiency.values = c(0))
eff.optimization.data <- data.frame()

hist(coal.data$Percent.SO2.reduction)
#TODO factor in prediction error: make an error vector as an argument to npv function


#TODO: make efficiency a parameter to dry scrub npv function
for(dry.scrub.efficiency in efficiencies){
  for(threshold in thresholds){
    npv.vec <- vector(mode = "numeric", length = opt.sim.count)
    for(i in 0:(opt.sim.count-1)){
      prices.sim <- random.walk.ea(i)
      #simulate 60 quarters
      npv.vec[i] <- dry.scrub.npv(prices.sim, threshold)
    }
    #average of each simulation
    optimization.data$npv[optimization.data$threshold.values == threshold] <- mean(npv.vec)
    #organized by efficiency level
    optimization.data$efficiency.values[optimization.data$threshold.values == threshold] <- dry.scrub.efficiency
  }
  #append new data for different efficiency level
  eff.optimization.data <- rbind(eff.optimization.data, optimization.data)
}


ggplot(eff.optimization.data, aes(x=threshold.values, y=npv, color=efficiency.values)) + geom_point()

#View(optimization.data)


```