# Type Curve Economics

In **Part 1** we downloaded production data and manipulated it. In **Part 2**, we learned how to create type-curves. Now, let’s create a quick single-well economic model, and I suppose play around with creating **functions** for the first time.

**Follow My Blog**

**Functions**

Not as bad as you think. Generally the rule of thumb is that if you write the same bit of code more than **3** times, then you should be using a **function**. Of course, if you’ve been following along with the tutorials so far then you’ve already been using functions. The pre-built R packages are really just a lot of bundled functions.

But now, we’ll just build one ourself.

**Revisiting aRpsDCA**

You will recall our single-segment **Arps Modified Hyperbolic Decline** curve from ** Part 2**.

```
library(aRpsDCA)
qiGas <- 20000 #mcfd
bGas <- 1 #B-factor
DiGas <- 0.9 #Effective Annual Decline (Initial)
DfGas <- 0.08 #Effective Annual Decline (Terminal)
curtail <- 5 #Months of curtailment
wellLife <- 30 #Years
gasFcst = curtailed.q(arps.decline(
qiGas*365.25, as.nominal(DiGas), bGas,
as.nominal(DfGas)),
curtail/12.0,seq(1/12, wellLife, by= (1/12)))/12
```

There are actually three functions here, **curtailed.q**, **arps.decline**, and **as.nominal**. All are defined within the package documentation for ** aRpsDCA**. These allow me to forecast well production. What I want to do is create a single function that will spit out a dataframe of oil, gas, and water production in monthly increments.

**Our Function**

There are two things I need to do:

- Determine the variables within
**function()** - Create the calculation after
**function()**

So I will just use the above code as my guide.

```
prodFunc <- function(wellLife,
qiGas, DiGas, bGas, DfGas, curtailGas,
qiOil, DiOil, bOil, DfOil, curtailOil,
qiWater, DiWater, bWater, DfWater, curtailWater) {
df1 <- data.frame(
#Months of Prouction
monthsOn = seq(1, wellLife*12, 1),
#Define Gas
gas = curtailed.q(arps.decline(
qiGas*365.25, as.nominal(DiGas), bGas,
as.nominal(DfGas)),
curtailGas/12.0,seq(1/12, wellLife, by= (1/12)))/12,
#Define Oil
oil = curtailed.q(arps.decline(
qiOil*365.25, as.nominal(DiOil), bOil,
as.nominal(DfOil)),
curtailOil/12.0,seq(1/12, wellLife, by= (1/12)))/12,
#Define Water
water = curtailed.q(arps.decline(
qiWater*365.25, as.nominal(DiWater), bWater,
as.nominal(DfWater)),
curtailWater/12.0,seq(1/12, wellLife, by= (1/12)))/12
)
return(df1)
}
```

**Calculation**

Now that our function is build, let’s test it out. To avoid confusion, I’m going to throw a **1** behind each **variable** I define.

```
qiGas1 <- 2000 #mcfd
bGas1 <- 1.2 #B-factor
DiGas1 <- 0.8 #Effective Annual Decline (Initial)
DfGas1 <- 0.08 #Effective Annual Decline (Terminal)
curtailGas1 <- 3 #Months of curtailment
qiOil1 <- 1000 #bopd
bOil1 <- 1 #B-factor
DiOil1 <- 0.9 #Effective Annual Decline (Initial)
DfOil1 <- 0.08 #Effective Annual Decline (Terminal)
curtailOil1 <- 1 #Months of curtailment
qiWater1 <- 1000 #bwpd
bWater1 <- 1 #B-factor
DiWater1 <- 0.9 #Effective Annual Decline (Initial)
DfWater1 <- 0.08 #Effective Annual Decline (Terminal)
curtailWater1 <- 0 #Months of curtailment
wellLife1 <- 30 #Years
#Create a dataframe with our function
prod <- prodFunc(wellLife = wellLife1,
qiGas = qiGas1, DiGas = DiGas1, bGas = bGas1,
DfGas = DfGas1, curtailGas = curtailGas1,
qiOil = qiOil1, DiOil = DiOil1, bOil = bOil1,
DfOil = DfOil1, curtailOil = curtailOil1,
qiWater = qiWater1, DiWater = DiWater1, bWater = bWater1,
DfWater = DfWater1, curtailWater = curtailWater1)
head(prod)
```

I get this.

monthsOn | gas | oil | water |
---|---|---|---|

1 | 60875 | 30438 | 25537 |

2 | 60875 | 25537 | 21996 |

3 | 60875 | 21996 | 19317 |

4 | 53756 | 19317 | 17220 |

5 | 48244 | 17220 | 15534 |

6 | 43839 | 15534 | 14148 |

If we sum each column, that should give us our **EUR (Estimated Ultimate Recovery)**.

```
sum(prod$oil)/1000 #Convert Oil EUR from BBL to MBO
```

**666.5 MBO** for this example.

**Putting it All Together**

**Prices**

