library(lubridate); library(ggplot2)
library(ggsci); library(magrittr);
library(dplyr); library(ggnewscale);
library(knitr) 

theme_set(theme_minimal())

Data on the Strategic Petroeum Reserve (SPR) levels was obtained from the U.S. Government website, https://www.eia.gov/dnav/pet/hist_xls/WCSSTUS1w.xls

Data for fuel prices was obtained from the U.S. Government website, https://www.eia.gov/dnav/pet/hist/LeafHandler.ashx?n=PET&s=RBRTE&f=M

Both files were downloaded on 14 November 2022.

SPR = read.csv("SPR.csv")
SPR$End %<>% mdy()
SPR$diff = c( -diff(SPR$Reserve), 0)

Original plot with color for diffs.

SPR$Decline = FALSE
SPR$Decline[SPR$diff < -10 ] = TRUE

p = SPR %>% 
  ggplot(aes(x = End, y = Reserve, color = Decline)) +
  geom_point() + scale_color_aaas() + 
  theme(legend.position = "none") + 
  scale_y_continuous(name = "Reserve, ('000 Barrels)",
                     labels = scales::comma) + 
  ggtitle("US Strategic Petroleum Reserves") +
  xlab("")

p

Collect the release data

SPR$Rel = 0
release = 1
Inrelease = TRUE
for(i in 1:dim(SPR)[1]){
  if(Inrelease){
    if(SPR$diff[i] < 0){
    #Case where we're in a release and it's ongoing
    SPR$Rel[i] = release #Current one
    } else{
    SPR$Rel[i] = 0
    Inrelease = FALSE
    release = release + 1 #increment here
    }
    }else{
    if(SPR$diff[i] < 0){
      SPR$Rel[i] = release
      Inrelease = TRUE
    }else{
      SPR$Rel[i] = 0
    }
  }
}

Tidy this up a bit:

SPR %>% filter(Rel > 0) %>%
  group_by(Rel) %>% 
  summarize(Start = min(End),
            End = max(End),
            Release = -sum(diff)) -> SPR2

Let’s consider the ‘gaps’ in releases. Specifically, if two releases are within 65 days, we will consider them the ‘same’ release.

daysBetween = 65
SPR %>% filter(Rel > 0) %>% 
  group_by(Rel) %>% 
  summarize(start = min(End), 
            end = max(End)) ->
  Stpdates
Stpdates$CorRel = 0
Stpdates$days = 0
this.rel = 1
for(i in 1:(dim(Stpdates)[1]-1)){
  int = (Stpdates$start[i] - 
           Stpdates$end[i+1]) %>% 
    as.numeric()
  Stpdates$days[i] = int
  if(int <daysBetween){
    Stpdates$CorRel[i] = this.rel
  }else{
    this.rel = this.rel + 1
    Stpdates$CorRel[i] = this.rel
  }
}

Bring in the brent crude data, as well as the gasoline prices (Both from source data)

brent = read.csv("Brent.csv")
brent$Date %<>% dmy()
brent %<>% filter(US.Brent > 0)
brent$Day = as.numeric(brent$Date - min(brent$Date))
gas = read.csv("Gas.csv")
gas$Date %<>% dmy()
gas %<>% filter(NY > 0)
gas$Day = as.numeric(gas$Date - min(gas$Date))
brent$isrelease = 0
for(i in 1:dim(SPR2)[1]){
  brent$isrelease[brent$Date %in% seq(SPR2$Start[i],
                                      to = SPR2$End[i],
                                      by = "day")] = 1
}
gas$isrelease = 0
for(i in 1:dim(SPR2)[1]){
  gas$isrelease[gas$Date %in% seq(SPR2$Start[i],
                                  to = SPR2$End[i],  
                                  by = "day")] = 1
}

Let’s do prices by release.

brent$relno = 0
brent$relamt = 0
for(i in 1:dim(SPR2)[1]){
  brent$relno[brent$Date %in% seq(SPR2$Start[i],
                                  to = SPR2$End[i],
                                  by = "day")] = SPR2$Rel[i]

  brent$relamt[brent$Date %in%
               seq(SPR2$Start[i], 
               to = SPR2$End[i],by = "day")] = 
    SPR2$Release[i]
  }

gas$relno = 0
gas$relamt = 0
for(i in 1:dim(SPR2)[1]){
  gas$relno[gas$Date %in% seq(SPR2$Start[i],
                                  to = SPR2$End[i],
                                  by = "day")] = SPR2$Rel[i]
  gas$relamt[gas$Date %in%
               seq(SPR2$Start[i], 
               to = SPR2$End[i],by = "day")] = 
    SPR2$Release[i]
  }
