Oil Price & Trend

chart_oil_price.r

#!/usr/bin/env Rscript
# Crude Oil (WTI) - Log-Trend & Super-Cycle Analysis
# Charts: 1_ Price+Trend, 2_ Deviation+Supertrend, 3_ Normalized
# Data: FRED DCOILWTICO daily spot price, monthly last observation

library(tidyverse)
library(ggrepel)
library(zoo)

sf <- 1.5
width <- 600 * sf
height <- 335 * sf
options(vsc.dev.args = list(width = width, height = height, res = 72 * sf))

LABEL <- "Crude Oil WTI Spot Price"
PREFIX <- "chart_crude_oil_wti"
START <- as.Date("1986-01-01")
FRED_SERIES <- "DCOILWTICO"
FRED_URL <- paste0("https://fred.stlouisfed.org/graph/fredgraph.csv?id=", FRED_SERIES)
LOCAL_CSV <- paste0(FRED_SERIES, ".csv")

RED <- "#C62828"; ORANGE <- "#E65100"; AMBER <- "#F9A825"
GREEN <- "#2E7D32"; BLUE <- "#1565C0"; GRAY <- "#757575"; TEAL <- "#00695C"

SUPER_CYCLE_YEARS_DEFAULT <- 12
LOESS_SPAN_MULT <- 1.8
LOESS_MIN_SPAN <- 0.18
LOESS_MAX_SPAN <- 0.98
EDGE_DAMPEN_YEARS <- 4
EDGE_MIN_WEIGHT <- 0.35
Z_BAND <- 1.96

if (file.exists(LOCAL_CSV)) {
  cat("-> Reading local FRED series", LOCAL_CSV, "...\n")
  raw <- read_csv(LOCAL_CSV, col_types = cols(.default = "c"), show_col_types = FALSE)
} else {
  cat("-> Downloading FRED series", FRED_SERIES, "...\n")
  raw <- read_csv(FRED_URL, col_types = cols(.default = "c"), show_col_types = FALSE)
}

df <- raw |>
  rename(date = observation_date, close = all_of(FRED_SERIES)) |>
  mutate(
    date = as.Date(date),
    close = suppressWarnings(as.numeric(close))
  ) |>
  filter(date >= START, !is.na(close), close > 0) |>
  mutate(month = as.yearmon(date)) |>
  group_by(month) |>
  slice_max(date, n = 1, with_ties = FALSE) |>
  ungroup() |>
  transmute(date = as.Date(month), close) |>
  arrange(date)

if (nrow(df) < 60) {
  stop("Not enough usable price history after parsing FRED data.")
}

cat(" Final:", nrow(df), "months,", format(min(df$date)), "to", format(max(df$date)), "\n")

df <- df |> mutate(date_num = as.numeric(date), log_close = log(close))
fit <- lm(log_close ~ date_num, data = df)
cagr <- (exp(coef(fit)[2] * 365.25) - 1) * 100

df <- df |>
  mutate(
    trend = exp(predict(fit, newdata = pick(everything()))),
    pct_dev = (close - trend) / trend * 100
  )

cat(" Trend CAGR:", round(cagr, 2), "%/yr\n")

data_years <- as.numeric(difftime(max(df$date), min(df$date), units = "days")) / 365.25
cycle_years <- if (data_years >= 20) SUPER_CYCLE_YEARS_DEFAULT else max(5, round(data_years * 0.35))
k_months <- as.integer(cycle_years * 12)
n <- nrow(df)
span <- min(LOESS_MAX_SPAN, max(LOESS_MIN_SPAN, LOESS_SPAN_MULT * k_months / n))
t_idx <- seq_len(n)

edge_n <- min(as.integer(EDGE_DAMPEN_YEARS * 12), floor(n / 2))
w <- rep(1, n)
if (edge_n >= 2) {
  ramp <- seq(EDGE_MIN_WEIGHT, 1, length.out = edge_n)
  w[seq_len(edge_n)] <- ramp
  w[(n - edge_n + 1):n] <- rev(ramp)
}

