flaringABM_validation.Rmd

---
title: "FlaringABM validation"
author: "Nick Willems"
---

## Distributions of operator production volumes

```{r}
market_shares <- fread("./inputs/processed/firm_market_shares.csv")

lookup_sql <- "SELECT %s FROM string_lookup WHERE column_name='%s' AND %s"

sql <- "SELECT oil_prod, gas_prod, model, RunID, time, firmID FROM agent_states 
        WHERE time IN (SELECT MIN(time) FROM agent_states UNION SELECT 0 UNION SELECT MAX(time) FROM agent_states)"

prod_vol <- dbGetQuery(db, sql)

setDT(prod_vol)

prod_vol[, "model":= as.character(model)]
prod_vol[, "model":= as.character(dbGetQuery(db, sprintf(lookup_sql,
                                    "string_key", "model", paste0("integer_key=", .BY)))), by=model]

print(head(prod_vol))
```

```{r, fig.width=9, fig.height=12}
ggplot(prod_vol, aes(log(oil_prod))) +
    geom_histogram(aes(y=after_stat(density)/uniqueN(prod_vol$RunID), color=as.factor(RunID)), bins=nclass.scott(prod_vol$oil_prod)) +
    stat_function(fun=dnorm, args=market_shares[, .("mean"=mean(log(scaled_oil_BBL)), "sd"=sd(log(scaled_oil_BBL)))]) +
    facet_grid(model~time) +
    labs(title="Distribution of firm oil production over time", color="RunID", y="density")
```

```{r, fig.width=9, fig.height=12}
ggplot(prod_vol, aes(log(gas_prod))) +
    geom_histogram(aes(y=after_stat(density)/uniqueN(prod_vol$RunID), color=as.factor(RunID)), bins=nclass.scott(prod_vol$gas_prod)) +
    stat_function(fun=dnorm, args=market_shares[scaled_gas_MCF>0, .("mean"=mean(log(scaled_gas_MCF)), "sd"=sd(log(scaled_gas_MCF)))]) +
    facet_grid(model~time) +
    labs(title="Distribution of firm gas production over time", color="RunID", y="density")
```

```{r}
agent_states[time==-1,
            {"h_data"= hist(log(gas_flared/oil_output), plot=FALSE);
            plot(exp(h_data$mids), h_data$density, log='x', type='h', lwd=10, lend=2,
            main="Distribution of firm flaring intensity", xlab="Flaring intensity (MCF/BBL)", ylab="Density"); NULL}]
```

## Agent behaviors

```{r}
ggplot(agent_states[behavior!="flaring", .N, keyby=.(model, RunID, time, behavior)]) +
    geom_step(aes(x=time, y=N, group=interaction(model,RunID), color=model)) +
    facet_grid(behavior~., scales="free_y")
```

```{r}
ggplot(agent_states[(behavior!="flaring") & (model=="complete"), .N, keyby=.(RunID, time, behavior)]) +
    geom_step(aes(x=time, y=N, group=RunID)) +
    facet_grid(behavior~., scales="free_y")
```

## Comparison of model components

```{r}
ggplot(agent_states[, sum(gas_flared) / sum(gas_output), by=.(model, RunID, time)]) +
    stat_summary(fun="mean", geom="line", aes(x=time, y=V1, color=model),
                    data=function(x) {x[!((time<0) & (model=="complete"))]}) +
    stat_summary(fun="mean", geom="line", aes(x=time, y=V1), color="black",
                    data=function(x) {x[((time<0) & (model=="complete"))]}) +
    geom_hline(yintercept=0.05, lty=2) +
    geom_hline(yintercept=0.03, lty=2) +
    labs(title="Relative impact of model components", x="Time", y="Total Flaring Intensity, gas basis", color="") +
    theme(legend.position="bottom")
```

