##-----------------------------------------------------------------------------------------------# ## Date: 16-09-2024 ## Author: Giacomo Benini ## Institution: ERE, NHH - ERE, Stanford ## ## Project: "Flaring & Venting" ## R-codes: "Results Analysis" ##-----------------------------------------------------------------------------------------------# rm(list = ls(all=TRUE)) options('scipen' = 100, 'digits' = 4) ##-----------------------------------------------------------------------------------------------# ## Load required packages ##-----------------------------------------------------------------------------------------------# library(tidyverse) library(ggplot2) library(fuzzyjoin) library(scales) library(viridis) library(usmap) library(maps) library(fmsb) library(patchwork) devtools::install_github("ricardo-bion/ggradar", dependencies = TRUE) library(ggradar) library(scales) #library(ggvanced) library(ggpubr) ##-----------------------------------------------------------------------------------------------# ## 1. Balanced vs Unbalanced F&V Results ##-----------------------------------------------------------------------------------------------# panel_fv <- 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/Second_Stage_Results.csv', sep = ',', header = TRUE, stringsAsFactors = FALSE, encoding = 'UTF-8') # 1.1 Balance versus Unbalance panel_fv <- panel_fv[!duplicated(panel_fv[c('RE_ID', 'Year')]),] # No Duplicates length(unique(panel_fv$RE_ID)) # 556 Oilfields length(unique(panel_fv$Year)) # 16 Years (2005-2020) balance_panel_fv <- panel_fv %>% group_by(RE_ID) %>% filter(n_distinct(Year) == n_distinct(.$Year)) %>% ungroup length(unique(balance_panel_fv$RE_ID)) # 13 Oilfields length(unique(balance_panel_fv$Year)) # 16 Years (2005-2020) # Balance dataset: 574 * 16 = 9184 # Actual dataset: 4940 / 9184 * 100 = 53.79% # Fraction of balanced oilfields: 13/574*100 = 2.26% # 1.2 Summary Statistics by Geological Type panel_unique <- panel_fv %>% select(RE_ID, Oil_Gas, Geo) %>% distinct() # 1476 table(panel_unique$Oil_Gas) # Gasfields = 0, Oilfields = 574 table(panel_unique$Geo, panel_unique$Oil_Gas) # Gas Oil # Heavy & Extra Heavy 0 4 # Light & Medium 0 145 # Other Oil 0 17 # Shale & Tight Oil 0 408 ##-----------------------------------------------------------------------------------------------# ## 2. Data Coverage -> Rystad (non NA Flaring) vs EIA data coverage ##-----------------------------------------------------------------------------------------------# # 2.1 Dataset vs EIA Production Data agg_fv <- panel_fv %>% 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_fv geo_prod <- panel_fv %>% 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_fv) %>% 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 Heavy & Extra Heavy 0.386 0.0681 # 2 Light & Medium 39.3 44.0 # 3 Other Oil 2.90 9.10 # 4 Shale & Tight Oil 60.3 50.2 # 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') # doble 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') # doble 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_fv, 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) # 36.88% mean(coverage$Gas_Coverage) # 11.70% rm(list = ls(all=TRUE)) ##-----------------------------------------------------------------------------------------------# ## 3. Energy Savings 2005-2020 ##-----------------------------------------------------------------------------------------------# panel_fv <- 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/Second_Stage_Results.csv', sep = ',', header = TRUE, stringsAsFactors = FALSE, encoding = 'UTF-8') energy_waste <- panel_fv %>% group_by(RE_ID) %>% summarise(Oil_MMB_Day = mean(Oil_B_Day, na.rm = TRUE) / 1000000, TotGas_MMB_Day = mean(TotGas_B_Day, na.rm = TRUE) / 1000000, Flared_Gas_MMB_Day = mean(Flared_Gas_B_Day, na.rm = TRUE) / 1000000, Non_Routine_Flaring_MMB_Day = mean(Non_Routine_Flaring_B_Day, na.rm = TRUE) / 1000000, Routine_Flaring_MMB_Day = mean(Routine_Flaring_B_Day , na.rm = TRUE) / 1000000, Vented_Gas_MMB_Day = mean(Vented_Gas_B_Day, na.rm = TRUE) / 1000000, IVented_Gas_MMB_Day = mean(IVented_Gas_B_Day, na.rm = TRUE) / 1000000, UVented_Gas_MMB_Day = mean(UVented_Gas_B_Day, na.rm = TRUE) / 1000000) %>% ungroup() %>% mutate(Energy_MMB_Day = Oil_MMB_Day + TotGas_MMB_Day, Energy_Waste_MMB_Day_without_reform = Non_Routine_Flaring_MMB_Day + Routine_Flaring_MMB_Day + IVented_Gas_MMB_Day + UVented_Gas_MMB_Day, Energy_Waste_MMB_Day_with_reform = Non_Routine_Flaring_MMB_Day + UVented_Gas_MMB_Day, Energy_Savings_MMB_Day_with_reform = Routine_Flaring_MMB_Day + IVented_Gas_MMB_Day) # 3.1 Wasted Energy sum(energy_waste$Energy_Waste_MMB_Day_without_reform) / sum(energy_waste$Energy_MMB_Day) * 100 # 2.78% of all energy extracted wasted sum(energy_waste$Energy_Waste_MMB_Day_without_reform) / sum(energy_waste$Oil_MMB_Day) * 100 # 4.44% of the oil energy extracted wasted sum(energy_waste$Energy_Waste_MMB_Day_without_reform) / sum(energy_waste$TotGas_MMB_Day) * 100 # 7.09% of the natural gas extracted wasted sum(energy_waste$Energy_Waste_MMB_Day_without_reform) # 0.1966 MM BOE Day aggregate energy waste mean(energy_waste$Energy_Waste_MMB_Day_without_reform) # 0.00035 MM BOE Day average energy waste median(energy_waste$Energy_Waste_MMB_Day_without_reform) # 0.00014 MM BOE Day median energy waste sum(energy_waste$Energy_Waste_MMB_Day_with_reform) / sum(energy_waste$Energy_MMB_Day) * 100 # 1.09% of all energy extracted wasted sum(energy_waste$Energy_Waste_MMB_Day_with_reform) / sum(energy_waste$Oil_MMB_Day) * 100 # 1.74% of the oil energy extracted wasted sum(energy_waste$Energy_Waste_MMB_Day_with_reform) / sum(energy_waste$TotGas_MMB_Day) * 100 # 2.92% of the natural gas extracted wasted sum(energy_waste$Energy_Waste_MMB_Day_with_reform) # 0.07715 MM BOE Day aggregate energy savings mean(energy_waste$Energy_Waste_MMB_Day_with_reform) # 0.0001388 MM BOE Day average energy waste median(energy_waste$Energy_Waste_MMB_Day_with_reform) # 0.00001819 MM BOE Day median energy waste sum(energy_waste$Energy_Waste_MMB_Day_with_reform) / sum(energy_waste$Energy_Waste_MMB_Day_without_reform) * 100 # 39.25% savings from the reform # 3.2 Recovered Energy sum(energy_waste$Routine_Flaring_MMB_Day) # 0.08602 MM BOE Day sum(energy_waste$IVented_Gas_MMB_Day) # 0.03338 MM BOE Day # 3.3 Unrecoverable Energy sum(energy_waste$Non_Routine_Flaring_MMB_Day) # 0.001196 MM BOE Day sum(energy_waste$UVented_Gas_MMB_Day) # 0.075950 MM BOE Day # 3.4 Wasted Natural Gas: Flaring & Venting sum(energy_waste$Energy_Waste_MMB_Day_without_reform) / sum(energy_waste$TotGas_MMB_Day) * 100 # average waste 7.088% of all energy natural gas sum(energy_waste$Non_Routine_Flaring_MMB_Day) / sum(energy_waste$Energy_MMB_Day) * 100 # average non.routine flaring 0.04% of the total natural gas extracted sum(energy_waste$Flared_Gas_MMB_Day) / sum(energy_waste$Energy_MMB_Day) * 100 # average flaring 3.30% of the total natural gas extracted sum(energy_waste$Vented_Gas_MMB_Day) / sum(energy_waste$Energy_MMB_Day) * 100 # average venting 3.788% of the total natural gas extracted # 3.5 Wasted Natural Gas: Intentional & Unintentional Venting sum(energy_waste$IVented_Gas_MMB_Day) / sum(energy_waste$Energy_MMB_Day) * 100 # average intentional venting 1.263% of total natural gas extracted sum(energy_waste$UVented_Gas_MMB_Day) / sum(energy_waste$Energy_MMB_Day) * 100 # average unintentional venting 2.873% of total natural gas extracted # 3.6 Wasted Natural Gas: Intentional-Unintentional Ratio sum(energy_waste$IVented_Gas_MMB_Day) / sum(energy_waste$Vented_Gas_MMB_Day) * 100 # average intentional venting 16.12% of total venting sum(energy_waste$UVented_Gas_MMB_Day) / sum(energy_waste$Vented_Gas_MMB_Day) * 100 # average unintentional venting 83.88% of total venting rm(list = ls(all=TRUE)) ##-----------------------------------------------------------------------------------------------# ## 4. CO2eq Savings 2005-2020 ##-----------------------------------------------------------------------------------------------# panel_fv <- 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/Second_Stage_Results.csv', sep = ',', header = TRUE, stringsAsFactors = FALSE, encoding = 'UTF-8') panel_fv <- panel_fv %>% mutate(Non_Routine_Flaring_TCO2e_Day = (Non_Routine_Flaring_B_Day * 0.30180539 * 0.98) + (Non_Routine_Flaring_B_Day * 3.95833079 * 0.02), Routine_Flaring_TCO2e_Day = (Routine_Flaring_B_Day * 0.30180539 * 0.98) + (Routine_Flaring_B_Day * 3.95833079 * 0.02), Flared_TCO2e_Day = Non_Routine_Flaring_TCO2e_Day + Routine_Flaring_TCO2e_Day, IVented_TCO2e_Day = IVented_Gas_B_Day * 3.95833079, UVented_TCO2e_Day = UVented_Gas_B_Day * 3.95833079, Vented_TCO2e_Day = IVented_TCO2e_Day + UVented_TCO2e_Day, Methane_TCO2e_Day = (Non_Routine_Flaring_B_Day * 3.95833079 * 0.02) + (Routine_Flaring_B_Day * 3.95833079 * 0.02) + Vented_TCO2e_Day) co2e_emissions <- panel_fv %>% group_by(RE_ID) %>% summarise(Flared_MMTCO2e_Day = mean(Flared_TCO2e_Day, na.rm = TRUE) / 1000000, Non_Routine_Flaring_MMTCO2e_Day = mean(Non_Routine_Flaring_TCO2e_Day, na.rm = TRUE) / 1000000, Routine_Flaring_MMTCO2e_Day = mean(Routine_Flaring_TCO2e_Day, na.rm = TRUE) / 1000000, Vented_MMTCO2e_Day = mean(Vented_TCO2e_Day, na.rm = TRUE) / 1000000, IVented_MMTCO2e_Day = mean(IVented_TCO2e_Day, na.rm = TRUE) / 1000000, UVented_MMTCO2e_Day = mean(UVented_TCO2e_Day, na.rm = TRUE) / 1000000, Methane_MMTCO2e_Day = mean(Methane_TCO2e_Day, na.rm = TRUE) / 1000000) %>% ungroup() %>% mutate(Emissions_MMTCO2e_Day_without_reform = Non_Routine_Flaring_MMTCO2e_Day + Routine_Flaring_MMTCO2e_Day + IVented_MMTCO2e_Day + UVented_MMTCO2e_Day, Emissions_MMTCO2e_Day_with_reform = Non_Routine_Flaring_MMTCO2e_Day + UVented_MMTCO2e_Day, Emissions_Savings_MMTCO2e_Day_with_reform = Routine_Flaring_MMTCO2e_Day + IVented_MMTCO2e_Day) # 4.1 Emissions sum(co2e_emissions$Emissions_MMTCO2e_Day_without_reform) # 0.4655 MMTCO2e/Day sum(co2e_emissions$Emissions_MMTCO2e_Day_with_reform) # 0.3011 MMTCO2e/Day sum(co2e_emissions$Emissions_Savings_MMTCO2e_Day_with_reform) # 0.1644 MMTCO2e/Day # 4.2 Saved Emissions sum(co2e_emissions$Routine_Flaring_MMTCO2e_Day) # 0.0322 MMTCO2e/Day sum(co2e_emissions$IVented_MMTCO2e_Day) # 0.1321 MMTCO2e/Day # 4.3 Unsaved Emissions sum(co2e_emissions$Non_Routine_Flaring_MMTCO2e_Day) # 0.0004 MMTCO2e/Day sum(co2e_emissions$UVented_MMTCO2e_Day) # 0.3006 MMTCO2e/Day # 4.4. Savable Emissions (sum(co2e_emissions$Routine_Flaring_MMTCO2e_Day) + sum(co2e_emissions$IVented_MMTCO2e_Day)) / sum(co2e_emissions$Emissions_MMTCO2e_Day_without_reform) * 100 # 35.31% # 4.5 Savable CO2 savings sum(co2e_emissions$Routine_Flaring_MMTCO2e_Day) / (sum(co2e_emissions$Routine_Flaring_MMTCO2e_Day) + sum(co2e_emissions$Non_Routine_Flaring_MMTCO2e_Day)) * 100 # 98.62% # 4.5 Savable CH4 savings sum(co2e_emissions$IVented_MMTCO2e_Day) / sum(co2e_emissions$Methane_MMTCO2e_Day) * 100 # 30.05% # 4.6 Savings with 91% efficiency of the flaring stuck panel_fv <- panel_fv %>% mutate(Non_Routine_Flaring_TCO2e_Day = (Non_Routine_Flaring_B_Day * 0.30180539 * 0.91) + (Non_Routine_Flaring_B_Day * 3.95833079 * 0.09), Routine_Flaring_TCO2e_Day = (Routine_Flaring_B_Day * 0.30180539 * 0.91) + (Routine_Flaring_B_Day * 3.95833079 * 0.09), Flared_TCO2e_Day = Non_Routine_Flaring_TCO2e_Day + Routine_Flaring_TCO2e_Day, IVented_TCO2e_Day = IVented_Gas_B_Day * 3.95833079, UVented_TCO2e_Day = UVented_Gas_B_Day * 3.95833079, Vented_TCO2e_Day = IVented_TCO2e_Day + UVented_TCO2e_Day, Methane_TCO2e_Day = (Non_Routine_Flaring_B_Day * 3.95833079 * 0.02) + (Routine_Flaring_B_Day * 3.95833079 * 0.02) + Vented_TCO2e_Day) co2e_emissions <- panel_fv %>% group_by(RE_ID) %>% summarise(Flared_MMTCO2e_Day = mean(Flared_TCO2e_Day, na.rm = TRUE) / 1000000, Non_Routine_Flaring_MMTCO2e_Day = mean(Non_Routine_Flaring_TCO2e_Day, na.rm = TRUE) / 1000000, Routine_Flaring_MMTCO2e_Day = mean(Routine_Flaring_TCO2e_Day, na.rm = TRUE) / 1000000, Vented_MMTCO2e_Day = mean(Vented_TCO2e_Day, na.rm = TRUE) / 1000000, IVented_MMTCO2e_Day = mean(IVented_TCO2e_Day, na.rm = TRUE) / 1000000, UVented_MMTCO2e_Day = mean(UVented_TCO2e_Day, na.rm = TRUE) / 1000000, Methane_MMTCO2e_Day = mean(Methane_TCO2e_Day, na.rm = TRUE) / 1000000) %>% ungroup() %>% mutate(Emissions_MMTCO2e_Day_without_reform = Non_Routine_Flaring_MMTCO2e_Day + Routine_Flaring_MMTCO2e_Day + IVented_MMTCO2e_Day + UVented_MMTCO2e_Day, Emissions_MMTCO2e_Day_with_reform = Non_Routine_Flaring_MMTCO2e_Day + UVented_MMTCO2e_Day, Emissions_Savings_MMTCO2e_Day_with_reform = Routine_Flaring_MMTCO2e_Day + IVented_MMTCO2e_Day) # 4.1 Emissions sum(co2e_emissions$Emissions_MMTCO2e_Day_without_reform) # 0.4515 MMTCO2e/Day sum(co2e_emissions$Emissions_MMTCO2e_Day_with_reform) # 0.3333 MMTCO2e/Day sum(co2e_emissions$Emissions_Savings_MMTCO2e_Day_with_reform) # 0.1182 MMTCO2e/Day # 4.2 Saved Emissions sum(co2e_emissions$Routine_Flaring_MMTCO2e_Day) # 0.05427 MMTCO2e/Day sum(co2e_emissions$IVented_MMTCO2e_Day) # 0.06391 MMTCO2e/Day # 4.3 Unsaved Emissions sum(co2e_emissions$Non_Routine_Flaring_MMTCO2e_Day) # 0.0007547 MMTCO2e/Day sum(co2e_emissions$UVented_MMTCO2e_Day) # 0.3325 MMTCO2e/Day # 4.4. Savable Emissions (sum(co2e_emissions$Routine_Flaring_MMTCO2e_Day) + sum(co2e_emissions$IVented_MMTCO2e_Day)) / sum(co2e_emissions$Emissions_MMTCO2e_Day_without_reform) * 100 # 26.18% # 4.5 Savable CO2 savings sum(co2e_emissions$Routine_Flaring_MMTCO2e_Day) / (sum(co2e_emissions$Routine_Flaring_MMTCO2e_Day) + sum(co2e_emissions$Non_Routine_Flaring_MMTCO2e_Day)) * 100 # 98.62% # 4.5 Savable CH4 savings sum(co2e_emissions$IVented_MMTCO2e_Day) / sum(co2e_emissions$Methane_MMTCO2e_Day) * 100 # 15.84% # 0.1182 - 0.0961 # line 302 - line 255 rm(list = ls(all=TRUE)) ##-----------------------------------------------------------------------------------------------# ## 5. Taxation Structure ##-----------------------------------------------------------------------------------------------# panel_fv <- 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/Second_Stage_Results_NEW.csv', sep = ',', header = TRUE, stringsAsFactors = FALSE, encoding = 'UTF-8') summary(panel_fv) estimate_mode <- function(x) { d <- density(x) d$x[which.max(d$y)] } #estimate_mode(panel_fv$GOR_B_B) # 5.1 Average across all the observed years agg_ivent <- panel_fv %>% group_by(RE_ID) %>% summarise(IVented_Gas_B_Day = mean(IVented_Gas_B_Day)) %>% ungroup() # 5.2 Total IVent Emissions sum(agg_ivent$IVented_Gas_B_Day) # 33378 BOE/Day # 5.3 IVent under different saving scenarios IVented_Gas_B_Day_100 <- 0 # 0.00 BOE/Day IVented_Gas_B_Day_90 <- 0.1 * sum(agg_ivent$IVented_Gas_B_Day) # 3337.82 BOE/Day IVented_Gas_B_Day_80 <- 0.2 * sum(agg_ivent$IVented_Gas_B_Day) # 6675.65 BOE/Day IVented_Gas_B_Day_70 <- 0.3 * sum(agg_ivent$IVented_Gas_B_Day) # 10013.47 BOE/Day IVented_Gas_B_Day_60 <- 0.4 * sum(agg_ivent$IVented_Gas_B_Day) # 13351.30 BOE/Day IVented_Gas_B_Day_50 <- 0.5 * sum(agg_ivent$IVented_Gas_B_Day) # 16689.12 BOE/Day IVented_Gas_B_Day_40 <- 0.6 * sum(agg_ivent$IVented_Gas_B_Day) # 20026.95 BOE/Day IVented_Gas_B_Day_30 <- 0.7 * sum(agg_ivent$IVented_Gas_B_Day) # 23364.77 BOE/Day IVented_Gas_B_Day_20 <- 0.8 * sum(agg_ivent$IVented_Gas_B_Day) # 26702.60 BOE/Day IVented_Gas_B_Day_10 <- 0.9 * sum(agg_ivent$IVented_Gas_B_Day) # 30040.42 BOE/Day IVented_Gas_B_Day_0 <- 1 * sum(agg_ivent$IVented_Gas_B_Day) # 33378.25 BOE/Day # Energy Savings sum(agg_ivent$IVented_Gas_B_Day) - (0 * sum(agg_ivent$IVented_Gas_B_Day)) sum(agg_ivent$IVented_Gas_B_Day) - (0.1 * sum(agg_ivent$IVented_Gas_B_Day)) sum(agg_ivent$IVented_Gas_B_Day) - (0.2 * sum(agg_ivent$IVented_Gas_B_Day)) sum(agg_ivent$IVented_Gas_B_Day) - (0.3 * sum(agg_ivent$IVented_Gas_B_Day)) sum(agg_ivent$IVented_Gas_B_Day) - (0.4 * sum(agg_ivent$IVented_Gas_B_Day)) sum(agg_ivent$IVented_Gas_B_Day) - (0.5 * sum(agg_ivent$IVented_Gas_B_Day)) sum(agg_ivent$IVented_Gas_B_Day) - (0.6 * sum(agg_ivent$IVented_Gas_B_Day)) sum(agg_ivent$IVented_Gas_B_Day) - (0.7 * sum(agg_ivent$IVented_Gas_B_Day)) sum(agg_ivent$IVented_Gas_B_Day) - (0.8 * sum(agg_ivent$IVented_Gas_B_Day)) sum(agg_ivent$IVented_Gas_B_Day) - (0.9 * sum(agg_ivent$IVented_Gas_B_Day)) sum(agg_ivent$IVented_Gas_B_Day) - (1 * sum(agg_ivent$IVented_Gas_B_Day)) # List of possible Natural Gas Rebates Max_Sub_Dollar_B <- unique(max(agg_ivent$IVented_Gas_B_Day) / (panel_fv$delta2_hat + panel_fv$delta5_hat)) Sub_Dollar_B <- seq(Max_Sub_Dollar_B, 0, 0.1) Sub_Dollar_B_75 <- panel_fv %>% group_by(RE_ID) %>% summarise(IVented_Gas_B_Day = mean(IVented_Gas_B_Day), delta2_hat = mean(delta2_hat), delta5_hat = mean(delta5_hat)) %>% ungroup() %>% crossing(., Sub_Dollar_B) %>% mutate(NewIVented_Gas_B_Day = ifelse(IVented_Gas_B_Day - ((delta2_hat + delta5_hat) * Sub_Dollar_B) > 0, IVented_Gas_B_Day - ((delta2_hat + delta5_hat) * Sub_Dollar_B), 0)) %>% group_by(Sub_Dollar_B) %>% summarise(NewIVented_Gas_B_Day = sum(NewIVented_Gas_B_Day, na.rm = TRUE)) %>% ungroup() %>% slice(which.min(abs(NewIVented_Gas_B_Day - IVented_Gas_B_Day_90))) summary(Sub_Dollar_B_75) # 100% -> 78.63 $/BOE # 90% -> 20.73 $/BOE # 80% -> 17.53 $/BOE # 70% -> 14.93 $/BOE # 60% -> 12.60 $/BOE # 50% -> 10.43 $/BOE # 40% -> 8.33 $/BOE # 30% -> 6.24 $/BOE # 20% -> 4.13 $/BOE # 10% -> 2.03 $/BOE # 0% -> 0.00 $/BOE # Oil Tax, Energy Savings, Routine Flaring, Intentional Venting, Agg CO2e Savings panel_econ_75 <- panel_fv %>% mutate(GOR_B_B = TotGas_B_Day/Oil_B_Day, Av_Gas_Price_Dollar_B = mean(HH_Dollars_B), Av_Oil_Price_Dollar_B = mean(WTI_Dollars_B)) %>% group_by(RE_ID) %>% summarise(Geo = unique(Geo), GOR_B_B = mean(GOR_B_B), Oil_B_Day = mean(Oil_B_Day), TotGas_B_Day = mean(TotGas_B_Day), Sold_Gas_B_Day = mean(Sold_Gas_B_Day), Non_Routine_Flaring_B_Day = mean(Non_Routine_Flaring_B_Day), Routine_Flaring_B_Day = mean(Routine_Flaring_B_Day), IVented_Gas_B_Day = mean(IVented_Gas_B_Day), UVented_Gas_B_Day = mean(UVented_Gas_B_Day), Insitu_Gas_B_Day = mean(Insitu_Gas_B_Day), Re_Gas_B_Day = mean(Re_Gas_B_Day), alpha0_hat = mean(alpha0_hat), alpha1_hat = mean(alpha1_hat), delta2_hat = mean(delta2_hat), delta4_hat = mean(delta4_hat), delta5_hat = mean(delta5_hat), Av_Gas_Price_Dollar_B = mean(Av_Gas_Price_Dollar_B), Av_Oil_Price_Dollar_B = mean(Av_Oil_Price_Dollar_B), Longitude = mean(Longitude), Latitude = mean(Latitude)) %>% crossing(., Sub_Dollar_B) %>% mutate(NewIVented_Gas_B_Day = ifelse(IVented_Gas_B_Day - ((delta2_hat + delta5_hat) * Sub_Dollar_B) > 0, IVented_Gas_B_Day - ((delta2_hat + delta5_hat) * Sub_Dollar_B), 0), TrueNewIVent_Gas_B_Day = ifelse(IVented_Gas_B_Day - (delta2_hat * Sub_Dollar_B) > 0, IVented_Gas_B_Day - (delta2_hat * Sub_Dollar_B), 0)) %>% filter(Sub_Dollar_B == Sub_Dollar_B_75$Sub_Dollar_B) %>% mutate(Gas_Tax_Dollar_B = Sub_Dollar_B_75$Sub_Dollar_B, Oil_Tax_Dollar_B = -(Gas_Tax_Dollar_B * GOR_B_B * (1 - (delta4_hat / GOR_B_B))), NewGas_Price_Dollar_B = Av_Gas_Price_Dollar_B + ((-1) * Gas_Tax_Dollar_B), NewRoutine_Flaring_B_Day = ifelse(alpha0_hat + (alpha1_hat * NewGas_Price_Dollar_B) > 0, alpha0_hat + (alpha1_hat * NewGas_Price_Dollar_B), 0), Savings_Flaring_B_Day = ifelse(Routine_Flaring_B_Day - NewRoutine_Flaring_B_Day > 0, Routine_Flaring_B_Day - NewRoutine_Flaring_B_Day , 0), Savings_IVented_B_Day = IVented_Gas_B_Day - NewIVented_Gas_B_Day, Savings_Tot_B_Day = Savings_Flaring_B_Day + Savings_IVented_B_Day) %>% select(RE_ID, Geo, Gas_Tax_Dollar_B, Oil_Tax_Dollar_B, Savings_Tot_B_Day, Sold_Gas_B_Day, Insitu_Gas_B_Day, Re_Gas_B_Day, Non_Routine_Flaring_B_Day, Routine_Flaring_B_Day, NewRoutine_Flaring_B_Day, Savings_Flaring_B_Day, IVented_Gas_B_Day, NewIVented_Gas_B_Day, Savings_IVented_B_Day, UVented_Gas_B_Day, GOR_B_B, Oil_B_Day, TotGas_B_Day, Av_Gas_Price_Dollar_B, Av_Oil_Price_Dollar_B, TrueNewIVent_Gas_B_Day) %>% drop_na() # Summary Taxation summary(panel_econ_75$Oil_Tax_Dollar_B) estimate_mode(panel_econ_75$Oil_Tax_Dollar_B) sd(panel_econ_75$Oil_Tax_Dollar_B, na.rm = T) # Summary Flaring & Venting Savings summary(panel_econ_75$Savings_Flaring_B_Day) sd(panel_econ_75$Savings_Flaring_B_Day) summary(panel_econ_75$Savings_IVented_B_Day) sd(panel_econ_75$Savings_IVented_B_Day) # Total Savings sum(panel_econ_75$Savings_Flaring_B_Day) sum(panel_econ_75$Savings_IVented_B_Day) sum(panel_econ_75$Savings_Flaring_B_Day) + sum(panel_econ_75$Savings_IVented_B_Day) (sum(panel_econ_75$Savings_Flaring_B_Day) + sum(panel_econ_75$Savings_IVented_B_Day)) / (sum(panel_econ_75$Oil_B_Day) + sum(panel_econ_75$TotGas_B_Day)) * 100 # Number of Oilfields #num_oilfields_zeros <- sum(panel_econ_75$NewIVented_Gas_B_Day == 0) #total <- nrow(panel_econ_75) #(num_oilfields_zeros / total) * 100 panel_econ_75 <- panel_econ_75 %>% mutate(Compensation_gas = (TotGas_B_Day - NewRoutine_Flaring_B_Day - TrueNewIVent_Gas_B_Day) * Gas_Tax_Dollar_B, Compensation_oil = Oil_B_Day * Oil_Tax_Dollar_B) avg_com <- panel_econ_75 %>% summarise(comp_gas = sum(Compensation_gas), Comp_oil = sum(Compensation_oil)) -1 * avg_com$comp_gas / avg_com$Comp_oil * 100 # Environmental Savings 98% Flaring Efficiency panel_env_75_98 <- panel_econ_75 %>% mutate(Non_Routine_Flaring_TCO2e_Day = (Non_Routine_Flaring_B_Day * 0.30180539 * 0.98) + (Non_Routine_Flaring_B_Day * 3.95833079 * 0.02), Routine_Flaring_TCO2e_Day = (Routine_Flaring_B_Day * 0.30180539 * 0.98) + (Routine_Flaring_B_Day * 3.95833079 * 0.02), NewRoutine_Flaring_TCO2e_Day = (NewRoutine_Flaring_B_Day * 0.30180539 * 0.98) + (NewRoutine_Flaring_B_Day * 3.95833079 * 0.02), IVented_Gas_TCO2e_Day = IVented_Gas_B_Day * 3.95833079, NewIVented_Gas_TCO2e_Day = NewIVented_Gas_B_Day * 3.95833079, UVented_Gas_TCO2e_Day = UVented_Gas_B_Day * 3.95833079) %>% select(RE_ID, Gas_Tax_Dollar_B, Oil_Tax_Dollar_B, Non_Routine_Flaring_TCO2e_Day, Routine_Flaring_TCO2e_Day, NewRoutine_Flaring_TCO2e_Day, IVented_Gas_TCO2e_Day, NewIVented_Gas_TCO2e_Day, UVented_Gas_TCO2e_Day) %>% mutate(Savings_Flaring_TCO2e_Day = ifelse(Routine_Flaring_TCO2e_Day - NewRoutine_Flaring_TCO2e_Day > 0, Routine_Flaring_TCO2e_Day - NewRoutine_Flaring_TCO2e_Day, 0), Savings_IVented_TCO2e_Day = IVented_Gas_TCO2e_Day - NewIVented_Gas_TCO2e_Day, Savings_TCO2e_Day = Savings_Flaring_TCO2e_Day + Savings_IVented_TCO2e_Day, Emissions_TCO2e_Day = Non_Routine_Flaring_TCO2e_Day + Routine_Flaring_TCO2e_Day + IVented_Gas_TCO2e_Day + UVented_Gas_TCO2e_Day) # Summary Flaring & Venting Savings summary(panel_env_75_98$Savings_Flaring_TCO2e_Day) estimate_mode(panel_env_75_98$Savings_Flaring_TCO2e_Day) sd(panel_env_75_98$Savings_Flaring_TCO2e_Day) summary(panel_env_75_98$Savings_IVented_TCO2e_Day) estimate_mode(panel_env_75_98$Savings_IVented_TCO2e_Day) sd(panel_env_75_98$Savings_IVented_TCO2e_Day) # Total Savings sum(panel_env_75_98$Savings_Flaring_TCO2e_Day) + sum(panel_env_75_98$Savings_IVented_TCO2e_Day) # Percentage Savings (sum(panel_env_75_98$Savings_Flaring_TCO2e_Day) + sum(panel_env_75_98$Savings_IVented_TCO2e_Day)) / (sum(panel_env_75_98$Non_Routine_Flaring_TCO2e_Day) + sum(panel_env_75_98$Routine_Flaring_TCO2e_Day) + sum(panel_env_75_98$IVented_Gas_TCO2e_Day) + sum(panel_env_75_98$UVented_Gas_TCO2e_Day)) * 100 # 23.41% # Environmental Savings 91% Flaring Efficiency panel_env_75_91 <- panel_econ_75 %>% mutate(Non_Routine_Flaring_TCO2e_Day = (Non_Routine_Flaring_B_Day * 0.30180539 * 0.91) + (Non_Routine_Flaring_B_Day * 3.95833079 * 0.09), Routine_Flaring_TCO2e_Day = (Routine_Flaring_B_Day * 0.30180539 * 0.91) + (Routine_Flaring_B_Day * 3.95833079 * 0.09), NewRoutine_Flaring_TCO2e_Day = (NewRoutine_Flaring_B_Day * 0.30180539 * 0.91) + (NewRoutine_Flaring_B_Day * 3.95833079 * 0.09), IVented_Gas_TCO2e_Day = IVented_Gas_B_Day * 3.95833079, NewIVented_Gas_TCO2e_Day = NewIVented_Gas_B_Day * 3.95833079, UVented_Gas_TCO2e_Day = UVented_Gas_B_Day * 3.95833079) %>% select(RE_ID, Gas_Tax_Dollar_B, Oil_Tax_Dollar_B, Non_Routine_Flaring_TCO2e_Day, Routine_Flaring_TCO2e_Day, NewRoutine_Flaring_TCO2e_Day, IVented_Gas_TCO2e_Day, NewIVented_Gas_TCO2e_Day, UVented_Gas_TCO2e_Day) %>% mutate(Savings_Flaring_TCO2e_Day = ifelse(Routine_Flaring_TCO2e_Day - NewRoutine_Flaring_TCO2e_Day > 0, Routine_Flaring_TCO2e_Day - NewRoutine_Flaring_TCO2e_Day, 0), Savings_IVented_TCO2e_Day = IVented_Gas_TCO2e_Day - NewIVented_Gas_TCO2e_Day, Savings_TCO2e_Day = Savings_Flaring_TCO2e_Day + Savings_IVented_TCO2e_Day, Emissions_TCO2e_Day = Non_Routine_Flaring_TCO2e_Day + Routine_Flaring_TCO2e_Day + IVented_Gas_TCO2e_Day + UVented_Gas_TCO2e_Day) (sum(panel_env_75_91$Savings_Flaring_TCO2e_Day) - sum(panel_env_75_98$Savings_Flaring_TCO2e_Day)) / 1000000 # 0.01142 MMTCO2e/Day # Savings from 91% efficiency for the 100% scenario: 0.02816 MMTCO2e/Day # Savings from 98% efficiency for the 100% scenario: 0.01673 MMTCO2e/Day rm(list = ls(all=TRUE)) ##-----------------------------------------------------------------------------------------------# ## 6. Figure 1 ##-----------------------------------------------------------------------------------------------# panel_fv <- 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/Second_Stage_Results.csv', sep = ',', header = TRUE, stringsAsFactors = FALSE, encoding = 'UTF-8') panel_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/F&V_Dataset.csv', sep = ',', header = TRUE, stringsAsFactors = FALSE, encoding = 'UTF-8') cross_gas <- panel_gas %>% group_by(RE_ID) %>% summarise(Sold_Gas_B_Day = mean(Sold_Gas_B_Day, na.rm = TRUE), Insitu_Gas_B_Day = mean(Insitu_Gas_B_Day, na.rm = TRUE), Re_Gas_B_Day = mean(Re_Gas_B_Day, na.rm = TRUE)) %>% ungroup() # 6.1 Cross-Section cross_fv <- panel_fv %>% drop_na(IVented_Gas_B_Day) %>% group_by(RE_ID) %>% summarise(Longitude = mean(Longitude), Latitude = mean(Latitude), Flared_Gas_B_Day = mean(Flared_Gas_B_Day, na.rm = TRUE), Vented_Gas_B_Day = mean(Vented_Gas_B_Day, na.rm = TRUE), FVR_B_B = atan(Vented_Gas_B_Day / Flared_Gas_B_Day), Wasted_Energy_B_Day = Flared_Gas_B_Day + Vented_Gas_B_Day, Oil_B_Day = mean(Oil_B_Day, na.rm = TRUE)) %>% ungroup() cross_fv <- left_join(cross_fv, cross_gas) cross_fv_plot <- cross_fv %>% drop_na(Longitude, Latitude) %>% filter(Latitude < 50) # 6.2 US Map map_fv <- ggplot() + geom_polygon(data = map_data('state'), aes(x = long, y = lat, group = group), colour = 'gray', fill = 'white', alpha = 0.3, size = 0.7/.pt) + geom_point(data = cross_fv_plot, aes(x = Longitude, y = Latitude, color = FVR_B_B, size = Wasted_Energy_B_Day), alpha = 0.3) + scale_color_continuous(low = '#FF6C00', high = '#0070FF', na.value = 'white', name = '', limits = c(0, 1.6), breaks = c(0, 1.6), labels = c('%Flaring','%Venting')) + # limits of FVR_B_B #scale_size_continuous(range = c(0.4, 4), breaks = seq(0, 10.70, by = 2.5), name = 'BOE/Day') # min(Wasted_Energy_MMB_Day) - max(Wasted_Energy_MMB_Day) scale_size_continuous(range = c(0.4, 4), breaks = seq(0, 4000, by = 1000), name = 'Natural Gas Waste (BOE/Day)') + # min(Wasted_Energy_MMB_Day) - max(Wasted_Energy_MMB_Day) theme_void() + theme(legend.direction = 'horizontal', legend.position = 'bottom', legend.title = element_text(size = 12), legend.text = element_text(size = 12), text = element_text(size = 12, family = 'helvetica'), legend.spacing.x = unit(1.2, 'cm')) + ggtitle('(A)') + 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')) + coord_map() map_fv # 6.3 Environmental Merit Base Curve cross_fv <- cross_fv[order(cross_fv$Wasted_Energy_B_Day),] cross_fv <- cross_fv %>% mutate(Flared_Gas_TCO2e_Day = (Flared_Gas_B_Day * 0.30180539 * 0.98) + (Flared_Gas_B_Day * 3.95833079 * 0.02), Vented_Gas_TCO2e_Day = Vented_Gas_B_Day * 3.95833079, Emissions_TCO2e_Day = Flared_Gas_TCO2e_Day + Vented_Gas_TCO2e_Day, Cum_Prod_B_Day = cumsum(Oil_B_Day), Lag_Cum_Prod_B_Day = Cum_Prod_B_Day - Oil_B_Day, Cum_Emissions_TCO2e_Day = (cumsum(Oil_B_Day * Emissions_TCO2e_Day)) / Cum_Prod_B_Day, Cum_Prod_per = Cum_Prod_B_Day / max(Cum_Prod_B_Day) * 100, Lag_Cum_Prod_per = Lag_Cum_Prod_B_Day / max(Lag_Cum_Prod_B_Day) * 100) figure_1 <- map_fv + ggplot(data = cross_fv, aes(ymin = 0)) + geom_rect(aes(xmin = Lag_Cum_Prod_per, xmax = Cum_Prod_per, ymax = Wasted_Energy_B_Day, fill = FVR_B_B), size = 1) + # fill = Geo FVR_B_B scale_fill_continuous(low = '#FF6C00', high = '#0070FF', na.value = 'white', name = '', limits = c(0, 1.6), breaks = c(0, 1.6), labels = c('Flaring','Venting')) + # limits of FVR_B_B theme_minimal() + labs(x = 'Cumulative Production (%)', y = 'Wasted Energy (BOE/Day)') + theme(legend.direction = 'horizontal', legend.position = 'none', text = element_text(size = 12, family = 'helvetica'), legend.spacing.x = unit(1.2, 'cm'), panel.background = element_rect(colour = 'white', size = 0.5), axis.text = element_text(size = 12), axis.ticks = element_line(color = 'black')) + ggtitle('(B)') + guides(fill = guide_legend(title.position = 'top', title.hjust = 0.5, label.position = 'top')) figure_1 ggsave("C:/Users/s15535/Desktop/2 - Energy Economics/2 - Oil & Natural Gas/2.5 - Flaring & Venting/1 - American FV&L/0 - Manuscript/Figure_1.svg", height = 15, width = 30, units = c('cm'), dpi = 400) rm(list = ls(all=TRUE)) ##-----------------------------------------------------------------------------------------------# ## 7. Figure 2 ##-----------------------------------------------------------------------------------------------# panel_fv <- 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/Second_Stage_Results.csv', sep = ',', header = TRUE, stringsAsFactors = FALSE, encoding = 'UTF-8') #devtools::install_github("ricardo-bion/ggradar", dependencies = TRUE) Sub_Dollar_B_seq <- seq(-80, 0, 1) panel_econ <- panel_fv %>% mutate(GOR_B_B = TotGas_B_Day/Oil_B_Day, Av_Gas_Price_Dollar_B = mean(HH_Dollars_B), Av_Oil_Price_Dollar_B = mean(WTI_Dollars_B)) %>% group_by(RE_ID) %>% summarise(GOR_B_B = mean(GOR_B_B), Oil_B_Day = mean(Oil_B_Day), TotGas_B_Day = mean(TotGas_B_Day), Non_Routine_Flaring_B_Day = mean(Non_Routine_Flaring_B_Day), Routine_Flaring_B_Day = mean(Routine_Flaring_B_Day), IVented_Gas_B_Day = mean(IVented_Gas_B_Day), UVented_Gas_B_Day = mean(UVented_Gas_B_Day), alpha0_hat = mean(alpha0_hat), alpha1_hat = mean(alpha1_hat), delta2_hat = mean(delta2_hat), delta4_hat = mean(delta4_hat), delta5_hat = mean(delta5_hat), Av_Gas_Price_Dollar_B = mean(Av_Gas_Price_Dollar_B), Av_Oil_Price_Dollar_B = mean(Av_Oil_Price_Dollar_B), Longitude = mean(Longitude), Latitude = mean(Latitude)) %>% crossing(., Sub_Dollar_B_seq) %>% mutate(NewIVented_Gas_B_Day = ifelse(IVented_Gas_B_Day - ((delta2_hat + delta5_hat) * Sub_Dollar_B_seq) > 0, IVented_Gas_B_Day - ((delta2_hat + delta5_hat) * Sub_Dollar_B_seq), 0)) %>% rename(Gas_Tax_Dollar_B = Sub_Dollar_B_seq) %>% mutate(Oil_Tax_Dollar_B = -(Gas_Tax_Dollar_B * GOR_B_B * (1 - (delta4_hat / GOR_B_B))), NewGas_Price_Dollar_B = Av_Gas_Price_Dollar_B + ((-1) * Gas_Tax_Dollar_B), NewRoutine_Flaring_B_Day = ifelse(alpha0_hat + (alpha1_hat * NewGas_Price_Dollar_B) > 0, alpha0_hat + (alpha1_hat * NewGas_Price_Dollar_B), 0), Savings_Flaring_B_Day = ifelse(Routine_Flaring_B_Day - NewRoutine_Flaring_B_Day > 0, Routine_Flaring_B_Day - NewRoutine_Flaring_B_Day , 0), Savings_IVented_B_Day = IVented_Gas_B_Day - NewIVented_Gas_B_Day, Savings_Tot_B_Day = Savings_Flaring_B_Day + Savings_IVented_B_Day, Savings_Flaring_TCO2e_Day = (Savings_Flaring_B_Day * 0.30180539 * 0.98) + (Savings_Flaring_B_Day * 3.95833079 * 0.02), Savings_IVented_TCO2e_Day = Savings_IVented_B_Day * 3.95833079, Savings_TCO2e_Day = Savings_Flaring_TCO2e_Day + Savings_IVented_TCO2e_Day) %>% drop_na(Oil_Tax_Dollar_B) %>% filter(Gas_Tax_Dollar_B %in% c(-78, -20, -10, -2)) %>% mutate(savings_per = ifelse(Gas_Tax_Dollar_B == -78, '100% (78.63 $/BOE)', ifelse(Gas_Tax_Dollar_B == -20, '90% (20.73 $/BOE)', ifelse(Gas_Tax_Dollar_B == -10, '50% (10.43 $/BOE)', '10% (2.03 $/BOE)'))), savings_per = factor(savings_per, levels = c('10% (2.03 $/BOE)', '50% (10.43 $/BOE)', '90% (20.73 $/BOE)', '100% (78.63 $/BOE)')), GOR_B_B = atan(GOR_B_B), FVR_TCO2e_TCO2e = atan(Savings_IVented_TCO2e_Day/Savings_Flaring_TCO2e_Day)) %>% drop_na(FVR_TCO2e_TCO2e) %>% drop_na(GOR_B_B) #%>% # filter(Oil_Tax_Dollar_B <= 250) %>% # filter(Savings_TCO2e_Day <= 600) viol1 <- ggplot(data = panel_econ, aes(x = '', color = GOR_B_B)) + geom_violin(aes(y = Oil_Tax_Dollar_B), color = '#c2bebe', size = 1, alpha = 0.1, scale = 'width', draw_quantiles = c(0.5)) + geom_point(aes(y = Oil_Tax_Dollar_B), position = position_jitterdodge(0.3, dodge.width = 0.05), alpha = 0.5) + scale_color_gradient(high = '#0000FF', low = '#FF0000', limits = c(0, 1.6), breaks = c(0, 1.6), labels = c('%Oil','%Gas')) + xlab('') + ylab('Oil Tax ($/BBL)') + theme_minimal() + labs(color = '') + # GOR theme(legend.direction = 'horizontal', legend.position = 'bottom', legend.text = element_text(size = 12), text = element_text(size = 12, family = 'helvetica'), legend.spacing.x = unit(2, 'cm'), legend.margin = margin(-50, 0, 0, 0), panel.background = element_rect(colour = 'black', size = 0.5), axis.ticks = element_line(color = 'black')) + 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')) + ggtitle('(A)') + facet_grid(~savings_per) viol1 gor_legend <- get_legend(viol1) gor_legend <- as_ggplot(gor_legend) viol2 <- ggplot(data = panel_econ, aes(x = '', color = FVR_TCO2e_TCO2e)) + geom_violin(aes(y = Savings_TCO2e_Day), color = '#c2bebe', size = 1, alpha = 0.1, scale = 'width', draw_quantiles = c(0.5)) + geom_point(aes(y = Savings_TCO2e_Day), position = position_jitterdodge(0.3, dodge.width = 0.05), alpha = 0.5) + scale_color_gradient(high = '#0070FF', low = '#FF6C00', limits = c(0, 1.6), breaks = c(0, 1.6), labels = c('%Flaring','%Venting')) + xlab('Methane Savings (decline in the natural gas tax)') + ylab('Savings (TCO2e/Day)') + theme_minimal() + labs(color = '') + theme(legend.direction = 'horizontal', legend.position = 'bottom', legend.text = element_text(size = 12), text = element_text(size = 12, family = 'helvetica'), legend.spacing.x = unit(2, 'cm'), legend.margin = margin(-50, 0, 0, 0), strip.background = element_blank(), strip.text = element_blank(), panel.background = element_rect(colour = 'black', size = 0.5), axis.ticks = element_line(color = 'black'), axis.title.x = element_text(vjust = -4.5)) + 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')) + ggtitle('(B)') + facet_grid(~savings_per) viol2 panel_100 <- panel_econ %>% filter(savings_per == '100% (78.63 $/BOE)') scatter2 <- ggplot(data = panel_100, aes(y = Savings_IVented_TCO2e_Day, x = Savings_Flaring_TCO2e_Day, color = atan(GOR_B_B), size = Oil_Tax_Dollar_B)) + geom_point(alpha = 0.5) + scale_color_gradient(high = '#0000FF', low = '#FF0000') + xlab('Routine Flaring Savings Ton CO2e/Day') + ylab('Venting Savings Ton CO2e/Day') + theme_minimal() + theme(legend.position = 'bottom') + facet_grid(~savings_per) scatter2 scatter2_colour <- ggplot_build(scatter2)$data[[1]]$colour panel_radar_100 <- panel_100 %>% select(RE_ID, Savings_IVented_TCO2e_Day, Savings_Flaring_TCO2e_Day, Gas_Tax_Dollar_B, Oil_Tax_Dollar_B) %>% mutate_at(vars(-RE_ID), rescale) %>% mutate(Gas_Tax_Dollar_B = median(Gas_Tax_Dollar_B)) panel_radar_100 <- panel_100 %>% select(RE_ID, Oil_Tax_Dollar_B, Savings_IVented_TCO2e_Day, Savings_Flaring_TCO2e_Day) %>% mutate_at(vars(-RE_ID), rescale) # Assign Colors to GOR panel_100 <- panel_100 %>% mutate(color_gor = scatter2_colour) tax_plt <- ggradar(panel_radar_100, group.line.width = 0.25, group.point.size = 0.5, group.colours = panel_100$color_gor, axis.labels = c('Oil Tax', 'Venting', 'Flaring'), font.radar = 'helvetica', background.circle.colour = 'white', axis.label.size = 5, grid.label.size = 5) + ggtitle('(C)') + theme(legend.position = 'none', plot.title = element_text(size = 15)) + gor_legend + plot_layout(ncol = 1, heights = c(8, 1)) tax_plt figure_2 <- (viol1/viol2)|tax_plt figure_2 ggsave("C:/Users/s15535/Desktop/2 - Energy Economics/2 - Oil & Natural Gas/2.5 - Flaring & Venting/1 - American FV&L/0 - Manuscript/Figure_2.svg", height = 15, width = 30, units = c('cm'), dpi = 400) rm(list = ls(all=TRUE)) ##-----------------------------------------------------------------------------------------------# ## 8. Data Anonimization ##-----------------------------------------------------------------------------------------------# panel_fv <- 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/Second_Stage_Results.csv', sep = ',', header = TRUE, stringsAsFactors = FALSE, encoding = 'UTF-8') panel_ano <- panel_fv %>% select(Geo, Flared_Gas_B_Day, Non_Routine_Flaring_B_Day, Routine_Flaring_B_Day, Vented_Gas_B_Day, IVented_Gas_B_Day, UVented_Gas_B_Day, GOR_SCF_B) %>% sample_frac(0.1) ##-----------------------------------------------------------------------------------------------# ##-----------------------------------------------------------------------------------------------#