lo <- loess(df$pct_dev ~ t_idx, span = span, degree = 1,
            family = "symmetric", weights = w,
            control = loess.control(surface = "direct"))
df$supertrend <- as.numeric(predict(lo, data.frame(t_idx = t_idx)))
df$norm_dev <- df$pct_dev - df$supertrend

cat(" Supertrend:", cycle_years, "yr cycle, span =", round(span, 3), "\n")

decade_labels <- function(dates, values, threshold = 60) {
  tbl <- tibble(date = dates, val = values) |>
    mutate(decade = floor(lubridate::year(date) / 10) * 10) |>
    group_by(decade) |>
    summarize(
      max_d = date[which.max(val)], max_v = max(val),
      min_d = date[which.min(val)], min_v = min(val),
      .groups = "drop"
    )
  bind_rows(
    tbl |> filter(abs(max_v) >= threshold) |> transmute(date = max_d, val = max_v),
    tbl |> filter(abs(min_v) >= threshold) |> transmute(date = min_d, val = min_v)
  ) |>
    distinct(date, .keep_all = TRUE)
}

last_str <- format(max(df$date), "%b %Y")
cur <- tail(df, 1)

cat("-> Chart 1: Price + Trend\n")
extrema1 <- bind_rows(
  decade_labels(df$date, df$pct_dev, 65),
  tibble(date = cur$date, val = cur$pct_dev)
) |>
  distinct(date, .keep_all = TRUE) |>
  left_join(df |> select(date, close), by = "date")

p1 <- ggplot(df, aes(x = date)) +
  geom_ribbon(aes(ymin = pmin(close, trend), ymax = trend), fill = GREEN, alpha = 0.15) +
  geom_ribbon(aes(ymin = trend, ymax = pmin(pmax(close, trend), trend * 1.75)), fill = AMBER, alpha = 0.12) +
  geom_ribbon(aes(ymin = trend * 1.75, ymax = pmax(close, trend * 1.75)), fill = RED, alpha = 0.18) +
  geom_line(aes(y = close), color = BLUE, linewidth = 0.5) +
  geom_line(aes(y = trend), color = "black", linetype = "dashed", linewidth = 0.5) +
  geom_point(data = extrema1, aes(y = close), color = ORANGE, size = 2) +
  geom_text_repel(data = extrema1, aes(y = close, label = sprintf("%+.0f%%", val)),
                  color = BLUE, size = 2.8, max.overlaps = 15, seed = 42) +
  scale_x_date(date_breaks = "5 year", date_labels = "%Y") +
  scale_y_continuous(trans = "log2", labels = scales::label_dollar()) +
  labs(
    title = paste0(LABEL, " (Log Scale)"),
    subtitle = paste0("Since ", format(min(df$date), "%Y"), " | Log-trend CAGR ~ ", round(cagr, 1), "%/yr"),
    x = NULL, y = "USD per barrel",
    caption = "Source: FRED DCOILWTICO. Monthly values use the last daily observation in each month."
  ) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 30, hjust = 1), legend.position = "none")

ggsave(paste0(PREFIX, "_1_price.png"), p1,
       width = width / 72 / sf, height = height / 72 / sf, dpi = 72 * sf)

cat("-> Chart 2: Deviation + Supertrend\n")
lbl2 <- bind_rows(
  decade_labels(df$date, df$pct_dev, 65),
  tibble(date = cur$date, val = cur$pct_dev)
) |>
  distinct(date, .keep_all = TRUE)