Now that I can create a production function, I can start tacking on additional functionality. For this example, I am going to use a forward strip price from the **Wall Street Journal **for my forecast, as I covered in **Scraping Strip with R (Oil & Gas Coding Series)**. But this time I’ll make it a **function**!

```
library(rvest)
library(httr)
library(dplyr)
library(zoo)
priceFunc <- function() {
crude = 'https://quotes.wsj.com/futures/CRUDE%20OIL%20-%20ELECTRONIC/contracts'
webpage <- read_html(crude)
tbls <- rvest::html_nodes(webpage, 'table')
tbls_ls <- webpage %>%
rvest::html_nodes('table') %>%
.[1] %>%
rvest::html_table(fill = TRUE)
wti <- tbls_ls[[1]] %>% select(MONTH, SETTLEMENT)
wti <- wti %>% filter(MONTH != 'Front Month')
wti$MONTH <- gsub('Crude Oil', '', wti$MONTH)
wti$MONTH <- trimws(wti$MONTH, which = 'both')
wti$Year <- substr(wti$MONTH, 5, 8)
wti$Month <- substr(wti$MONTH, 1, 3)
wti$Date <- as.Date(paste0(wti$Month, '/01/', wti$Year), format = '%b/%d/%Y')
wti <- wti %>% select(DATE = Date, WTI = SETTLEMENT)
gas = 'https://quotes.wsj.com/futures/NATURAL%20GAS/contracts'
webpage <- read_html(gas)
tbls <- rvest::html_nodes(webpage, 'table')
tbls_ls <- webpage %>%
rvest::html_nodes('table') %>%
.[1] %>%
rvest::html_table(fill = TRUE)
hh <- tbls_ls[[1]]%>% select(MONTH, SETTLEMENT) %>% filter(MONTH != 'Front Month')
hh$MONTH <- gsub('Natural Gas', '', hh$MONTH)
hh$MONTH <- trimws(hh$MONTH, which = 'both')
hh$Year <- substr(hh$MONTH, 5, 8)
hh$Month <- substr(hh$MONTH, 1, 3)
hh$Date <- as.Date(paste0(hh$Month, '/01/', hh$Year), format = '%b/%d/%Y')
hh <- hh %>% select(DATE = Date, HH = SETTLEMENT)
wti <- wti %>% left_join(hh)
wti <- wti %>% filter(!is.na(WTI)) %>% filter(!is.na(HH))
return(wti)
}
```

Then, to call my price file, I just run:

```
prices <- priceFunc()
```

So, now if I want to rerun my price file, I can just run **priceFunc()** and then join it to my production dataframe.

You could potentially throw this function into our original **prodFunc()**, but that is not preferable. You will eventually be running a lot of different variations of type-curves in the future, so it is preferable to only run the **priceFunc()** at the beginning of your session. Otherwise you will be trying to scrape prices each time you change a type-curve, and that is **SLOW**.

Below would be the preferred workflow.

```
prod <- prodFunc()
prices <- priceFunc()
#Add similar monthsOn column to prices
prices <- prices %>% mutate(monthsOn = seq(1, n(), 1))
#Join prices to production file on monthsOn column
prod <- prod %>% left_join(prices)
#Drag missing values downwards
prod$WTI <- zoo::na.locf(prod$WTI)
prod$HH <- zoo::na.locf(prod$HH)
#View first rows
head(prod)
```

monthsOn | gas | oil | water | DATE | WTI | HH |
---|---|---|---|---|---|---|

1 | 60875 | 30438 | 25537 | 2020-10-01 | 38.05 | 2.406 |

2 | 60875 | 25537 | 21996 | 2020-11-01 | 38.41 | 2.873 |

3 | 60875 | 21996 | 19317 | 2020-12-01 | 38.87 | 3.241 |

4 | 53756 | 19317 | 17220 | 2021-01-01 | 39.34 | 3.361 |

5 | 48244 | 17220 | 15534 | 2021-02-01 | 39.78 | 3.323 |

6 | 43839 | 15534 | 14148 | 2021-03-01 | 40.18 | 3.201 |

**Revenue**

Now I will just add in a quick revenue model based on absolutely nothing but probably close enough. What will I need for this?

**Working Interest/Net Revenue Interest****Shrink (or as I define it Gas Sold as % of Gas Produced)****NGL Yield (BBL/MMCF of wellhead gas)****Diff’s (I’ll just use simple percents here but can be per-unit)**

```
wi <- 1
nri <- 0.75 #To the 100% of WI
shrink <- 0.75
nglYield <- 100
oilDiff <- 0.95
gasDiff <- 0.8
nglDiff <- 0.3 #%WTI
prod$ngl <- prod$gas*nglYield/1000
prod$oilRevenue <- prod$oil*wi*nri*prod$WTI*oilDiff
prod$gasRevenue <- prod$gas*wi*nri*prod$HH*gasDiff*shrink
prod$nglRevenue <- prod$ngl*wi*nri*prod$WTI*nglDiff
prod$revenue <- prod$oilRevenue+prod$gasRevenue+prod$nglRevenue
```