SigRel = 2500 #2.5 million barrel threshold  
brent %>% filter(relno > 0) %>% 
  filter(relamt > SigRel) %>% 
  group_by(relno) %>%
  mutate(relday = Day - min(Day)) %>%
  ungroup() -> pdat
pdat$isOne = 0
pdat$isOne[pdat$relno == 1] = 1
p = ggplot() + geom_point(data = pdat %>% filter(isOne == 0),
                          aes(x = relday, y= US.Brent, color = (relno))) +
  new_scale_color() + 
  geom_point(data = pdat %>% filter(isOne == 1),
             aes(x = relday, y= US.Brent, color = relno)) +
  scale_color_gradient(low = "red", high = "red") +
  theme(legend.position = "none") +
  scale_y_continuous(labels = scales::dollar,
                     name = "US Brent Crude Price") +
  xlab("Days since Release Start") + 
  ggtitle("US Brent Crude Prices by Release Day")

p

gas %>% filter(relno > 0) %>% 
  filter(relamt > SigRel) %>% 
  group_by(relno) %>%
  mutate(relday = Day - min(Day)) %>%
  ungroup() -> gdat
gdat$isOne = 0
gdat$isOne[gdat$relno == 1] = 1
p = ggplot() + 
  geom_point(data = gdat %>% filter(isOne == 0),
             aes(x = relday, y= NY, color = (relno))) +
  new_scale_color() +
  geom_point(data = gdat %>% filter(isOne == 1),
             aes(x = relday, y= NY, color = relno)) +
  scale_color_gradient(low = "orange", high = "orange") +
  theme(legend.position = "none") +
  scale_y_continuous(labels = scales::dollar,
                     name = "New York Gallon Price") + 
  xlab("Days since Release Start") + 
  ggtitle("Gasoline Price by SPR Release")
p

Another View

p = pdat %>% 
  group_by(relno) %>% 
  summarize(Days = as.numeric(max(Date) -
                                min(Date)), 
            amt = max(relamt), 
            PriceDiff = max(US.Brent) -
              min(US.Brent), 
            isOne = max(isOne), 
            start = year(min(Date))) %>%
  mutate(Impact = 1000*PriceDiff / amt) %>%
  ggplot(aes(x = amt, y= PriceDiff,
             label = start)) + 
  geom_label(check_overlap = TRUE) +
  scale_x_continuous(name = 
                       "Release Amount (total, '000)",
                     labels = scales::comma) +
  scale_y_continuous(name = "Price differential",
                     labels = scales::dollar)
p

alternate graph:

 p = pdat %>% 
  group_by(relno) %>%
  summarize(Days =
              as.numeric(max(Date) -
                           min(Date)),
            amt = max(relamt), 
            PriceDiff = max(US.Brent) -
              min(US.Brent), 
            isOne = max(isOne), 
            start = year(min(Date))) %>%
  mutate(Impact = 1000*PriceDiff / amt)  %>% 
  ggplot(aes(x = start,
             y = -PriceDiff, fill = Days))  + 
  geom_col(position = "dodge") + 
  scale_color_continuous(low = "blue", high = "red") + 
  theme(legend.position = "none") + 
  scale_y_continuous(name =
                       "Price differential",
                     labels = scales::dollar) +
  xlab("Release Start Date") + scale_x_continuous(position = "top")

p

Data Formats

Below are samples of the data formats downloaded from the eia websites.

SPR

read.csv("SPR.csv") %>% head() %>% kable()
End Reserve
8/1/2022 445057
7/1/2022 468006
6/1/2022 493324
5/1/2022 523109
4/1/2022 547866
3/1/2022 566061

Crude and Gasoline Prices

read.csv("Gas.csv") %>% head() %>% kable()
Date NY GC
2-Jun-86 0.468 0.445
3-Jun-86 0.436 0.418
4-Jun-86 0.418 0.398
5-Jun-86 0.431 0.415
6-Jun-86 0.421 0.403
9-Jun-86 0.409 0.398 
read.csv("brent.csv") %>% head() %>% kable()
Date US.Brent EUR.Brent
2-Jan-86 25.56 NA
3-Jan-86 26.00 NA
6-Jan-86 26.53 NA
7-Jan-86 25.85 NA
8-Jan-86 25.87 NA
9-Jan-86 26.03 NA