p2 <- ggplot(df, aes(date, pct_dev)) +
  geom_ribbon(aes(ymin = pmin(pct_dev, 0), ymax = 0), fill = GREEN, alpha = 0.15) +
  geom_ribbon(aes(ymin = 0, ymax = pmin(pmax(pct_dev, 0), 75)), fill = AMBER, alpha = 0.12) +
  geom_ribbon(aes(ymin = 75, ymax = pmax(pct_dev, 75)), fill = RED, alpha = 0.18) +
  geom_hline(yintercept = 0, linetype = "dashed", linewidth = 0.4) +
  geom_line(color = BLUE, linewidth = 0.4, alpha = 0.6) +
  geom_line(aes(y = supertrend), color = TEAL, linewidth = 0.8) +
  geom_point(data = lbl2, aes(y = val), color = ORANGE, size = 1.8) +
  geom_text_repel(
    data = lbl2 |> mutate(lbl = if_else(date == max(date), paste0(last_str, ": ", sprintf("%+.0f%%", val)), sprintf("%+.0f%%", val))),
    aes(y = val, label = lbl),
    color = BLUE, size = 2.8, fontface = "bold",
    nudge_y = 5, box.padding = 0.5, max.overlaps = 15, seed = 42
  ) +
  scale_x_date(date_breaks = "5 year", date_labels = "%Y") +
  scale_y_continuous(labels = \(x) paste0(x, "%")) +
  labs(
    title = paste0(LABEL, ": Deviation from Log-Trend + Supertrend"),
    subtitle = paste0("Trend CAGR ~ ", round(cagr, 1), "%/yr | Teal = ", cycle_years, "yr supertrend"),
    x = NULL, y = "Deviation from trend (%)",
    caption = "Source: FRED DCOILWTICO. Supertrend: robust LOESS with edge down-weighting."
  ) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 30, hjust = 1), legend.position = "none")

ggsave(paste0(PREFIX, "_2_deviation.png"), p2,
       width = width / 72 / sf, height = height / 72 / sf, dpi = 72 * sf)

cat("-> Chart 3: Double-Normalized\n")
norm_sd <- sd(df$norm_dev, na.rm = TRUE)
band <- Z_BAND * norm_sd

lbl3 <- bind_rows(
  decade_labels(df$date, df$norm_dev, 45),
  tibble(date = cur$date, val = cur$norm_dev)
) |>
  distinct(date, .keep_all = TRUE)

p3 <- ggplot(df, aes(date, norm_dev)) +
  geom_ribbon(aes(ymin = pmin(norm_dev, 0), ymax = 0), fill = GREEN, alpha = 0.15) +
  geom_ribbon(aes(ymin = 0, ymax = pmin(pmax(norm_dev, 0), band)), fill = AMBER, alpha = 0.12) +
  geom_ribbon(aes(ymin = band, ymax = pmax(norm_dev, band)), fill = RED, alpha = 0.18) +
  geom_hline(yintercept = 0, linetype = "dashed", linewidth = 0.4) +
  geom_hline(yintercept = c(-band, band), color = GRAY, linetype = "dotted", linewidth = 0.35) +
  geom_line(color = BLUE, linewidth = 0.5) +
  geom_point(data = lbl3, aes(y = val), color = ORANGE, size = 1.8) +
  geom_text_repel(
    data = lbl3 |> mutate(lbl = if_else(date == max(date), paste0(last_str, ": ", sprintf("%+.0f%%", val)), sprintf("%+.0f%%", val))),
    aes(y = val, label = lbl),
    color = BLUE, size = 2.8, fontface = "bold",
    nudge_y = 3, box.padding = 0.4, max.overlaps = 20, seed = 42
  ) +
  scale_x_date(date_breaks = "5 year", date_labels = "%Y") +
  scale_y_continuous(labels = \(x) paste0(x, "%")) +
  labs(
    title = paste0(LABEL, ": Double-Normalized (Log + Supertrend)"),
    subtitle = paste0("Deviation minus ", cycle_years, "yr supertrend | Dotted = +/-", Z_BAND, " SD"),
    x = NULL, y = "Normalized deviation (%)",
    caption = "Source: FRED DCOILWTICO. Cyclical over/undervaluation after removing secular trend."
  ) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 30, hjust = 1), legend.position = "none")

ggsave(paste0(PREFIX, "_3_normalized.png"), p3,
       width = width / 72 / sf, height = height / 72 / sf, dpi = 72 * sf)

cat("\nSaved:", paste0(PREFIX, c("_1_price.png", "_2_deviation.png", "_3_normalized.png")), "\n")