**Expenses**

Expenses are usually based on the **Working Interest** ownership, so that’s what I will use here. There is a lot you can do, but for my analysis I am completely making it up (or maybe it’s influenced by actuals). For this model I will assume:

- Fixed Cost at $10,000 $/month until total fluid rate goes below 100 barrels per day, or 3,000 per month. $5,000 thereafter.
- Variable Expenses of $0.5/bbl oil, $0.2/mcf wellhead gas, $1/bbl water.

You could potentially run assumed Workover Expenses as well, but I usually would lump them into fixed cost assumptions for ease.

```
fixed1 <- 10000
switch <- 100*30.45
fixed2 <- 5000
varOil <- 0.5
varGas <- 0.2
varWtr <- 1
prod$expense <- fixed1 + prod$oil*wi*varOil+prod$gas*wi*varGas + prod$water*wi*varWtr
prod$expense[prod$oil + prod$water < switch] <- prod$expense[prod$oil + prod$water < switch] - fixed2
```

**Taxes**

I am going to use **Texas** severance taxes, assuming an oil well, and then assume **2.5%** of **revenue **for **Ad Valorem**.

```
stxOil <- 0.046 #Texas Oil Severance
stxGas <- 0.075 #Texas Gas/NGL Severance
atx <- 0.025 #Ad Valorem
prod$tax <- prod$oilRevenue*stxOil + (prod$gasRevenue+prod$nglRevenue)*stxGas + prod$revenue*atx
```

**Cash Flow**

Simply enough, I will just calculate total cash flow. However, I am going to assume that the well is shut in on the last positive cash flow month.

```
prod$cf <- prod$revenue - prod$expense - prod$tax
last1 <- which(prod$cf > 0)
last1 <- last1[length(last1)]
prod <- prod[1:last1,]
```

**Capex/PNA**

A somewhat lazy assumption I will make is that **P&A** costs are incurred at the last month of the data frame. It is not that big of an impact from a PV-perspective that far out into the future, though you COULD defer it if it so pleases you. I also need to define capex and then how long we go from **SPUD** to **first production date**.

```
drillCapex <- 4000000 #Drill Cost
complCapex <- 6000000 #Capex + Facilities
pna <- 50000 #P&A
spudToProd <- 3 #Production occurs on month 4.
prod <- prod %>% subset(select = -c(DATE)) #Remove useless Date column
prod$capex <- 0 #Create capex column
prod$capex[nrow(prod$capex)] <- pna #Add P&A to last month
capex1 <- prod[1:spudToProd,]
capex1[1:nrow(capex1),1:length(capex1)] <- 0
capex1$capex[1] <- drillCapex
capex1$capex[spudToProd] <- complCapex
prod <- rbind(capex1, prod)
prod$monthsOn <- seq(0, nrow(prod)-1, 1) #New monthly column
```

**Cash Flow/NPV-10**

Now that we have the cash flows built out, we can have a quick peek at single-well economics. I am going to show discount rates assuming annual compounding (ie **PV = CF / (1 + r/n) ^{t*n}**). Obviously, if you change the

**discount rate**, that changes the

**PV**.

```
prod$fcf <- prod$cf - prod$capex
r <- .1
n <- 1
prod$pv <- prod$fcf/((1+r/n)^(prod$monthsOn/12*n))
library(scales)
scales::dollar(sum(prod$pv))
```

Based on my price assumption (strip as of 9/9/2020), then the **NPV-10** is **$3.7 Million US**.

**Figure 1: Cumulative Net Present Value (10%) – Example Well**

**IRR**

Instead of an **IRR package** in R, I’m going to build out an **IRR** **function** that loops quickly until we find a value, or returns 0 if it’s negative.

```
IRRcalc <- function(cf, months){ #Function inputs are our Free Cash Flow profile and our months
if(sum(cf) > 0){
IRR1 <- 3 #Start at 300%, because, come on, like it'll ever even be near this. Just starting high.
loop <- 1 #Initiate loop at 1
while(sum(cf/((1+IRR1)^(months/12))) #While with our NPV Formula
< 0){
IRR1 <- IRR1 - 0.01 #If not true, knock off 1%
loop = loop + 1 #Next Loop
}
}else {
IRR1 <- 0 #If total cf is less than 0, return 0
}
return(IRR1)
}
print(scales::percent(IRRcalc(prod$fcf, prod$monthsOn)))
```

In this scenario, I get a** 24% IRR**.

**Summary**

That’s probably enough brain-hurt for this session. If you want to test yourself, try and combine everything into a single function and see how it goes. I will be doing a pseudo-version when I publish **Part 2** of the **Shiny** **Version** later this week. Thanks for tuning in!

Pingback: Oil & Gas Coding with Python (part 3) - Shale Insights