```{r}
ggplot(agent_states[, sum(gas_flared) / sum(oil_output), by=.(model, RunID, time)]) +
    stat_summary(fun="mean", geom="line", aes(x=time, y=V1, color=model),
                    data=function(x) {x[!((time<0) & (model=="complete"))]}) +
    stat_summary(fun="mean", geom="line", aes(x=time, y=V1), color="black",
                    data=function(x) {x[((time<0) & (model=="complete"))]}) +
    geom_hline(yintercept=0.1, lty=2) +
    labs(title="Relative impact of model components", x="Time", y="Total Flaring Intensity, oil basis", color="") +
    theme(legend.position="bottom")
```

```{r}
ggplot(agent_states[, sum(gas_flared) / sum(oil_output), by=.(model, RunID, time)]) +
    stat_summary(fun="mean", geom="line", aes(x=time, y=V1, color=model),
                    data=function(x) {x[!((time<0) & (model=="complete"))]}) +
    stat_summary(fun="mean", geom="line", aes(x=time, y=V1), color="black",
                    data=function(x) {x[((time<0) & (model=="complete"))]}) +
    geom_hline(yintercept=0.1, lty=2) +
    labs(title="Relative impact of model components", x="Time", y="Total Flaring Intensity, oil basis", color="") +
    theme(legend.position="bottom") +
    coord_cartesian(ylim=c(0, 0.25))
```

```{r}
ggplot(agent_states[, sum(gas_flared) / sum(oil_output), by=.(model, RunID, time)]) +
    geom_smooth(aes(x=time, y=V1, color=model)) +
    labs(title="Relative impact of model components", x="Time", y="Total Flaring Intensity", color="") +
    theme(legend.position="bottom") +
    coord_cartesian(ylim=c(0, 0.25))
```

```{r}
ggplot(agent_states[model=="complete", sum(gas_flared) / sum(oil_output), by=.(RunID, time)]) +
    geom_smooth(aes(x=time, y=V1), color="black") +
    geom_vline(xintercept=0, lty=2) +
    annotate(x=0, y=Inf, "text", label="Warm start", hjust=1.1, vjust=1.5) +
    annotate(x=0, y=Inf, "text", label="Forecast", hjust=-0.1, vjust=1.5) +
    labs(title="Market evolution", x="Time", y="Total Flaring Intensity", color="") +
    theme_bw() + theme(plot.title=element_text(hjust=0.5)) +
    coord_cartesian(ylim=c(0, 0.25))
```

```{r}
flaring_intensity <- agent_states[(model=="complete"), .("oil"= sum(gas_flared)/sum(oil_output),
                                                        "gas"=sum(gas_flared)/sum(gas_output)), by=.(RunID, time)]

ggplot(melt(flaring_intensity, measure.vars=c("gas","oil"), variable.name="basis")) +
    stat_summary(fun="mean", geom="line", aes(x=time, y=value)) +
    geom_vline(xintercept=0, lty=2) +
    annotate(x=0, y=Inf, "text", label="Warm start", hjust=1.1, vjust=1.5) +
    annotate(x=0, y=Inf, "text", label="Forecast", hjust=-0.1, vjust=1.5) +
    labs(title="Market evolution", x="Time", y="Total Flaring Intensity", color="") +
    facet_wrap(.~basis, scales="free_y", labeller=as_labeller(function(x) paste(x, "basis"))) +
    theme_bw() + theme(plot.title=element_text(hjust=0.5), legend.position="bottom")
```

```{r}
ggplot(agent_states[model=="complete", sum(gas_flared) / sum(oil_output), by=.(RunID, time)]) +
    geom_smooth(aes(x=time, y=V1), color="black") +
    geom_vline(xintercept=0, lty=2) +
    annotate(x=0, y=Inf, "text", label="Warm start", hjust=1.1, vjust=1.5) +
    annotate(x=0, y=Inf, "text", label="Forecast", hjust=-0.1, vjust=1.5) +
    labs(title="Market evolution", x="Time", y="Total Flaring Intensity", color="") +
    theme_bw() + theme(plot.title=element_text(hjust=0.5))
```

