flaringABM_validation_init.Rmd

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

```{r}
library(ggplot2)
library(data.table)

# better break points for adjusted log scale
log1p_breaks <- function(x, n) { replace(axisTicks(c(0, log10(max(x))), log=TRUE, n=n), 1, 0) }
```

# Lease expenses
## Capital costs  
Gas leases should have only one well
```{r}
print(leases_emp[OIL_GAS_CODE=="G", summary(N)])
```

Expect a capital expenditures to be about [$7/BOE](http://graphics.wsj.com/oil-barrel-breakdown/)
```{r}
leases[, cat("All leases:", sum(capEx) / sum(EUR), "\t Assigned leases:", .SD[!is.na(firmID), sum(capEx) / sum(EUR)])]
```

## Operating costs per MCF of gas at oil and gas leases  
```{r, fig.cap="Gas lease operating expenditures are skewed left (cheaper) while oil leases are skewed right"}
ggplot(leases[gas_MCF+csgd_MCF>0], aes(color=ifelse(OIL_GAS_CODE=="O", "Oil", "Gas"))) +
    geom_density(aes(x=opEx_pMCF + ifelse(csgd_MCF==0, 0, capEx_csgd / csgd_MCF / lifetime))) +
    theme_bw() + theme(plot.title=element_text(hjust=0.5), legend.position=c(0.85, 0.8), legend.background=element_rect(color="black")) +
    labs(color="Lease type", x="Average Cost per MCF", title="Gas operating costs") +
    xlim(0, 5)
```

## Operating costs per BBL of oil at oil leases  
```{r}
fed_data <- fread(text="area,           year, Mean,  Min, Max
                        Eagle Ford,     2017, 29,    10,  50
                        Eagle Ford,     2021, 16.58, 7,   35
                        Delaware Basin, 2017, 33,    15,  55
                        Delaware Basin, 2021, 26.19, 7,   50
                        Midland Basin,  2017, 24,    10,  55
                        Midland Basin,  2021, 27.05, 6,   60")

ggplot(leases[oil_BBL>0], aes(color=area)) +
    geom_density(data= function(x) {x[,"area":= ifelse(area=="Spraberry", "Midland Basin", area)]},
                aes(x=cost_oil/oil_BBL, color=area)) +
    geom_point(aes(x=Mean, y=(0.5 + (year==2021))), data=fed_data) +
    geom_errorbarh(aes(xmin=Min, xmax=Max, y=(0.5 + (year==2021))), alpha=0.5, lty=2, data=fed_data, show.legend=FALSE) +
    annotate("text", x=fed_data[year==2017, max(Mean)], y=0.5, label="2017", alpha=0.5, vjust=-1) +
    annotate("text", x=fed_data[year==2021, min(Mean)], y=1.5, label="2021", alpha=0.5, vjust=-1) +
    scale_x_log10() +
    theme_bw() + theme(plot.title=element_text(hjust=0.5), legend.position=c(0.85, 0.8), legend.background=element_rect(color="black")) +
    labs(x="$ / barrel", title="Oil operating costs at oil leases", color="")
```

```{r, echo=FALSE}
knitr::kable(fed_data, col.names=c("Region", "Year", "Mean", "Min", "Max"),
                caption=paste("Break-even operating expenses ($/barrel) reported by Dallas Fed",
                                "([2017](https://www.dallasfed.org/research/surveys/des/2017/1701.aspx#tab-questions),",
                                "[2021](https://www.dallasfed.org/research/surveys/des/2021/2101.aspx#tab-questions)):"))
```

# Lease assignment summary  
### Number of assigned leases   
```{r}
print(leases[!is.na(firmID), .N, by=OIL_GAS_CODE])
```

### Number of leases left to be discovered  
```{r}
print(leases[is.na(firmID), .N, by=OIL_GAS_CODE])
```

### Lifetime estimation
Goodness of fit among historical leases
```{r}
leases_emp[(expiration<201912) & (start>199301),
            .("R2"= summary(lm(log(m_rem)~log(q_i/q_curve)))$r.squared), by=area]
```

```{r}
leases_emp[(expiration<201912) & (start>199301),
            with(summary(lm(log(m_rem)~log(q_i/q_curve))),
                setNames(.(r.squared, coefficients[2,1]), c("R2", "slope"))), by=area]
```

```{r}
ggplot(leases_emp, aes(x=q_i/EUR, y=lifetime)) +
    geom_point(aes(color=ifelse((expiration<201912) & (start>199301),"Observed","Predicted")), alpha=0.3) +
    scale_x_log10() + scale_y_log10() +
    theme_bw() + theme(plot.title=element_text(hjust=0.5)) +
    facet_grid(.~area) +
    labs(title="Lease lifetimes based on production profile", color="")
```

### Expiring & new production
```{r}
ggplot(leases[!is.na(firmID), sum(oil_BBL + cond_BBL), by=.("ex"=lifetime + t_found - Params$t0)]) +
    geom_col(aes(x=ex, y=V1)) +
    geom_hline(aes(yintercept=median(V1, na.rm=TRUE)), lty=2) +
    geom_text(data=function(x) {x[,.(median(V1, na.rm=TRUE))]},
         aes(y=V1, label=sprintf("Median: %.02f", V1)), x=-Inf, hjust=-0.01, vjust=-1) +
    labs(title="Quantity of expiring oil production", y="Oil + Condensate (BBL)", x="Time (Months)") +
    coord_cartesian(xlim=c(0, Params$tf - Params$t0))
```

```{r}
cat("Potential new production\n")
print(leases[is.na(firmID), summary(oil_BBL + cond_BBL)])
```

### Percent of operators with no gas leases  
```{r}
# Empirical
print(market_shares[, sum(scaled_gas_MCF==0)/.N])
# Generated
print(leases[, sum(gas_MCF), by=firmID][, sum(V1==0)/.N])
```

# Distributions of operator production volumes
```{r}
par(mfrow=c(1,2))
hist(log(market_shares$scaled_oil_BBL), breaks="scott", freq=FALSE,
    main="Empirical", xlab="log(Oil (BBL))", xlim=c(0,20), ylim=c(0,0.2))
market_shares[, {curve(dnorm(x, mean(log(scaled_oil_BBL)), sd(log(scaled_oil_BBL))), col="blue", add=TRUE); NULL}]
hist(log(leases[!is.na(firmID), sum(oil_BBL), by=firmID]$V1), breaks="scott", freq=FALSE,
    main="Generated", xlab="log(Oil (BBL))", xlim=c(0,20), ylim=c(0,0.2))
market_shares[, {curve(dnorm(x, mean(log(scaled_oil_BBL)), sd(log(scaled_oil_BBL))), col="blue", add=TRUE); NULL}]
```

```{r}
par(mfrow=c(1,2))
hist(log(market_shares$scaled_gas_MCF), breaks="scott", freq=FALSE,
    main="Empirical", xlab="log(Gas (MCF))", xlim=c(0,20), ylim=c(0,0.2))
market_shares[scaled_gas_MCF>0, .(log(scaled_gas_MCF))][, {curve(dnorm(x, mean(V1), sd(V1)), col="blue", add=TRUE); NULL}]
hist(log(leases[!is.na(firmID), sum(gas_MCF), by=firmID]$V1), breaks="scott", freq=FALSE,
    main="Generated", xlab="log(Gas (MCF))", xlim=c(0,20), ylim=c(0,0.2))
market_shares[scaled_gas_MCF>0, .(log(scaled_gas_MCF))][, {curve(dnorm(x, mean(V1), sd(V1)), col="blue", add=TRUE); NULL}]
```

Joint density of oil and (casinghead) gas outputs
```{r, fig.width=10}
ggplot(melt(market_shares, id.vars="scaled_oil_BBL", measure.vars=patterns("_MCF$"), variable.name="type")) +
    stat_density_2d(geom="raster", aes(x=scaled_oil_BBL, y=value, fill=after_stat(density)), contour=FALSE) +
    stat_bin_2d(geom="point", aes(x=V1, y=V2, size=after_stat(count), color="red", fill=NULL), alpha=0.5,
        data=leases[!is.na(firmID), .(sum(oil_BBL), c(sum(gas_MCF), sum(csgd_MCF)),
                                    "type"=c("scaled_gas_MCF","scaled_csgd_MCF")), by=firmID]) +
    scale_size_area(guide=FALSE) +
    scale_color_manual("", labels="Generated", values="red", guide=guide_legend(override.aes=list("alpha"=1))) +
    scale_fill_viridis_c("Empirical") +
    scale_x_continuous(trans=scales::log1p_trans(), breaks= function(x) {log1p_breaks(x, 5)}) +
    scale_y_continuous(trans=scales::log1p_trans(), breaks= function(x) {log1p_breaks(x, 5)}) +
    facet_grid(.~type, labeller=as_labeller(function(x) ifelse(x=="scaled_gas_MCF", "Unassociated Gas", "Casinghead Gas"))) +
    theme_bw() + theme(plot.title=element_text(hjust=0.5)) + coord_cartesian(expand=FALSE) +
    labs(title="Joint distribution of firm oil and gas production", x="Production (BBL)", y="Production (MCF)")
```

# Market conditions
## Oil
### Supply curve
```{r}
ggplot(leases[!is.na(firmID) & oil_BBL>0][order(opEx_pBBL)]) +
    geom_line(aes(x=cumsum(oil_BBL), y=opEx_pBBL)) +
    scale_y_log10() +
    labs(x="Q", y="P", title="Oil supply curve")
```

## Natural gas supply and demand
### Supply curve
```{r}
ggplot(leases[class=="developed"][order(opEx_pMCF)]) +
    geom_line(aes(x=cumsum(gas_MCF), y=opEx_pMCF, color="Only un-associated gas")) +
    geom_line(aes(x=cumsum(csgd_MCF + gas_MCF), y=opEx_pMCF, color="Including casinghead gas")) +
    theme(legend.position=c(0.2, 0.8), legend.background=element_blank()) +
    labs(x="Q (MCF / month)", y="Marginal Cost ($)", title="Natural gas supply curve", color="")
```

### Demand curve
```{r, fig.cap="Sample demand curves using single months. Model demand curves use samples of one or more months"}
ggplot(demand$historical_market, aes(color=factor(month, levels=month.abb))) +
    geom_point(aes(x=q, y=p), shape=4) +
    geom_line(data= function(x) {x[sample(.N, 5), {"shift"= q - (q_TX*frac*0.95);
                                                    cbind(shift, demand$new_schedule(0, sample_set=.BY)(x$q - shift))},
                                    by=.(year, month)][p_grey>0]},
            aes(x=q+shift, y=p_grey, group=interaction(year, month)), show.legend=FALSE) +
    geom_text(aes(q, p, label=substr(year,3,4)), size=3, hjust=0) +
    labs(title="Sampling of representative demand curves", x="Quantity (MCF)", y="Price ($/MCF)", color="month")
```

### Emergent market
```{r}
ggplot(leases[!is.na(firmID)][order(opEx_pMCF)]) +
    geom_line(aes(x=cumsum(gas_MCF), y=opEx_pMCF, color="Only un-associated gas")) +
    geom_line(aes(x=cumsum(csgd_MCF + gas_MCF), y=opEx_pMCF, color="Including casinghead gas"),
                data=leases[class=="developed"][order(opEx_pMCF)]) +
    geom_line(aes(x=cumsum(csgd_MCF + gas_MCF), y=opEx_pMCF, color="Assuming 100%\ncasinghead gas capture")) +
    geom_line(data= function(x) {rbindlist(lapply(1:5, function(z) demand$new_schedule(0)(x[, cumsum(gas_MCF+csgd_MCF)])), idcol=TRUE)},
                aes(x=q, y=p_grey, group=.id, color="Sample demand")) +
    theme(legend.position=c(0.85, 0.15), legend.background=element_blank()) +
    ylim(0,7) +
    labs(x="Q (MCF / month)", y="Price ($)", title="Natural gas market", color="")
```

```{r}
ggplot(leases[firms[behavior!="flaring"], on="firmID"][order(opEx_pMCF)]) +
    geom_line(aes(x=cumsum(gas_MCF), y=opEx_pMCF, color="Only un-associated gas")) +
    geom_line(aes(x=cumsum(csgd_MCF + gas_MCF), y=opEx_pMCF, color="Assuming 100%\ncasinghead gas capture")) +
    geom_line(data= function(x) {rbindlist(lapply(1:5, function(z) demand$new_schedule(0.1)(x[, cumsum(gas_MCF+csgd_MCF)])), idcol=TRUE)},
                aes(x=q, y=p_green, group=.id, color="Sample green demand")) +
    theme(legend.position=c(0.85, 0.15), legend.background=element_blank()) +
    ylim(0,7) +
    labs(x="Q (MCF / month)", y="Price ($)", title="Green natural gas market", color="")
```

### Flaring intensity floor
```{r}
leases[, sum(sopf_MCF) / sum(oil_BBL + cond_BBL)]
leases[!is.na(firmID), sum(sopf_MCF) / sum(oil_BBL + cond_BBL)]
leases[is.na(firmID), sum(sopf_MCF) / sum(oil_BBL + cond_BBL)]
```