##-----------------------------------------------------------------------------------------------# ## Date: 22-08-2023 ## Author: Giacomo Benini ## Institution: ERE, NHH - ERE, Stanford ## ## Project: "The Free lunch of Climate Change Mitigation" ## R-codes: "Dataset Analysis" ##-----------------------------------------------------------------------------------------------# rm(list = ls(all=TRUE)) options('scipen' = 100, 'digits' = 4) ##-----------------------------------------------------------------------------------------------# ## Load required packages ##-----------------------------------------------------------------------------------------------# library(tidyverse) library(ggplot2) library(fuzzyjoin) library(lubridate) ##-----------------------------------------------------------------------------------------------# ## 1. Geology and Volumes of Production ##-----------------------------------------------------------------------------------------------# panel <- read.csv('C:/Users/s15535/Desktop/2 - Energy Economics/2 - Oil & Natural Gas/2.5 - Flaring & Venting/1 - American FV&L/2 - Codes/1 - Upstream Emissions/F&V_Dataset.csv', sep = ',', header = TRUE, stringsAsFactors = FALSE, encoding = 'UTF-8') # 1.1 Time frame Start <- 2005 End <- 2020 panel <- panel %>% filter(Year >= Start & Year <= End) panel <- panel[!duplicated(panel[c('RE_ID', 'Year')]),] # No Duplicates length(unique(panel$RE_ID)) # 1464 Oil & Gas Fields length(unique(panel$Year)) # 16 Years (2005-2020) # 1.2 Balanced vs Unbalanced Data set library(BMisc) panel_b <- makeBalancedPanel(panel, idname = 'RE_ID', tname = 'Year', return_data.table = FALSE) length(unique(panel_b$RE_ID)) # 692 Fields length(unique(panel_b$RE_ID)) / length(unique(panel$RE_ID)) * 100 # 47.27% of observations are balanced # 1.3 Oil versus Gas fields panel_unique <- panel %>% select(RE_ID, Oil_Gas, Geo) %>% distinct() # 1464 table(panel_unique$Oil_Gas) # Gasfields = 373, Oilfields = 1091 table(panel_unique$Geo, panel_unique$Oil_Gas) # Gas Oil # Coalbed Methane 43 0 # Heavy & Extra Heavy 0 32 # Light & Medium 0 325 # Other Oil 0 203 # Shale & Tight Gas 330 0 # Shale & Tight Oil 0 531 # 1.4 Aggregate Production Analysis tapply(panel$Oil_B_Day, panel$Geo, summary) tapply(panel$Oil_B_Day, panel$Geo, sd) tapply(panel$TotGas_B_Day, panel$Geo, summary) tapply(panel$TotGas_B_Day, panel$Geo, sd) # 1.5 Map panel_plot_energy <- panel %>% filter(Year == 2019 & Latitude < 50) %>% # 2019 mutate(GOR_B_B = atan(TotGas_B_Day / Oil_B_Day), Energy_TB_Day = (TotGas_B_Day + Oil_B_Day) / 1000) map_og <- ggplot() map_og <- map_og + geom_polygon(data = map_data('state'), aes(x = long, y = lat, group = group), colour = 'gray', fill = 'white', alpha = 0.3, size = 0.7/.pt) map_og <- map_og + geom_point(data = panel_plot_energy, aes(x = Longitude, y = Latitude, color = GOR_B_B, size = Energy_TB_Day), alpha = 0.3) map_og <- map_og + scale_color_continuous(low = '#d1351d', high = '#204b91', na.value = 'white', name = '', limits = c(0, 1.6), breaks = c(0, 1.6), labels=c('Oil','Natural Gas')) map_og <- map_og + scale_size_continuous(range = c(0.4, 4), breaks = seq(0, 10.70, by = 2.5), name = 'Thousand BOE Day') map_og <- map_og + theme_void() map_og <- map_og + theme(legend.direction = 'horizontal', legend.position = 'bottom', legend.title = element_text(size = 8), legend.text = element_text(size = 8), text = element_text(size = 16, family = 'serif'), legend.spacing.x = unit(1.2, 'cm')) map_og <- map_og + guides(color = guide_colourbar(barwidth = 6, barheight = 0.25, title.position = 'top', title.hjust = 0.5, label.position = 'top'), size = guide_legend(title.position = 'top', title.hjust = 0.5, label.position = 'top')) map_og <- map_og + coord_map() map_og ggsave("C:/Users/s15535/Desktop/2 - Energy Economics/2 - Oil & Natural Gas/2.5 - Flaring & Venting/1 - American FV&L/0 - Manuscript/Figure1_Appendix_Map.pdf", height = 14, width = 18, units = c('cm')) rm(map_og, panel_b, panel_plot_energy, panel_unique) ##-----------------------------------------------------------------------------------------------# ## 3. Flaring Dataset ##-----------------------------------------------------------------------------------------------# # 3.1 Definition of Routine vs Non-Routine Flaring panel_gas <- panel %>% filter(Geo == 'Coalbed Methane' | Geo == 'Shale & Tight Gas') %>% mutate(Flaring_Rate = (Flared_Gas_B_Day / TotGas_B_Day), Flaring_Rate_per = Flaring_Rate*100) %>% filter(Flaring_Rate <= 0.0516289279) tapply(panel_gas$Flared_Gas_B_Day, panel_gas$Geo, summary) tapply(panel_gas$Flared_Gas_B_Day, panel_gas$Geo, sd) quantile(panel_gas$Flaring_Rate_per, probs = seq(0, 1, 0.05), na.rm = TRUE, names = TRUE) mean(panel_gas$Flaring_Rate_per) library(RColorBrewer) nb.cols <- 6 # Define the number of colors you want mycolors <- colorRampPalette(brewer.pal(11, "RdBu"))(nb.cols) dens_flare <- ggplot(data = panel_gas) dens_flare <- dens_flare + geom_density(aes(x = Flaring_Rate_per, color = Geo), size = 1) dens_flare <- dens_flare + theme_minimal() dens_flare <- dens_flare + xlim(-0.01, 1) dens_flare <- dens_flare + geom_vline(xintercept = 0.054093745, linetype = 'dotted', color = 'blue', size = 1) dens_flare <- dens_flare + scale_color_manual(values = mycolors) dens_flare <- dens_flare + labs(x = 'Percentage of Natural Gas Flared', y = 'Empirical Density') dens_flare <- dens_flare + theme(legend.direction = "horizontal", legend.position = "bottom", axis.text.x = element_text(angle = 90, hjust = 1)) dens_flare <- dens_flare + guides(fill = guide_legend(title.position = 'top', title.hjust = 0.5, label.position = 'top')) dens_flare ggsave("C:/Users/s15535/Desktop/2 - Energy Economics/2 - Oil & Natural Gas/2.5 - Flaring & Venting/1 - American FV&L/0 - Manuscript/Figure2_Appendix_GasfieldsFlareDensity.pdf", height = 15, width = 20, units = c('cm')) # 3.2 Flaring Dataset panel_oil <- panel %>% filter(Geo == 'Heavy & Extra Heavy' | Geo == 'Light & Medium' | Geo == 'Other Oil' | Geo == 'Shale & Tight Oil') panel_oil <- panel_oil %>% filter(GOR_SCF_B < 100000) # Threshold defined by the EIA definition of an oilfield panel_prob <- panel %>% filter(Geo == 'Heavy & Extra Heavy' | Geo == 'Light & Medium' | Geo == 'Other Oil' | Geo == 'Shale & Tight Oil') %>% filter(GOR_SCF_B > 100000) length(unique(panel_prob$RE_ID)) # 287 Fields rm(panel_prob) panel_oil <- panel_oil %>% filter(!is.na(Flared_Gas_B_Day)) %>% mutate(Non_Routine_Flaring_B_Day = ifelse(quantile(panel_gas$Flaring_Rate, probs = 0.5, na.rm = TRUE) * TotGas_B_Day <= Flared_Gas_B_Day, quantile(panel_gas$Flaring_Rate, probs = 0.5, na.rm = TRUE) * TotGas_B_Day, Flared_Gas_B_Day), Routine_Flaring_B_Day = ifelse(Flared_Gas_B_Day - Non_Routine_Flaring_B_Day <= 0, 0, Flared_Gas_B_Day - Non_Routine_Flaring_B_Day), Dummy = ifelse(Routine_Flaring_B_Day == 0, 0, 1)) %>% distinct() length(unique(panel_oil$RE_ID)) # 627 Fields panel_prob <- panel_oil %>% filter(GOR_SCF_B > 100000) # EIA definition of an oilfield length(unique(panel_oil$RE_ID)) # 599 Fields rm(dens_flare) tapply(panel_oil$Oil_B_Day, panel_oil$Geo, summary) tapply(panel_oil$Oil_B_Day, panel_oil$Geo, sd) tapply(panel_oil$TotGas_B_Day, panel_oil$Geo, summary) tapply(panel_oil$TotGas_B_Day, panel_oil$Geo, sd) ##-----------------------------------------------------------------------------------------------# ## 2. Production Analysis -> Rystad vs EIA data coverage ##-----------------------------------------------------------------------------------------------# # 2.1 Dataset vs EIA Production Data agg_prod <- panel %>% mutate(Year = as.numeric(Year)) %>% group_by(Year) %>% summarise(Oil_MMB_Day = sum(Oil_B_Day, na.rm = TRUE) / 1000000, Gas_MMB_Day = sum(TotGas_B_Day, na.rm = TRUE) / 1000000) %>% ungroup() agg_prod geo_prod <- panel %>% mutate(Year = as.numeric(Year)) %>% group_by(Geo, Year) %>% summarise(Oil_MMB_Day_Geo = sum(Oil_B_Day, na.rm = TRUE) / 1000000, Gas_MMB_Day_Geo = sum(TotGas_B_Day, na.rm = TRUE) / 1000000) %>% ungroup() %>% left_join(., agg_prod) %>% group_by(Geo) %>% summarise(Oil_share = mean(Oil_MMB_Day_Geo / Oil_MMB_Day * 100), Gas_share = mean(Gas_MMB_Day_Geo / Gas_MMB_Day * 100)) %>% ungroup() geo_prod # Geo Oil_share Gas_share # # 1 Coalbed Methane 0.0829 2.29 # 2 Heavy & Extra Heavy 8.36 0.232 # 3 Light & Medium 36.3 28.3 # 4 Other Oil 9.43 9.86 # 5 Shale & Tight Gas 5.57 47.5 # 6 Shale & Tight Oil 40.3 11.8 # 2.2 EIA Production Data # EIA Oil Production (30/09/2021) # https://www.eia.gov/dnav/pet/pet_crd_crpdn_adc_mbblpd_a.htm eia_oil <- read.csv('C:/Users/s15535/Desktop/2 - Energy Economics/2 - Oil & Natural Gas/2.5 - Flaring & Venting/1 - American FV&L/2 - Codes/1 - Upstream Emissions/Aggregate US Production/US_Oil_Production.csv', sep = ',', header = TRUE, stringsAsFactors = FALSE, encoding = 'UTF-8') # double check of the totals eia_oil_check <- eia_oil %>% group_by(Date) %>% mutate(Tot_TB_Day = mean(U.S..Field.Production.of.Crude.Oil..Thousand.Barrels.per.Day., na.rm = TRUE), Sum_TB_Day = sum(East.Coast..PADD.1..Field.Production.of.Crude.Oil..Thousand.Barrels.per.Day., Midwest..PADD.2..Field.Production.of.Crude.Oil..Thousand.Barrels.per.Day., Gulf.Coast..PADD.3..Field.Production.of.Crude.Oil..Thousand.Barrels.per.Day., Rocky.Mountain..PADD.4..Field.Production.of.Crude.Oil..Thousand.Barrels.per.Day., West.Coast..PADD.5..Field.Production.of.Crude.Oil..Thousand.Barrels.per.Day., na.rm = TRUE)) %>% select(Date, Tot_TB_Day, Sum_TB_Day) # aggregation eia_oil <- eia_oil %>% mutate(EIA_Tot_Oil_MMB_Day = U.S..Field.Production.of.Crude.Oil..Thousand.Barrels.per.Day. / 1000, EIA_Offshore_Oil_MMB_Day = (Federal.Offshore..Gulf.of.Mexico.Field.Production.of.Crude.Oil..Thousand.Barrels.per.Day. / 1000) + (Federal.Offshore.PADD.5.Field.Production.of.Crude.Oil..Thousand.Barrels.per.Day./ 1000), EIA_Onshore_Oil_MMB_Day = EIA_Tot_Oil_MMB_Day - EIA_Offshore_Oil_MMB_Day) %>% rename(Year = Date) %>% select(Year, EIA_Tot_Oil_MMB_Day, EIA_Onshore_Oil_MMB_Day, EIA_Offshore_Oil_MMB_Day) # EIA Gas Production (30/09/2021) # https://www.eia.gov/dnav/ng/ng_prod_sum_a_EPG0_FGW_mmcf_a.htm eia_gas <- read.csv('C:/Users/s15535/Desktop/2 - Energy Economics/2 - Oil & Natural Gas/2.5 - Flaring & Venting/1 - American FV&L/2 - Codes/1 - Upstream Emissions/Aggregate US Production/US_Gas_Production.csv', sep = ',', header = TRUE, stringsAsFactors = FALSE, encoding = 'UTF-8') # double check of the totals eia_gas_check <- eia_gas %>% rowwise() %>% mutate(Sum_Off = sum(US..State.Offshore.Natural.Gas.Gross.Withdrawals..MMcf., US..Federal.Offshore.Natural.Gas.Gross.Withdrawals..MMcf., na.rm = TRUE)) %>% ungroup() %>% select(Date, U.S..Natural.Gas.Gross.Withdrawals..MMcf., U.S..Natural.Gas.Gross.Withdrawals.Offshore..MMcf., Sum_Off) # aggregation eia_gas <- eia_gas %>% mutate(EIA_Tot_Gas_MMB_Day = (U.S..Natural.Gas.Gross.Withdrawals..MMcf. / 6000) / 365, EIA_Offshore_Gas_MMB_Day = (U.S..Natural.Gas.Gross.Withdrawals.Offshore..MMcf. / 6000) / 365, EIA_Onshore_Gas_MMB_Day = EIA_Tot_Gas_MMB_Day - EIA_Offshore_Gas_MMB_Day) %>% rename(Year = Date) %>% select(Year, EIA_Tot_Gas_MMB_Day, EIA_Onshore_Gas_MMB_Day, EIA_Offshore_Gas_MMB_Day) # US Coverage Panel Fit coverage <- left_join(agg_prod, eia_oil) %>% left_join(., eia_gas) %>% group_by(Year) %>% summarise(Oil_Coverage = Oil_MMB_Day / EIA_Onshore_Oil_MMB_Day * 100, Gas_Coverage = Gas_MMB_Day / EIA_Onshore_Gas_MMB_Day * 100) %>% ungroup() coverage mean(coverage$Oil_Coverage) # 75.72 % mean(coverage$Gas_Coverage) # 77.93 % rm(agg_prod, coverage, eia_gas_check, geo_prod, eia_oil_check, panel_unique) ##-----------------------------------------------------------------------------------------------# ## 4. Analysis of the Flaring Dataset ##-----------------------------------------------------------------------------------------------# # 4.1 Flaring Dataset vs EIA Production Data agg_prod <- panel_oil %>% mutate(Year = as.numeric(Year)) %>% group_by(Year) %>% summarise(Oil_MMB_Day = sum(Oil_B_Day, na.rm = TRUE) / 1000000, Gas_MMB_Day = sum(TotGas_B_Day, na.rm = TRUE) / 1000000) %>% ungroup() agg_prod geo_prod <- panel_oil %>% mutate(Year = as.numeric(Year)) %>% group_by(Geo, Year) %>% summarise(Oil_MMB_Day_Geo = sum(Oil_B_Day, na.rm = TRUE) / 1000000, Gas_MMB_Day_Geo = sum(TotGas_B_Day, na.rm = TRUE) / 1000000) %>% ungroup() %>% left_join(., agg_prod) %>% group_by(Geo) %>% summarise(Oil_share = mean(Oil_MMB_Day_Geo / Oil_MMB_Day * 100), Gas_share = mean(Gas_MMB_Day_Geo / Gas_MMB_Day * 100)) %>% ungroup() geo_prod # Geo Oil_share Gas_share # * # 1 Coalbed Methane 0.129 1.39 # 2 Heavy & Extra Heavy 0.0521 0.000214 # 3 Light & Medium 37.0 19.7 # 4 Other Oil 6.57 5.90 # 5 Shale & Tight Gas 14.4 65.4 # 6 Shale & Tight Oil 41.9 7.57 coverage <- left_join(agg_prod, eia_oil) %>% left_join(., eia_gas) %>% group_by(Year) %>% summarise(Oil_Coverage = Oil_MMB_Day / EIA_Onshore_Oil_MMB_Day * 100, Gas_Coverage = Gas_MMB_Day / EIA_Onshore_Gas_MMB_Day * 100) %>% ungroup() coverage mean(coverage$Oil_Coverage) # 39.42 % mean(coverage$Gas_Coverage) # 12.56 % rm(agg_prod, coverage, eia_gas, eia_oil, geo_prod, panel_oil_unique) # 4.2 Flared Gas BOE/Day library(ggpubr) dens_non_routine_flare <- ggplot(data = panel_oil) dens_non_routine_flare <- dens_non_routine_flare + geom_density(aes(x = Non_Routine_Flaring_B_Day, color = Geo), size = 1) dens_non_routine_flare <- dens_non_routine_flare + theme_minimal() dens_non_routine_flare <- dens_non_routine_flare + xlim(-1, 10) dens_non_routine_flare <- dens_non_routine_flare + scale_color_manual(values = mycolors) dens_non_routine_flare <- dens_non_routine_flare + labs(x = 'Non-Routine Flaring (BOE / Day)', y = 'Empirical Density') dens_non_routine_flare <- dens_non_routine_flare + theme(legend.direction = "horizontal", legend.position = "bottom", axis.text.x = element_text(angle = 90, hjust = 1)) dens_non_routine_flare <- dens_non_routine_flare + guides(fill = guide_legend(title.position = 'top', title.hjust = 0.5, label.position = 'top')) dens_non_routine_flare dens_flare <- ggplot(data = panel_oil) dens_flare <- dens_flare + geom_density(aes(x = Routine_Flaring_B_Day, color = Geo), size = 1) dens_flare <- dens_flare + theme_minimal() dens_flare <- dens_flare + xlim(-10, 100) dens_flare <- dens_flare + scale_color_manual(values = mycolors) dens_flare <- dens_flare + labs(x = 'Routine Flaring (BOE / Day)', y = 'Empirical Density') dens_flare <- dens_flare + theme(legend.direction = "horizontal", legend.position = "bottom", axis.text.x = element_text(angle = 90, hjust = 1)) dens_flare <- dens_flare + guides(fill = guide_legend(title.position = 'top', title.hjust = 0.5, label.position = 'top')) dens_flare figure <- ggarrange(dens_non_routine_flare, dens_flare, # labels = c("Non-Routine Flaring", "Routine Flaring"), common.legend = TRUE, legend = "bottom", ncol = 2, nrow = 1) figure summary(panel_oil$Non_Routine_Flaring_B_Day) sd(panel_oil$Non_Routine_Flaring_B_Day) tapply(panel_oil$Non_Routine_Flaring_B_Day, panel_oil$Geo, summary) tapply(panel_oil$Non_Routine_Flaring_B_Day, panel_oil$Geo, sd) tapply(panel_oil$Routine_Flaring_B_Day, panel_oil$Geo, summary) tapply(panel_oil$Routine_Flaring_B_Day, panel_oil$Geo, sd) ggsave("C:/Users/s15535/Desktop/2 - Energy Economics/2 - Oil & Natural Gas/2.5 - Flaring & Venting/1 - American FV&L/0 - Manuscript/Figure3_Appendix_RoutineFlareDensity.pdf", height = 15, width = 20, units = c('cm')) ##-----------------------------------------------------------------------------------------------# ## 4. Plot for Second Stage Estimation ##-----------------------------------------------------------------------------------------------# library(RColorBrewer) nb.cols <- 6 # Define the number of colors you want mycolors <- colorRampPalette(brewer.pal(11, "RdBu"))(nb.cols) # 4.1 Other Gas BOE/Day dens_other <- ggplot(data = panel_oil) dens_other <- dens_other + geom_density(aes(x = Other_Gas_B_Day, color = Geo), size = 1) dens_other <- dens_other + theme_minimal() dens_other <- dens_other + xlim(-1, 100) dens_other <- dens_other + scale_color_manual(values = mycolors) dens_other <- dens_other + labs(x = 'Other (BOE / Day)', y = 'Empirical Density') dens_other <- dens_other + theme(legend.direction = "horizontal", legend.position = "bottom", axis.text.x = element_text(angle = 90, hjust = 1)) dens_other <- dens_other + guides(fill = guide_legend(title.position = 'top', title.hjust = 0.5, label.position = 'top')) dens_other # 4.2 Flared Gas BOE/Day dens_flare <- ggplot(data = panel_oil) dens_flare <- dens_flare + geom_density(aes(x = Flared_Gas_B_Day, color = Geo), size = 1) dens_flare <- dens_flare + theme_minimal() dens_flare <- dens_flare + xlim(-1, 100) dens_flare <- dens_flare + scale_color_manual(values = mycolors) dens_flare <- dens_flare + labs(x = 'Flaring (BOE / Day)', y = 'Empirical Density') dens_flare <- dens_flare + theme(legend.direction = "horizontal", legend.position = "bottom", axis.text.x = element_text(angle = 90, hjust = 1)) dens_flare <- dens_flare + guides(fill = guide_legend(title.position = 'top', title.hjust = 0.5, label.position = 'top')) dens_flare # 4.3 Gas Price $/BOE flare_price <- ggplot(panel_oil, aes(x = R_HH_Dollars_B, y = Flared_Gas_B_Day)) flare_price <- flare_price + geom_point() flare_price <- flare_price + geom_smooth(method = 'lm') flare_price <- flare_price + theme_minimal() flare_price <- flare_price + labs(x = 'Flared Gas (B / Day)', y = 'Real Natural Gas Price ($ / BOE)') flare_price <- flare_price + theme(legend.direction = "horizontal", legend.position = "bottom", axis.text.x = element_text(angle = 90, hjust = 1)) flare_price <- flare_price + guides(fill = guide_legend(title.position = 'top', title.hjust = 0.5, label.position = 'top')) flare_price # 4.4 Insitu Gas BOE/Day dens_insitu <- ggplot(data = panel_oil) dens_insitu <- dens_insitu + geom_density(aes(x = Insitu_Gas_B_Day, color = Geo), size = 1) dens_insitu <- dens_insitu + theme_minimal() dens_insitu <- dens_insitu + xlim(-1, 100) dens_insitu <- dens_insitu + scale_color_manual(values = mycolors) dens_insitu <- dens_insitu + labs(x = 'Insitu (BOE / Day)', y = 'Empirical Density') dens_insitu <- dens_insitu + theme(legend.direction = "horizontal", legend.position = "bottom", axis.text.x = element_text(angle = 90, hjust = 1)) dens_insitu <- dens_insitu + guides(fill = guide_legend(title.position = 'top', title.hjust = 0.5, label.position = 'top')) dens_insitu # 4.5 Injection Gas BOE/Day tapply(panel_oil$Re_Gas_B_Day, panel_oil$Geo, summary) dens_re <- ggplot(data = panel_oil) dens_re <- dens_re + geom_density(aes(x = Re_Gas_B_Day, color = Geo), size = 1) dens_re <- dens_re + theme_minimal() dens_re <- dens_re + xlim(-1, 1) dens_re <- dens_re + scale_color_manual(values = mycolors) dens_re <- dens_re + labs(x = 'Injections (BOE / Day)', y = 'Empirical Density') dens_re <- dens_re + theme(legend.direction = "horizontal", legend.position = "bottom", axis.text.x = element_text(angle = 90, hjust = 1)) dens_re <- dens_re + guides(fill = guide_legend(title.position = 'top', title.hjust = 0.5, label.position = 'top')) dens_re panel_noshale <- panel_oil %>% filter(Geo != 'Shale & Tight Oil') %>% filter(Re_Gas_B_Day > 0) %>% distinct() dens_re <- ggplot(data = panel_noshale) dens_re <- dens_re + geom_density(aes(x = Re_Gas_B_Day, color = Geo), size = 1) dens_re <- dens_re + theme_minimal() dens_re <- dens_re + xlim(-1, 100) dens_re <- dens_re + scale_color_manual(values = mycolors) dens_re <- dens_re + labs(x = 'Injections (BOE / Day)', y = 'Empirical Density') dens_re <- dens_re + theme(legend.direction = "horizontal", legend.position = "bottom", axis.text.x = element_text(angle = 90, hjust = 1)) dens_re <- dens_re + guides(fill = guide_legend(title.position = 'top', title.hjust = 0.5, label.position = 'top')) dens_re rm(dens_flare, dens_insitu, dens_re, dens_other, flare_price, mycolors, nb.cols, other_insitu, other_re, wells_re) agg_prod <- panel_oil %>% mutate(Year = as.numeric(Year)) %>% group_by(Year) %>% summarise(Oil_MMB_Day = sum(Oil_B_Day, na.rm = TRUE) / 1000000, Gas_MMB_Day = sum(TotGas_B_Day, na.rm = TRUE) / 1000000) %>% ungroup() agg_prod geo_prod_per <- panel_oil %>% mutate(Year = as.numeric(Year)) %>% group_by(Geo, Year) %>% summarise(Oil_MMB_Day_Geo = sum(Oil_B_Day, na.rm = TRUE) / 1000000, Gas_MMB_Day_Geo = sum(TotGas_B_Day, na.rm = TRUE) / 1000000) %>% ungroup() %>% left_join(., agg_prod) %>% group_by(Geo) %>% summarise(Oil_share = mean(Oil_MMB_Day_Geo / Oil_MMB_Day * 100), Gas_share = mean(Gas_MMB_Day_Geo / Gas_MMB_Day * 100)) %>% ungroup() geo_prod_per # Geo Oil_share Gas_share # # 1 Heavy & Extra Heavy 8.78 0.447 # 2 Light & Medium 38.1 53.6 # 3 Other Oil 9.94 19.5 # 4 Shale & Tight Oil 43.2 26.4 geo_prod_vol <- panel_oil %>% mutate(Year = as.numeric(Year)) %>% group_by(Geo, Year) %>% summarise(Oil_MMB_Day_Geo = sum(Oil_B_Day, na.rm = TRUE) / 1000000, Gas_MMB_Day_Geo = sum(TotGas_B_Day, na.rm = TRUE) / 1000000) %>% ungroup() %>% left_join(., agg_prod) %>% group_by(Geo) %>% summarise(Oil = mean(Oil_MMB_Day_Geo), Gas = mean(Gas_MMB_Day_Geo)) %>% ungroup() geo_prod_vol rm(agg_prod, coverage, eia_gas, eia_oil, geo_prod, panel_oil_unique) rm(geo_prod_per, geo_prod_vol, panel_unique, End, Start) rm(panel_b)