```{r}
flaring_intensity <- agent_states[time>=-1, .("oil"= sum(gas_flared)/sum(oil_output),
                                            "gas"=sum(gas_flared)/sum(gas_output)), by=.(model, RunID, time)]

ggplot(melt(flaring_intensity, measure.vars=c("gas","oil"), variable.name="basis")) +
    stat_summary(fun="mean", geom="line", aes(x=time, y=value, color=model),
                    data=function(x) {x[!((time<0) & (model=="complete"))]}) +
    stat_summary(fun="mean", geom="line", aes(x=time, y=value), color="black",
                    data=function(x) {x[((time<0) & (model=="complete"))]}) +
    labs(title="Relative impact of model components", x="Time", y="Total Flaring Intensity", color="") +
    facet_wrap(.~basis, scales="free_y", labeller=as_labeller(function(x) paste(x, "basis"))) +
    theme_bw() + theme(plot.title=element_text(hjust=0.5), legend.position="bottom")
```


## Rate of asset retirement and acquisition
```{r}
sql <- "SELECT %1$s, model, RunID, %2$s FROM lease_info GROUP BY model, RunID, %2$s"

new_prod <- dbGetQuery(db, sprintf(sql, "SUM(oil_BBL+cond_BBL) AS new_BBL", "t_found"))
ret_prod <- dbGetQuery(db, sprintf(sql, sprintf("SUM(CASE WHEN status=(%s) THEN oil_BBL+cond_BBL ELSE 0 END) AS ret_BBL",
                                   sprintf(lookup_sql, "integer_key", "status", "string_key='retired'")), "t_last"))

setDT(new_prod)
setDT(ret_prod)
prod <- ret_prod[new_prod, on=c("model", "RunID", t_last="t_found")][t_last>=min(ret_prod$t_last)]

prod[, "model":= as.character(model)]
prod[, "model":= as.character(dbGetQuery(db, sprintf(lookup_sql,
                                    "string_key", "model", paste0("integer_key=", .BY)))), by=model]

setnames(prod, "t_last", "time")
setnafill(prod, cols=c("ret_BBL", "new_BBL"), fill=0)

prod[, "ratio":= new_BBL / ret_BBL]
prod[(new_BBL==0) & (ret_BBL==0), "ratio":= 0]

print(prod)
```

```{r, fig.height=7.5}
ggplot(prod[time>min(time)], aes(x=time, y=ratio, group=model)) +
    geom_line() +
    geom_hline(yintercept=1.23, lty=3, color="red") +
    stat_summary(aes(x=1, yintercept=..y..), fun="median", na.rm=TRUE, geom="hline", lty=2) +
    stat_summary(aes(x=10, label=sprintf("Median: %.02f",..y..)), fun="median", na.rm=TRUE, geom="text", vjust=-2.5) +
    labs(title="Oil production gained and lost from new and retiring wells", y="Production gained / Production lost") +
    facet_grid(model~.) +
    coord_cartesian(ylim=c(0, 10))
```

```{r, fig.height=7.5}
ggplot(prod[(time>min(time)) & (model=="complete")], aes(x=time, y=ratio, group=model)) +
    geom_line() +
    geom_hline(yintercept=1.23, lty=3, color="red") +
    stat_summary(aes(x=1, yintercept=..y..), fun="median", na.rm=TRUE, geom="hline", lty=2) +
    stat_summary(aes(x=10, label=sprintf("Median: %.02f",..y..)), fun="median", na.rm=TRUE, geom="text", vjust=-2.5) +
    labs(title="Oil production gained and lost from new and retiring wells", y="Production gained / Production lost") +
    facet_grid(model~.) +
    coord_cartesian(ylim=c(0, 10))
```

## Market shares
#### CRSP/Compustat Merged (Annual) for oil and gas firms
```{r}
CCM_OG <- fread("./inputs/processed/firm_finances.csv")
```

```{r}
ggplot(CCM_OG) +
    geom_histogram(aes(x=market_share)) +
    scale_x_log10() +
    facet_wrap(fyear~.)
```

### Joint distribution of market share and market value
```{r}
ggplot(CCM_OG) +
    geom_density2d_filled(aes(x=market_share, y=CI)) +
    scale_x_log10() +
    scale_y_log10()
```

