R. Попутный газ в Северной Дакоте

Вопрос дня: насколько успешно удается реализовывать попутный газ в месторождениях Северной Дакоты? Ниже я попытаюсь ответить на этот вопрос при помощи R, как всегда по шагам.

1. Загрузка библиотек

library(dplyr)
library(ggplot2)
library(curl)
library(stringr)
library(stringi)
library(data.table)
library(xts)
library(tidyr)
library(forecast) #ggplot for ts class
library(seasonal) #X11 seasonal decomposition

2. Загрузка данных из уже сформированной мной базы данных (БД).

Формат этой БД соответствует формату ежемесячных отчетов о добыче нефти и газа, публикуемых правительством Северной Дакоты в виде XLXS-файлов (https://www.dmr.nd.gov/oilgas/mprindex.asp), по-сути эта БД и представляет собой объедение этих Excel-файлов и аналогичных таблиц в PDF, относящихся к периоду до 2015 года. Для удобства сразу же преобразую таблицу в формат data.table

load("odb_mm_updated.RData")
odb_mm_dt <- data.table(odb_mm_updated)

3. Расширяем таблицу вычисляемыми полями и проводим фильтрацию данных.

Для каждой скважины сформируем уникальный идентификатор из некоторого набора описательных полей.

odb_mm_dt[ , u2 := paste(FileNo, County, Pool, Quarter, Section, Township, Range, WellName)]
setindex( odb_mm_dt,u2 )

Для каждой скважины вычисляем дату начала ее функционирования (точнее первый месяц предоставления отчетности).

##Минимальная дата
x1 <- odb_mm_dt[ , .(minDate = min(ReportDate)), by=.(u2)]
setindex(x1, u2)
odb_mm_dt$minDate <- x1[ .(odb_mm_dt$u2), minDate, on="u2"]

Данный датасет содержит некоторые проблемы, а именно: на один идентификатор могут попадать данные нескольких скважин. Подобные задвоения выявляются подсчетом числа строк для данного идентификатора и числа месяцев наблюдений, для незадвоенных данных эти числа должны совпадать.

u1 <- unique( odb_mm_dt$u2 )
print("Исходное число задвоений")

## [1] "Исходное число задвоений"

odb_mm_dt[ , .(n=.N), by=.(u2)][ , max(n)]

## [1] 132

print("Максимальное расхождение задвоений")

## [1] "Максимальное расхождение задвоений"

odb_mm_dt[ , .(ndoubles=(length(ReportDate)- length(unique(ReportDate)))), by=.(u2)][ , max(ndoubles)]

## [1] 4

#4
odb_mm_dt[ , .(ndoubles=(length(ReportDate)- length(unique(ReportDate)))), by=.(u2)][ ndoubles > 0 , .N ]

## [1] 45

#45

Задвоений оказывается немного на фоне десятков тысяч скважин - и можно удалить без ущерба для статистики. Большинство таких задвоений относится к скважинам, отмеченным как конфиденциальные, потому убираем их, а потом и другие задвоения.

w1 <- odb_mm_dt[ !is.na(minDate) & (Pool !="CONFIDENTIAL") & !is.na(Pool) &
                   (!is.na(Gas) & !is.na(GasSold)) & (minDate >="2009-01-01")]

print( "Number of doubles")

## [1] "Number of doubles"

w1[ , .(ndoubles=(length(ReportDate)- length(unique(ReportDate)))), by=.(u2)][ ndoubles > 0 , .N ]

## [1] 9

i1 <- w1[ , .(ndoubles=(length(ReportDate)- length(unique(ReportDate)))), by=.(u2)][ ndoubles > 0 ]
#отфильтровываем их 
w1 <- w1[ !( u2 %in% i1$u2)   ]

u1 <- unique(w1$u)
print( "Number of doubles")

## [1] "Number of doubles"

w1[ , .(ndoubles=(length(ReportDate)- length(unique(ReportDate)))), by=.(u2)][ ndoubles > 0 , .N ]

## [1] 0

save(w1, file="w1.RData")

4. Начинаем построение разных графиков. Для начала посмотрим как росла добыча газа и его продажи.

С 2009 года они увеличились почти в 100 раз! При этом доля продаж газа подросла с 65% до примерно 75-85% - отличный показатель утилизации.

x1 <- w1[  , .(gas_produced = sum(Gas), gas_sold=sum(GasSold)), by =.(ReportDate)]

x2 <- gather(x1, key="Type",value="V", "gas_sold", "gas_produced")

ggplot(x2, aes(x=ReportDate, y=log10(V), color=Type))+geom_point()+ggtitle("Gas sold and produces") + ylab("Log10(Volume,mcf/y)")

ggplot(x1, aes(x=ReportDate, y=gas_sold/gas_produced))+geom_path(color="Red")+ggtitle("Share of sold gas")

Доля продаваемого газа варьируется скорее год от года, чем в зависимости от общего объема добычи.

ggplot(x1, aes(x=gas_produced, y=gas_sold/gas_produced, color=as.factor(year(ReportDate))))+
  geom_point()+
  scale_color_discrete()+
  geom_smooth(aes(x=gas_produced, y=gas_sold/gas_produced, color="")) +
  ggtitle("Share of sold gas vs production volume")+
  labs(color="Report Year") +xlab("Gas produced, mcf/y")

5. В разрезе провинций (графств) добыча газа выглядит так. Есть четыре лидера, но больше всех добывает графство McKenzie.

d1 <- w1[year(ReportDate)==2018, .(gas_produced = sum(Gas), gas_sold = sum(GasSold)), by=County][order(gas_produced)]

d2 <- bind_rows( d1[ gas_produced >= 0.05*sum(gas_produced) ], d1[ gas_produced < 0.05*sum(gas_produced), .(County="Others", gas_produced=sum(gas_produced), gas_sold=sum(gas_sold)) ] )

ggplot(d2, aes(x="", fill=factor(County, levels = County[order(gas_produced)]),y=gas_produced))+
  geom_bar(stat="identity")+
  geom_text(stat = "identity", aes(y=gas_produced,label=paste0( County, "\n(",round(100*gas_produced/sum(gas_produced)),"%)")), hjust=0.5, vjust=0.5,position= position_stack(0.5), size=3)+ 
  coord_polar("y", start=0)+
  ggtitle( "North Dakota gas production by county, 2018") + labs(fill = "ND County code")+xlab("Volume, mcf/y")

Аналогичная картинка будет и относительно продаваемого газа. В целом все продают примерно одну и ту же долю своей добычи.

ggplot(w1[year(ReportDate)==2018, .(gas_produced= sum(Gas), gas_sold = sum(GasSold)), by=County], aes(x=County,y=gas_produced, fill="1"))+
  geom_bar(stat="identity", fill="red", show.legend=TRUE)+
  geom_bar(aes(y=gas_sold, fill="2"), stat="identity", fill="blue")+
  ylab("Volume produced, sold, mcf/y")+
  ggtitle("ND Gas produced, sold by county")

Провинция, которая производит больше нефти, производит и больше газа.

ggplot(w1[year(ReportDate)==2018, .(gs= sum(Gas, na.rm=TRUE), go = sum(Oil, na.rm=TRUE)), by=County][order(gs)], aes(x=go,y=gs))+
  geom_point(color="red", size=3)+
  geom_text(aes(label=County), hjust=1)+
  xlab("County oil poductionб b/y")+
  ylab("County gas production, mcf/y")+
  ggtitle("Gas vs. Oil production by ND county, 2018")

6. Выясним как соотносится добыча газа и его реализация по скважинам - для какого-нибудь года. Как меняется доля проданного газа в зависимости от объема добычи скважины?

x2018 <- w1[ (year(ReportDate)==2018)  , .(gas_produced = sum(Gas), gas_sold=sum(GasSold), FileNo), by =.(u2)][ (gas_produced > 0) ]

##Без групп скважин 
ggplot(x2018[ FileNo != "" ], aes(x=gas_produced, y=gas_sold))+geom_point(color="red")+xlab("Gas produced, mcf/y")+ylab("Gas Sold, mcf/y")+ggtitle("Gas produced and sold by wells in 2018")

Из приведенной картинки примерно понятно, что зависимости доли утилизации от объема производства нет, все продают примерно одинаково, но еще лучше это видно на картинке, отражающей статистику при их группировке скважин по объемам добычи (с шагом 50000 mcf/y).

ggplot(x2018[ FileNo != "" ], aes(x=as.factor( floor(gas_produced/5e4)), y=gas_sold/gas_produced))+geom_boxplot()+xlab("Gas produced by well /50000mcf/y")+ylab("Share of gas sold")+ggtitle("Share of gas sold vs. gas produced by wells in 2018")

Как меняется доля проданного газа в зависимости от года открытия скважины? В 2018г. у основной массы (более 75%) скважин доля проданного газа была более 50%, хотя вариативность тут достаточно сильная. Чем новее скважины тем в среднем они лучше утилизируют газ.

w2 <- w1[ (year(ReportDate)==2018), .(gas_produced = sum(Gas), gas_sold=sum(GasSold)), by =.(u2, ComplYear = year(minDate))][ (gas_produced > 0) ]

ggplot(w2, aes(x=as.factor(ComplYear),y=gas_sold/gas_produced))+
  geom_boxplot()+
  xlab("Year of completion")+
  ylab("Gas sold/Gas produced") + 
  ggtitle("Share of sold gas distribution (in 2018)")

ggplot(w2, aes(x=as.factor(ComplYear),y=gas_sold/gas_produced) )+
  geom_violin()+
  xlab("Year of completion")+
  ylab("Gas sold/Gas produced") + 
  ggtitle("Share of sold gas distribution (in 2018)")

В информационном плане эквивалентная картинка, но показывающая кумулятивную функцию распределения скважин по доле проданного газа. Опять-таки видим,что более 75% новых скважин продают более 75% добываемого каждой из них газа.

w3ecdf <- w2[  ,.( xx=seq(0,1,by=0.05), ecdf1 = ecdf(1-gas_sold/gas_produced)(seq(0,1,by=0.05))), by=.(ComplYear) ]

tt <- "Year of completeon"
ggplot(w3ecdf, aes(x=xx,y=ecdf1, color=as.factor(ComplYear), linetype=as.factor(ComplYear)) )+geom_line()+ xlab("Share of unsold gas")+ylab("Share of number of wells (unsold<x)")+labs(color=tt, linetype=tt) +
  ggtitle("Distribution of number of wells in respect to share of unsold gas\n(in 2018)")

Как меняется доля проданного газа в зависимости от текущего года? Попутный газ год от года утилизируется все лучше, в 2009 г. очень многие скважины этот газ вообще почти не продавали.

w3 <- w1[ , .(gas_produced = sum(Gas), gas_sold=sum(GasSold)), by =.(u2, ReportYear= year(ReportDate))][ (gas_produced > 0) ]

ggplot(w3, aes(x=as.factor(ReportYear),y=gas_sold/gas_produced))+
  geom_boxplot()+
  xlab("Report Year")+
  ylab("Gas sold/Gas produced") + 
  ggtitle("Share of sold gas distribution")

ggplot(w3, aes(x=as.factor(ReportYear),y=gas_sold/gas_produced) )+
  geom_violin()+
  xlab("Report Year")+
  ylab("Gas sold/Gas produced") + 
  ggtitle("Share of sold gas distribution")

7. Посмотрим сезонность добычи и реализации.

xts1 <- xts(log10(x1$gas_sold/x1$gas_produced), order.by = x1$ReportDate, frequency = 12)
xts1 <- xts1["2010-01-01/2019-01-01"]
ts1 <- ts(xts1, frequency=12, start=c(year(index(xts1)[1]), month(index(xts1)[1])))
ts1d <-decompose( ts1, type="add" ) 
autoplot( ts1d) + ggtitle("Decomposition of 'share of sold gas' series" )

xts1 <- xts(log(x1$gas_sold), order.by = x1$ReportDate, frequency=12)
xts1 <- xts1["2010-01-01/2019-01-01"]
ts1 <- ts(xts1, frequency=12, start=c(year(index(xts1)[1]), month(index(xts1)[1])))
ts1d <- decompose( ts1, type="additive" )
autoplot(ts1d) + ggtitle("Decomposition of 'log(Volume of sold gas, mcf/m') series" )

ts1d$seasonal[1:12]

##  [1] -0.026525928 -0.085062099  0.019801519 -0.008962222  0.032513311
##  [6] -0.004718137  0.025615591  0.022801115  0.001729001  0.039650526
## [11]  0.002933995 -0.019776673

На фоне прочей вариативности сезонные колебания доли продаваемого газа, да и его объема, оказываются малозаметными.

Для эксперимента пробуем теперь более продвинутый метод сезонной декомпозиции, используемый правительством США, - X11.

s1 <- seas( ts1, transform.function = "none")
summary(s1)

## 
## Call:
## seas(x = ts1, transform.function = "none")
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## Constant        0.03599    0.01268   2.839  0.00452 ** 
## LS2013.Dec     -0.20289    0.04175  -4.859 1.18e-06 ***
## LS2016.Dec     -0.22356    0.04198  -5.326 1.00e-07 ***
## AR-Seasonal-12  0.91485    0.03987  22.943  < 2e-16 ***
## MA-Seasonal-12  0.56344    0.07794   7.229 4.87e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## SEATS adj.  ARIMA: (0 1 0)(1 0 1)  Obs.: 108  Transform: none
## AICc: -337.6, BIC: -322.4  QS (no seasonality in final):    0  
## Box-Ljung (no autocorr.): 26.51   Shapiro (normality):  0.99

autoplot(s1) + ggtitle("X11 seasonal decomposition of 'log(volume of sold gas)' series" )

На качественном уровне результат тот же. Сезонные колебания “тонут” на фоне прочих.