library(lubridate); library(ggplot2)
library(ggsci); library(magrittr);
library(dplyr); library(reshape2)
library(inflection)
library(ggnewscale)

theme_set(theme_minimal())
SPR = read.csv("SPR.csv")
SPR$End %<>% mdy()

Find the diffs

SPR$diff = c( -diff(SPR$Reserve), 0)

Original plot with color for diffs.

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

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

Collect the release data

Note: This is important. Read this code!

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
    }
  }
}

We need to clean this up.

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

SPR2 %>% ggplot(aes(x = Start, y = Release)) +
  geom_col() + 
  ylab("Release Size ('000 Barrells)")+
  scale_y_continuous(labels = scales::comma)

Let’s consider the ‘gaps’ in releases. Specfically, 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]
  
  }
#The money plot
SigRel = 2500 #thousands of barrels 

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")

p

#The money plot

library(ggnewscale)

SigRel = 2500 #thousands of barrels 

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("Gas Price by SPR Release")

p

Parking Lot

Ok. Let’s think about the releases that are ‘significant’, by which we mean 1% of capacity:

thresh = 2500

brent %>% filter(relamt > thresh) %>% 
  group_by(relno) %>% 
  mutate(relday = Day - min(Day)) %>%
  ungroup()-> bigrel

Now, find the days to inflection point. {what is an inflection point}? Let’s plot these and see what we can see?

Nos = bigrel$relno %>% unique()

for(i in Nos){
  p = bigrel %>% filter(relno == i) %>% 
    ggplot(aes(x = relday, y = US.Brent)) + geom_point() + 
    ggtitle(i) + geom_smooth()
  
  print(p)
  
}

Ok, let’s find the rightmost inflection points using the ese method!

inflections = vector()

for(i in Nos[Nos>1]){
  
  bigrel %>% filter(relno == i) %>% 
    select(relday, US.Brent) %$% 
    ese(x = relday, y = US.Brent, index = 0) -> infle
  inflections[i] = infle[2]
  
}
Z = data.frame(inflex = inflections[!is.na(inflections)])

Z %>% ggplot(aes(x = 1, y = inflex)) + 
  geom_boxplot() + 
  geom_boxplot(aes(x = 1, y = 129, color = "red")) +
  theme(legend.position = "none") +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) +
  ylab("Days to Inflection Point")

Current Inflection point

bigrel %>% filter(relno == 1) %>% 
    select(relday, US.Brent) %$% 
    ese(x = relday, y = US.Brent, index = 0) -> infle
   infle[2]
## [1] 129