```{r}
ggplot(agent_states[time==max(time)]) +
    geom_density(aes(x=gas_flared/oil_output, color=model)) +
    scale_x_log10() +
    theme_bw() + theme(plot.title=element_text(hjust=0.5), legend.position="bottom")
```

```{r}
ggplot(agent_states[(time==max(time)) & (model=="complete")], aes(x=1+market_value, y=gas_flared/oil_output)) +
    geom_density_2d_filled(show.legend=FALSE) +
    geom_hline(yintercept=0.05, color="white", linetype=2) +
    scale_x_log10() + scale_y_log10() +
    annotate("text", x=-Inf, y=0.05, vjust=-1, label="Green market\nthreshold", color="white") +
    labs(title="Density of all firms", x="Market Value", y="Flaring Intensity") +
    theme_bw() + theme(plot.title=element_text(hjust=0.5)) + coord_cartesian(expand=FALSE)
```

```{r}
ggplot(agent_states[time==max(time)], aes(x=1+market_value, y=gas_flared/oil_output)) +
    geom_density_2d_filled(show.legend=FALSE) +
    geom_hline(yintercept=0.05, color="white", linetype=2) +
    scale_x_log10() + scale_y_log10() +
    facet_wrap(model~.) +
    annotate("text", x=0, y=0.05, vjust=-1, label="Green market\nthreshold", color="white") +
    labs(title="Density of all firms", x="Market Value", y="Flaring Intensity") +
    theme_bw() + theme(plot.title=element_text(hjust=0.5)) + coord_cartesian(expand=FALSE)
```

```{r}
ggplot(agent_states[(time==max(time)) & (behavior!="flaring") & (model=="complete")], aes(x=1+market_value)) +
    geom_density(aes(color=behavior)) +
    labs(title="Density of non flaring firms", x="Market Value", color="Behavior:") +
    scale_x_log10() +
    theme_bw() + coord_cartesian(expand=FALSE) +
    theme(plot.title=element_text(hjust=0.5), legend.position=c(0.8,0.8), legend.background=element_rect(fill="white",color="black"))
    #theme_bw() + theme(plot.title=element_text(hjust=0.5), legend.position=c(0.8,0.8)) + coord_cartesian(expand=FALSE)
```

## Market Dynamics
```{r}
sql <- "SELECT %1$s.gas_MCF, %1$s.csgd_MCF, %1$s.opEx_pMCF, %2$s.market, %2$s.RunID, %2$s.model FROM %2$s 
        INNER JOIN %1$s ON %2$s.leaseID=%1$s.leaseID AND %1$s.RunID=%2$s.RunID AND %1$s.model=%2$s.model 
            AND %2$s.time=(SELECT MAX(time) FROM %2$s) 
            AND %2$s.market IN (%3$s) 
            AND (%1$s.gas_MCF>0 OR %1$s.csgd_MCF>0) 
        ORDER BY opEx_pMCF ASC"

gas_supply <- dbGetQuery(db, sprintf(sql, "lease_info", "lease_states",
                               sprintf(lookup_sql, "integer_key", "market", "string_key IN ('green','grey')")))

setDT(gas_supply)
gas_supply[, "model":= as.character(model)]
gas_supply[, "model":= as.character(dbGetQuery(db, sprintf(lookup_sql,
                                    "string_key", "model", paste0("integer_key=", .BY)))), by=model]

gas_supply[, "market":= as.character(market)]
gas_supply[, "market":= as.character(dbGetQuery(db, sprintf(lookup_sql,
                                    "string_key", "market", paste0("integer_key=", .BY)))), by=market]

gas_supply[, "q_sup":= cumsum(gas_MCF+csgd_MCF), by=.(model, RunID, market)]
gas_supply[market_history, on=c("model","RunID"), "max_prop_green":= market_prop_green]

gas_demand <- unserialize(dbGetQuery(db,
        "SELECT * FROM demand_functions ORDER BY RANDOM()")$fun_bin[[1]])

gas_supply[, c("p_grey","p_green"):= gas_demand$new_schedule(first(max_prop_green))(q_sup)[, .(p_grey,p_green)], by=model]
```

