Uploaded Profit Assessment

SOE-response-memo/SOE_Profit.Rmd

@@ -0,0 +1,318 @@
+---
+title: "Profit_Assessment"
+author: "Geret DePiper"
+date: "`r Sys.Date()`"
+
+bibliography: ["Profit_MAFMC.bib"]
+link-citations: yes
+output:
+  bookdown::pdf_document2:
+    toc: False
+    fig_caption: yes
+    keep_tex: yes
+---
+
+```{r setup, include=FALSE}
+knitr::opts_chunk$set(echo = FALSE)
+knitr::opts_chunk$set(warning = FALSE)
+knitr::opts_chunk$set(comment = FALSE)
+knitr::opts_chunk$set(error = FALSE)
+knitr::opts_chunk$set(message = FALSE)
+
+library(dplyr)
+library(tidyr)
+library(dbplyr)
+library(ggplot2)
+library(here)
+library(data.table)
+library(plyr)
+library(reshape2)
+library(ggpubr)
+library(urca)
+library(kableExtra)
+load(file="F:/MAFMC Risk Assessment/Profit/Net_revenue_data_CAMS_LAND")
+```
+
+## Introduction
+
+The SSC asked for a refined indicator of net revenue to better proxy for profitability. In what follows I provide an overview of the analysis conducted, per SSC recommendations. Specifically, I present a net revenue assessment (gross revenue minus trip costs, for federally reporting trips), profitability ratio (gross revenue for all trips divided by trip costs), and cointegration analysis of gross revenue and fuel prices, by ecological production unit between the years 2000 - 2023. The aim of this document is to provide information sufficient for the SSC to select a preferred indicator to present within the MAFMC's State of the Ecosystem Report and Ecosystem Approach to Fishery Management Risk Assessment. All dollar values are reported in $2023 real values.
+
+Two important notes on the analysis: The underlying data does not currently include information for trips in which there is a valid trip report but no corresponding dealer report. In addition, the algorithm currently employed to assign trips to region fished differs slightly than the algorithm used in the formal State of the Ecosystem Report. Both of these issues will look to be addressed in the coming year. This analysis will therefore differ from the final indicators developed, but should suffice for decision-making with respect to the overarching methodologies to employ.
+
+## Cost Coverage
+
+As background, Figure \@ref(fig:coverage) presents an overview of gross revenue from federal trips, for which costs can be calculated, versus gross revenue for which trip costs cannot be computed. A number of issues are highlighted by the graphs. First, the coverage varies substantially by ecological production unit, with Georges Bank having highest percentage of gross revenue covered. The Mid-Atlantic Bight shows substantial interannual variability in coverage, with costs covering between 50\% and 75\% of total gross revenue generated in most recent years. Second, and in particular for the Mid-Atlantic, the gross revenue for federally permitted trips presents different trends than the the gross revenue by non-federally permitted trip. This suggests different production functions between the two fleet segments, and thus inappropriateness of imputing costs from a model based on federally permitted trips onto non-federally permitted trips. 
+
+```{r coverage, fig.cap="Gross revenue grouped by federally permitted and non-federally permitted trips, by region.", warnings=FALSE, comment=FALSE,error=FALSE, message=FALSE}
+
+Cost_Coverage <- CAMS_Trip_Revenue %>%
+  mutate(Missing= "Net Revenue Possible",
+    Missing= ifelse(is.na(Real_Cost),"Missing Costs", Missing)) %>%
+  group_by(EPU,Missing,YEAR) %>%
+  dplyr::summarise(Total_revenue=sum(Real_Revenue,na.rm = TRUE)) %>% ungroup()
+
+ggplot(Cost_Coverage,aes(x=YEAR,y=Total_revenue/1000000,colour=Missing))+
+  geom_line(stat="identity")+facet_wrap(vars(EPU), scales="free_y")+ expand_limits(y = 0)+ theme(axis.text.x = element_text(angle=90, hjust=1))+labs(y="Total Revenue ($2023 Million)",x="Year")
+```
+
+## Net Revenue
+
+Figure \@ref(fig:coverage) presents the summed net revenue estimates for federally permitted trips. Again, there is substantial heterogeneity and variability in net revenue across regions. Apart from a slight uptick between 2019 - 2020, there has been a steady decline in net revenue within the Mid-Atlantic since 2016, with 2023 representing the lowest net revenue estimate of the series.
+
+```{r netrev, echo=FALSE, fig.cap="Net revenue generated by federally permitted trips."}
+Net_revenue <- CAMS_Trip_Revenue %>% 
+  filter(!is.na(VTR_TRIPID)) %>%
+  group_by(YEAR, EPU) %>%
+    dplyr::summarise(Net_Revenue =sum(Net_Revenue, na.rm=TRUE)) %>% ungroup()
+
+#Useful plots
+
+ggplot(Net_revenue, aes(x=YEAR,y=Net_Revenue/1000000))+
+  geom_line()+facet_wrap(vars(EPU), scales="free_y")+ expand_limits( y = 0)+labs(y="Net Revenue ($2023 Million)",x="Year")
+
+```
+```{r netrev_violin, echo=FALSE, include=FALSE}
+ggplot(CAMS_Trip_Revenue %>%
+  filter(!is.na(VTR_TRIPID)), aes(x=as.factor(YEAR),y=Net_Revenue, fill=as.factor(YEAR)))+
+  geom_violin(show.legend = FALSE)+facet_wrap(vars(EPU), scales="free_y")+ theme(axis.text.x = element_text(angle=45, hjust=1))
+```
+Although net revenue is the most-oft assessed metric of profitability in fisheries management, the fact that 25-50\% of trip revenue is excluded from this analysis has been a serious deficiency pointed out by both the SSC and Council, particularly through its Ecosystem and Ocean Planning Committee. As such, two additional analyses are presented in what follows. First, we present a "profitabilty index", which in this case is a ratio of gross revenue to estimated costs. This indicator is common in the benchmarking literature [@balk_residual_2003], despite the ratio not being unique across combinations of revenue and costs. 
+
+In this application, we generate three separate indexes of performance: Revenue, Cost, and Profit. The Revenue index is the sum of gross revenue from both federally permitted and non-federally permitted trips by region, normalized by the arithmetic mean gross revenue across regions in 2000, the first year of the series. The Cost index is the geometric mean of the normalized trip-level costs by region, with normalization based on the arithmetic average cost across all federally permitted trips. The trip level normalization ensures the index can be compared across regions, abstracting away from the issue of scale and downweighting outliers, as is standard in the literature [e.g. @walden_profits_2022]. The final Profit index is the ratio of Revenue and Cost indices. The implicit assumption embedded in the Cost and Profit indices is that average costs in the federally permitted fleet adequately proxies for non-permitted fleet costs, a tenuous assumption given the trends seen in figure \@ref(fig:coverage). 
+
+Each of these indices are presented in figure \@ref(fig:Profitability). For the Revenue index, all three regions present index numbers below 1 in 2023, indicating that revenue is below the 2000 base year. In the Mid-Atlantic, the revenue index is at a historical low, consistent with the gross numbers presented in figure \@ref(fig:coverage). The cost index for all three regions measure betweeen 1.5 (50\% higher than 2000) in the Gulf of Maine to 2.25 (125\% higher than 2000) in Georges Bank. Overall, we see costs have been consistently higher than the 2000 benchmark since 2021.
+
+The Profit index across all three regions is near historic lows, driven by both historic low revenue in the Gulf of Maine and Mid-Atlantic, coupled by relatively high costs across all three regions.
+
+```{r Profitability, echo=FALSE}
+CI <- function(data, year, cost_ref) { #COST INDEX
+  data <- data[,c("YEAR","CAMSID","Real_Cost","EPU")] %>%
+    mutate(index=Real_Cost/as.numeric(cost_ref))
+
+  cost_index<- setDT(data)[,. (gmean=exp(mean(log(index), na.rm=TRUE))), by=c("YEAR","EPU")] 
+  cost_base=cost_index[(cost_index$YEAR==year),] %>% 
+  select(gmean,EPU)
+colnames(cost_base) <- c("cost_base","EPU")
+  cost_index <- cost_index %>%
+    left_join(cost_base, by="EPU") %>%
+    mutate(cindex=round(gmean/cost_base,3)) %>% select(YEAR,EPU,cindex) %>%
+    return()
+  
+ }
+
+RI<-function(data,base_rev,base_year){ #Revenue Index
+  
+  # sum_rev<- data %>% filter(YEAR>=base_year) %>%
+  #   group_by("YEAR","EPU") %>%
+  #   mutate(trev=sum(Real_Revenue, na.rm=TRUE)) %>%
+  #   ungroup()
+  
+  data<-data[(data$YEAR>=base_year),]
+  sum_rev<-setDT(data)[,.(trev=sum(Real_Revenue, na.rm=TRUE)),
+                             by=c("YEAR","EPU")]
+  colnames(sum_rev)[1]<-"YEAR"
+  colnames(sum_rev)[2]<-"EPU"
+   
+  sum_rev$index=round(sum_rev$trev/base_rev,3)
+  denom=as.numeric(unlist(sum_rev[(sum_rev$YEAR==base_year),"index"]))
+  sum_rev$rindex=round(sum_rev$index/denom,3)
+  sum_rev<-sum_rev[,c("YEAR","EPU","rindex")]
+  
+  return(sum_rev)
+}
+PI<-function(REV_I, COST_I){
+  Index_fin<-join(REV_I,COST_I, by=c("YEAR","EPU"), type="inner")
+  Index_fin$prof_index=Index_fin$rindex/Index_fin$cindex
+  Index_fin<-Index_fin[,c("YEAR","EPU","rindex","cindex","prof_index")]
+  
+  return(Index_fin)
+}
+PLOTS<-function(PI,AREA,START,END){
+    Index_for_graph<-melt(PI,id="YEAR")
+    
+    labels=c('Revenue Index', 'Cost Index', 'Profitability Index')
+    title1<-paste0(AREA," Revenue, Cost and Profitability Indices ", START," -",END)
+    subtitle1<-paste0("Base Year ",START)
+  
+    PI_PLOT<- ggplot(Index_for_graph, aes(x = YEAR, y = value)) +
+      geom_line(aes(color=variable))+
+      labs(title=title1,
+           subtitle=subtitle1,
+           x="Year", y="Index Value")+
+      theme(plot.title=element_text(hjust=0.5),
+            plot.subtitle=element_text(hjust=0.5))+
+      scale_color_discrete(labels=labels)
+    
+}
+
+PI_PLOT<-function(PI,START,END){
+  Index_for_graph<-melt(PI,id="YEAR")
+  
+  labels=c('GB', 'GOM', 'MAB')
+  title1<-paste0(" Profitability Indices by Year and Area ", START,"-",END)
+  subtitle1<-paste0("Base Year ",START)
+  
+  PI_PLOT<- ggplot(Index_for_graph, aes(x = YEAR, y = value)) +
+    geom_line(aes(color=variable))+
+    labs(title=title1,
+         subtitle=subtitle1,
+         x="Year", y="Index Value")+
+    theme(plot.title=element_text(hjust=0.5),
+          plot.subtitle=element_text(hjust=0.5))+
+    scale_color_discrete(labels=labels)
+  
+}
+
+start=min(CAMS_Trip_Revenue$YEAR)
+terminal=max(CAMS_Trip_Revenue$YEAR)
+ref_year_all<-subset(CAMS_Trip_Revenue, YEAR==start)
+# avg_cost_all <- CAMS_Trip_Revenue %>% filter(YEAR==start) %>%
+#   group_by(EPU) %>%
+#   dplyr::summarise(cost_ref=mean(Real_Cost, na.rm=TRUE)) %>% ungroup()
+
+avg_cost_all <- CAMS_Trip_Revenue %>% filter(YEAR==start) %>%
+  dplyr::summarise(cost_ref=mean(Real_Cost, na.rm=TRUE))
+
+#Cost Index CAMS_Trip_Revenue
+CostIndex<-CI(CAMS_Trip_Revenue%>%
+  filter(EPU!="Other") ,start,avg_cost_all)
+
+mean_rev <- CAMS_Trip_Revenue %>% filter(YEAR==start) %>%
+  filter(EPU!="Other") %>%
+  group_by(EPU) %>% 
+  mutate(trev=sum(Real_Revenue, na.rm=TRUE)) %>%
+  ungroup() %>%
+  dplyr::summarise(trev=mean(trev, na.rm=TRUE)) %>% 
+  unlist() %>% as.numeric()
+
+RevenueIndex <- RI(CAMS_Trip_Revenue %>%
+  filter(EPU!="Other") ,mean_rev,start)
+
+PI<-PI(RevenueIndex,CostIndex)
+
+p1 <- ggplot(PI,aes(x=YEAR,y=rindex, colour=EPU))+geom_line()
+p2 <- ggplot(PI,aes(x=YEAR,y=cindex, colour=EPU))+geom_line()
+p3 <- ggplot(PI,aes(x=YEAR,y=prof_index, colour=EPU))+geom_line()
+
+ggarrange(p1, p2, p3, ncol=1, common.legend = TRUE, legend="bottom")
+
+```
+In order to assess how well costs on federally permitted trips proxy for non-federally permitted trips, we conduct a cointegration analysis, which captures the extent to which two variables move concurrently across time. In this case, we look at whether the gross revenue from non-federally permitted trips move in syncopation with fuel prices in the region, particularly Diesel no. 2 spot prices from the port of New York, downloaded from the Department of Energy on February 7, 2025 (https://www.eia.gov/dnav/pet/hist/LeafHandler.ashx?n=PET&s=EER_EPD2F_PF4_Y35NY_DPG&f=W
+). We then compare this analysis with one conducted on federally permitted trips, to assess whether differences exist between the two fleet segments. 
+
+```{r stationarity, echo=FALSE}
+Diesel <- read.csv(file=here::here("Weekly_New_York_Harbor_No._2_Heating_Oil_Spot_Price_FOB.csv"),
+                     skip = 5, col.names=c("Date","Price"), header=FALSE) %>%
+  mutate(Date = as.Date(Date,"%m/%d/%Y"),
+         Week = strftime(Date, 
+                         format = "%V"),
+         YEAR = as.numeric(strftime(Date,
+                         format="%Y")))
+
+Cointegration <- CAMS_Trip_Revenue %>%
+  mutate(Missing= "Net Revenue Possible",
+    Missing= ifelse(is.na(Real_Cost),"Missing Costs", Missing),
+    Week = strftime(DATE_TRIP, 
+                         format = "%V")) %>%
+  left_join(Diesel) %>%
+  mutate(Price=Price/value) %>%
+   filter(!is.na(Price)) %>% 
+  group_by(Date,Missing) %>%
+  dplyr::summarise(Real_Revenue=sum(Real_Revenue,na.rm=TRUE),
+            Price=mean(Price, na.rm=TRUE)) %>% ungroup
+
+Missing_cost <- Cointegration %>% filter(Missing=="Missing Costs")
+Nonmissing_cost <- Cointegration %>% filter(Missing=="Net Revenue Possible")
+
+Stationarity1 <- ur.df(Missing_cost$Real_Revenue/100000000,
+                       type="none",
+                       lags = 15,
+                       selectlags="AIC")
+p1 <- ""
+  p1 <- ifelse(Stationarity1@teststat<Stationarity1@cval[3]," * ",p1)
+  p1 <- ifelse(Stationarity1@teststat<Stationarity1@cval[2]," * * ",p1)
+  p1 <- ifelse(Stationarity1@teststat<Stationarity1@cval[1]," * * * ",p1)
+#summary(Stationarity1)
+
+Stationarity2 <- ur.ers(Missing_cost$Real_Revenue/100000000,
+                       lag.max = 15)
+p2 <- ""
+  p2 <- ifelse(Stationarity2@teststat<Stationarity2@cval[3],
+               " * ",p2)
+  p2 <- ifelse(Stationarity2@teststat<Stationarity2@cval[2],
+               " * * ",p2)
+  p2 <- ifelse(Stationarity2@teststat<Stationarity2@cval[1],
+               " * * * ",p2)
+#summary(Stationarity2)
+
+Stationarity3 <- ur.df(Nonmissing_cost$Price,
+                       type="none",
+                       lags = 15,
+                       selectlags="AIC")
+p3 <- ""
+  p3 <- ifelse(Stationarity3@teststat<Stationarity3@cval[3],
+               " * ",p3)
+  p3 <- ifelse(Stationarity3@teststat<Stationarity3@cval[2],
+               " * * ",p3)
+  p3 <- ifelse(Stationarity3@teststat<Stationarity3@cval[1],
+               " * * * ",p3)
+#summary(Stationarity3)
+
+Stationarity4 <- ur.ers(Nonmissing_cost$Price,
+                       lag.max = 15)
+p4 <- ""
+  p4 <- ifelse(Stationarity4@teststat<Stationarity4@cval[3],
+               " * ",p4)
+  p4 <- ifelse(Stationarity4@teststat<Stationarity4@cval[2],
+               " * * ",p4)
+  p4 <- ifelse(Stationarity4@teststat<Stationarity4@cval[1],
+               " * * * ",p4)
+#summary(Stationarity4)
+
+Stationarity_results <- rbind(
+  c("Gross Revenue",[email protected],Stationarity1@lags,paste0(Stationarity1@teststat,p1, sep="")),
+c("Gross Revenue",[email protected],"15",paste0(Stationarity2@teststat,p2, sep="")))
+
+Stationarity_results <- rbind(Stationarity_results,
+c("Diesel Price",[email protected],Stationarity3@lags,paste0(Stationarity3@teststat,p3, sep="")))
+
+Stationarity_results <- rbind(Stationarity_results,
+c("Diesel Price",[email protected],"15",paste0(Stationarity4@teststat,p4, sep="")))
+
+knitr::kable(Stationarity_results, format = "latex", linesep = "",
+      col.names = c("Variable","Model","Lags","tau"),
+      caption = "Unit root tests for Gross Revenue and Diesel prices.",
+      #align = aca,
+      booktabs = T) %>%
+  row_spec(0, bold = TRUE) %>% 
+  add_footnote(" * * * \\textit{p}-value < 0.01,* * \\textit{p}-value < 0.05, * \\textit{p}-value < 0.1.",notation="symbol",threepart=TRUE)
+```
+
+```{r cointegration, echo=FALSE}
+Cointegration_1 <- ur.df(Missing_cost$Real_Revenue/100000000-Missing_cost$Price, 
+      lags = 15, 
+      selectlags = "AIC", 
+      type = "drift")
+summary(Cointegration_1)
+
+Cointegration_2 <-ur.ers(Missing_cost$Real_Revenue/100000000-Missing_cost$Price,
+       model = "constant", 
+       lag.max = 15)
+summary(Cointegration_2)
+
+
+Cointegration_3 <- ur.df(Nonmissing_cost$Real_Revenue/100000000-Nonmissing_cost$Price, 
+      lags = 15, 
+      selectlags = "AIC", 
+      type = "drift")
+summary(Cointegration_3)
+
+Cointegration_4 <-ur.ers(Missing_cost$Real_Revenue/100000000-Missing_cost$Price,
+       model = "constant", 
+       lag.max = 15)
+summary(Cointegration_4)
+
+```
+
+# References
+