```{r, fig.height=12}
ggplot(gas_supply, aes(x=q_sup, color=market)) +
    geom_line(aes(y=opEx_pMCF, group=RunID)) +
    geom_line(aes(y=ifelse(market=="green", p_green, p_grey), color="Sample demand"), linetype="dashed") +
    scale_color_manual(values=c("grey"="grey40", "green"="darkgreen", "Sample demand"="black"), guide="none") +
    facet_grid(model~market, scales="free_x") +
    labs(x="Q (MCF / month)", y="Price ($)", title="Natural gas market", color="")
```

```{r}
ggplot(gas_supply[model=="complete"][RunID==RunID[which.max(opEx_pMCF)]], aes(x=q_sup, color=market)) +
    geom_line(aes(y=opEx_pMCF), size=0.85) +
    geom_line(aes(y=p_grey, color="grey"), linetype="dashed", size=0.85) +
    geom_line(aes(y=replace(p_green, p_green==0, NA), color="green"), linetype="dashed", size=0.85) +
    scale_color_manual(values=c("grey"="grey40", "green"="darkgreen")) +
    scale_x_continuous(labels = function(x) format(x, scientific = TRUE)) +
    labs(x="Quantity (MCF / month)", y="Price ($)", title="Sample Natural Gas Market", color="Market")
```

### Emergent gas market prices
```{r, fig.height=7.5}
ggplot(market_history, aes(x=time, group=RunID)) +
    geom_line(aes(y=p_grey, color="grey"), alpha=0.4) +
    geom_line(aes(y=p_green, color="green"), alpha=0.4) +
    scale_color_manual(values=c("grey"="grey40", "green"="darkgreen")) +
    facet_grid(model~.) +
    labs(title="Gas market price", y="$ / MCF", color="Market")
```

Expect:
$$1.07 P_{grey} <= P_{green} <= 1.3 P_{grey}$$
```{r, fig.height=7.5}
ggplot(market_history[p_green>0], aes(x=time, group=RunID, color=as.factor(RunID))) +
    geom_point(aes(y=p_green/p_grey), alpha=0.4) +
    geom_rug(aes(y=p_green/p_grey), alpha=0.1, color="grey40", sides='r', outside=TRUE) +
    coord_cartesian(clip='off') +
    geom_hline(yintercept=c(1.07, 1.3), color="grey40", linetype="dashed") +
    facet_grid(model~.) +
    labs(title="Green gas premium", y=expression(P[green] / P[grey]), color="RunID") +
    theme_bw() + theme(plot.title=element_text(hjust=0.5), strip.switch.pad.grid=unit(0.3,'cm'), strip.placement='outside')
```

### Comparison with emprical values
```{r}
monthly_prices <- fread("./inputs/processed/historical_NG_demand.csv")
monthly_prices[, "Date":= as.Date(paste("01", month, year), format='%d %b %Y')]

ggplot(monthly_prices[year>2010], aes(x=Date, y=p)) +
    geom_line() +
    geom_hline(aes(yintercept=mean(p)), color="blue") +
    labs(title="Empirical monthly natural gas prices since 2010", y="$ per MCF")
```

```{r, fig.height=7.5}
ggplot(market_history, aes(x=time, y=p_grey, group=RunID, color=as.factor(RunID))) +
    geom_line() +
    geom_hline(aes(yintercept=mean(p_grey)), color="blue") +
    facet_grid(model~.) +
    labs(title="Gas market price", y="$ / MCF", color="RunID")
```

```{r, fig.width=10}
par(mgp=c(3, 1.5, 0))
boxplot(data.table(matrix(nrow=0, ncol=market_history[time<0, uniqueN(model) + 1])),
        boxfill=NA, border=NA, ylim=c(1,6), ylab="$ / MCF", xaxt="n")
title("Emprical & modeled distribution of conventional gas prices")
abline(h=1:6, col = "grey", lty = "dotted")
boxplot(monthly_prices[year>2010]$p, boxwex=1.6, add=TRUE, at=1)
axis(1, at=1, labels="empirical", col.axis="blue")
boxplot(p_grey~gsub("\\n", " ", model), data=market_history[time<0],
    add=TRUE, at=c(seq(market_history[time<0, uniqueN(model)]) + 1))
```

### Gas market size
```{r, fig.height=7.5}
ggplot(market_history, aes(x=time, group=RunID, color=as.factor(RunID))) +
    geom_point(aes(y=q_green/(q_grey+q_green)), alpha=0.4) +
    geom_line(aes(y=market_prop_green), linetype="dashed", color="black") +
    facet_grid(model~.) +
    labs(title="Green gas market size", y="% of total sales", color="RunID")
```

```{r, fig.height=7.5}
ggplot(market_history, aes(x=time, group=RunID, color=as.factor(RunID))) +
    geom_point(aes(y=q_green+q_grey), alpha=0.4) +
    facet_grid(model~.) +
    labs(title="Total gas sold", y="MCF", color="RunID")
```

```{r}
ggplot(market_history[model %in% c("complete","no differentiation")], aes(x=time, group=RunID, color=as.factor(RunID))) +
    geom_point(aes(y=q_green+q_grey), alpha=0.4) +
    facet_grid(model~.) +
    labs(title="Total gas sold", y="MCF", color="RunID")
```

Storage and withdrawls compared to scaled [Texas capacity](https://www.eia.gov/dnav/ng/ng_sum_lsum_dcu_STX_m.htm)
```{r, fig.height=7.5}
ggplot(market_history[!is.na(q_stored)][order(time),
                        .(time, "net"= cumsum(q_stored), "frac"= mean(frac)), by=.(model, RunID)]) +
    geom_line(aes(x=time, y=net, group=RunID, color=as.factor(RunID)), alpha=0.4) +
    geom_hline(aes(yintercept= mean(frac) * 400000 * 1000), color='darkred', linetype='dashed') +
    annotate("text", x=0, y=0.9*400000*1000*mean(market_history$frac),
            label="Scaled Working Gas Storage", color='darkred', hjust=1) +
    facet_grid(model~.) +
    labs(title="Net gas stored and withdrawn", y="MCF", color="RunID") +
    theme_bw() + theme(plot.title=element_text(hjust=0.5))
```

# Shareholder valuation
```{r}
agent_states[params, on=c("model","RunID","time"), "SRoR":= SRoR]

agent_states[(cost_M>0) & (time>0) & (model!="no shareholder\nvaluation"),
                        {MV0i= pmax(market_value - (cost_M*SRoR), 1e-10);
                        SRoR * (cost_M/MV0i) / cut(cost_M/MV0i, 10, labels=FALSE)}, by=.(RunID, time)][, summary(V1)]
```

Shareholder valuation ($\theta$) of CSR spending ($C$) should increase market value ($MV$) proportional to CSR rating ($R \in [1-10]$)  
$$ MV_0 (1 + \alpha R) \simeq MV_0 + \theta C$$  
[Luo and Bhattacharya, 2006](10.1509/jmkg.70.4.001)  
```{r}
ggplot(agent_states[model!="no shareholder\nvaluation",
                    cbind(.SD, "Ci" = cost_M*SRoR)][
                    (cost_M>0) & (time>0) & (market_value>Ci),
                            {"MV0i" = market_value - Ci;
                            Ci / MV0i / cut(cost_M/MV0i, 10, labels=FALSE)}, by=.(RunID, time)]) +
    geom_boxplot(aes(x=as.factor(RunID), y=V1, color=as.factor(RunID)), outlier.shape=NA, show.legend=FALSE) +
    geom_hline(yintercept=0.18*0.19, lty=2) +
    geom_hline(yintercept=0.11+(0.18*0.19), lty=2) +
    labs(title="Effect of shareholder valuation on market value", x="RunID", y="% Increase in Market Value") +
    theme_bw() + theme(plot.title=element_text(hjust=0.5)) + coord_cartesian(ylim=c(0,0.5))
```