This Research report delivers a comprehensive quantitative analysis of the United States energy complex. Using live EIA API data, merit-order dispatch modeling, supply–demand equilibria, and R-based simulations, it characterises pricing dynamics across electricity, natural gas, crude oil, LNG, nuclear, renewables, battery storage, and emerging hydrogen markets.
The U.S. bulk power system is divided between competitive ISO/RTO markets (serving ~70% of load) and vertically integrated bilateral markets in the Southeast and parts of the West. Seven federally-recognised ISOs/RTOs operate day-ahead and real-time energy auctions under FERC oversight.
Southeast Bilateral Markets (TVA, Duke Energy,
Southern Company, Dominion) operate under state-regulated integrated
resource planning. Power is traded bilaterally at negotiated prices,
with no centrally cleared energy auction. These markets cover ~30% of
U.S. load and are under increasing FERC scrutiny to adopt market
mechanisms.
Western Interconnect (outside
CAISO) is served by a patchwork of investor-owned utilities,
cooperatives, and federal power authorities (BPA, WAPA). The Western
Energy Imbalance Market (WEIM), administered by CAISO, provides a
voluntary 5-minute real-time balancing market across the region.
LIVE EIA DATA
Solving the “Other” problem. EIA’s operational-data
route returns ~45 fueltypeid codes at mixed
granularity — top-level aggregates (COW all coal,
PET all petroleum, BIO all biomass) sit
alongside their own sub-components. Summing everything double-counts; a
hand-maintained exclude list is brittle and silently leaks volume into
Other. Instead we request an explicit
mutually-exclusive, collectively-exhaustive (MECE) set
of 13 codes. The integrity check below dumps every code EIA returns and
confirms the curated stack reconciles to the reported all-fuels total —
verified live at a 0.00% residual for 2023 (4,183 TWh),
so Other is a genuine 0.2%, with pumped storage and other gases
(blast-furnace/refinery gas) broken out rather than hidden.
# EIA API v2 — annual net generation by fuel, all sectors, national.
#
# WHY A CURATED CODE LIST: this route exposes ~45 fueltypeid codes at MIXED
# granularity — top-level aggregates (COW=all coal, PET=all petroleum,
# BIO=all biomass) sit alongside their own components (BIT/SUB/LIG, PEL/PC,
# WWW/WAS/LFG…). Summing "everything except a hand-maintained exclude list" is
# fragile: it double-counted biomass (~84 vs 47 TWh) and petroleum (~28 vs 16),
# while pumped-storage/other-gases leaked into "Other". Requesting an explicit,
# mutually-exclusive (MECE) set instead is deterministic and self-reconciling.
#
# Verified live (2023): these 13 codes sum to 4,183.3 TWh — EXACTLY EIA's
# reported all-fuels total (residual 0). "Other" (OTH) is then a true 0.2%.
fuel_map <- c(
NG = "Natural Gas", OOG = "Other Gases",
NUC = "Nuclear", COW = "Coal",
HYC = "Hydro", HPS = "Pumped Storage",
WND = "Wind", SUN = "Solar",
GEO = "Geothermal", PET = "Petroleum",
WWW = "Biomass", WAS = "Biomass",
OTH = "Other"
)
gen_raw <- if (api_ok) {
fetch_eia(
route = "electricity/electric-power-operational-data",
frequency = "annual",
data_cols = "generation",
facets = list(location = "US", sectorid = "99",
fueltypeid = names(fuel_map)), # request only MECE codes
start = "2013",
n = 5000
)
} else NULL
# ── Integrity diagnostic: reconcile curated stack against EIA all-fuels total ──
if (api_ok) diagnose_fuel_stack(fuel_map, year = 2023)##
## ── Fuel-stack integrity check (2023) ──
## fueltypeid codes returned by EIA:
## [1] "ALL" "ANT" "AOR" "BIO" "BIS" "BIT" "COL" "COW" "DFO" "DPV" "FOS" "GEO"
## [13] "HPS" "HYC" "LFG" "LIG" "MLG" "MSB" "NG" "NGO" "NUC" "OB2" "OBW" "OOG"
## [25] "ORW" "OTH" "PC" "PEL" "PET" "RC" "REN" "RFO" "SPV" "STH" "SUB" "SUN"
## [37] "TPV" "TSN" "WAS" "WND" "WNS" "WNT" "WOC" "WOO" "WWW"
##
## EIA all-fuels total : 4183.3 TWh
## Curated MECE stack : 4183.3 TWh
## Residual ('Other' risk): 0.0 TWh (0.00% of total)# Clean & reshape (fall back to representative data if API unavailable)
if (!is.null(gen_raw) && nrow(gen_raw) > 0) {
gen_df <- gen_raw %>%
mutate(
fuel = dplyr::recode(fueltypeid, !!!fuel_map, .default = "Other"),
year = as.integer(period),
gen_twh = as.numeric(generation) / 1e3 # thousand MWh → TWh (unit fix)
) %>%
filter(!is.na(gen_twh)) %>%
group_by(year, fuel) %>%
summarise(gen_twh = sum(gen_twh, na.rm = TRUE), .groups = "drop")
data_source <- "EIA API v2 (live) — curated 13-code MECE stack"
} else {
# Representative EIA 2013-2023 national generation data (TWh)
gen_df <- tibble::tribble(
~year, ~fuel, ~gen_twh,
2023, "Natural Gas", 1755,
2023, "Coal", 669,
2023, "Nuclear", 775,
2023, "Wind", 425,
2023, "Solar", 238,
2023, "Hydro", 258,
2023, "Petroleum", 19,
2023, "Biomass", 55,
2023, "Geothermal", 17,
2022, "Natural Gas", 1688,
2022, "Coal", 773,
2022, "Nuclear", 772,
2022, "Wind", 435,
2022, "Solar", 177,
2022, "Hydro", 247,
2022, "Petroleum", 22,
2022, "Biomass", 57,
2022, "Geothermal", 17,
2021, "Natural Gas", 1574,
2021, "Coal", 899,
2021, "Nuclear", 778,
2021, "Wind", 379,
2021, "Solar", 130,
2021, "Hydro", 249,
2021, "Petroleum", 20,
2021, "Biomass", 57,
2021, "Geothermal", 17,
2020, "Natural Gas", 1617,
2020, "Coal", 774,
2020, "Nuclear", 790,
2020, "Wind", 338,
2020, "Solar", 91,
2020, "Hydro", 280,
2020, "Petroleum", 17,
2020, "Biomass", 56,
2020, "Geothermal", 17,
2019, "Natural Gas", 1582,
2019, "Coal", 966,
2019, "Nuclear", 809,
2019, "Wind", 295,
2019, "Solar", 72,
2019, "Hydro", 274,
2019, "Petroleum", 22,
2019, "Biomass", 60,
2019, "Geothermal", 17,
2015, "Natural Gas", 1331,
2015, "Coal", 1356,
2015, "Nuclear", 797,
2015, "Wind", 191,
2015, "Solar", 26,
2015, "Hydro", 249,
2015, "Petroleum", 24,
2015, "Biomass", 64,
2015, "Geothermal", 17,
2013, "Natural Gas", 1124,
2013, "Coal", 1581,
2013, "Nuclear", 789,
2013, "Wind", 168,
2013, "Solar", 9,
2013, "Hydro", 268,
2013, "Petroleum", 23,
2013, "Biomass", 65,
2013, "Geothermal", 17
)
data_source <- "EIA (representative historical data)"
}REPRESENTATIVE DATA
\[CF_i = \frac{\text{Actual Generation}_i}{\text{Nameplate Capacity}_i \times 8760\text{ hrs}}\]
Interpretation. “Single-phase” vs “dual-phase” is read here as the number of thermodynamic cycles: a simple-cycle gas turbine runs one open Brayton cycle and vents ~600 °C exhaust; a combined-cycle unit adds a second Rankine (steam) bottoming cycle via a Heat-Recovery Steam Generator (HRSG), recovering that exhaust to lift efficiency from ~40% to ~60%. The emerging supercritical-CO₂ (sCO₂) machine replaces the working fluid with CO₂ held above its critical point (31 °C, 73.8 bar) in a closed Brayton loop.
Efficiency ↔︎ heat rate: \[\eta_{\text{LHV}} = \frac{3412\ \text{Btu/kWh}}{\text{Heat Rate}} \qquad \text{Fuel Cost}\ (\$/\text{MWh}) = \frac{\text{Heat Rate}}{1000}\times P_{\text{gas}}\ (\$/\text{MMBtu})\]
A CCGT at 6,000 Btu/kWh burning $3.00 gas costs ~$18/MWh in fuel; a simple-cycle peaker at 10,500 Btu/kWh costs ~$31/MWh — the efficiency gap is the merit-order gap.
So “most efficient” is a single-shaft H/J-class CCGT; “most flexible” is a multi-shaft block or an aeroderivative. Desks care about both: efficiency sets the spark-spread breakeven, flexibility sets who captures the evening ramp.
Why sCO₂ matters for markets. Supercritical CO₂ is ~10× denser than steam, so the turbomachinery is dramatically smaller (a utility-scale sCO₂ turbine fits on a workbench-scale footprint), enabling fast thermal ramps and cheap dry cooling — attractive where water is scarce. The Allam-Fetvedt cycle burns gas in pure oxygen with an sCO₂ working fluid, so the combustion product is a high-pressure CO₂ stream that is pipeline-ready for sequestration or EOR at near-zero energy penalty — potentially decoupling firm gas generation from carbon exposure. Status is pre-commercial: DOE’s 10 MWe STEP pilot (San Antonio) and NET Power’s La Porte demo are the bellwethers; materials/recuperator cost and high-pressure sealing are the gating risks. For a trading desk, commercial sCO₂ would compress the spark-spread advantage peakers enjoy during scarcity while giving gas a low-carbon compliance pathway.
In competitive wholesale markets, electricity prices are set by the marginal plant — the most expensive generator needed to satisfy demand at any given moment. This is the foundational pricing mechanism in all ISO/RTO energy auctions.
Market-Clearing Price:
\[P_t = MC_{\text{marginal plant}}\]
where marginal cost is:
\[MC_i = \underbrace{\text{Fuel Cost}_i}_{\text{Heat Rate} \times P_{\text{fuel}}} + \text{Variable O\&M}_i + \underbrace{\text{Emissions Cost}_i}_{E_i \times P_{CO_2}}\]
Why does the most expensive plant set the price? Under uniform-price auction rules, every accepted bid — including cheap wind and nuclear — receives the same clearing price as the marginal (most expensive) plant dispatched. This creates an economic rent for inframarginal units and is the mechanism that makes zero-marginal-cost renewables profitable under energy-only market designs.
A wholesale power market is a uniform-price reverse (procurement) auction. In an ordinary auction buyers bid a price up; in a reverse auction, sellers (generators) offer prices down and the buyer — the ISO, acting for load — selects the cheapest offers until demand is met. Every cleared unit is then paid the single marginal clearing price, not its own offer (uniform-price, not pay-as-bid). All seven U.S. ISOs run this as a two-settlement system:
The auction is a constrained optimisation. SCUC minimises total as-offered cost:
\[\min_{g_i,\,u_i}\ \sum_i \Big( C_i^{\text{energy}}(g_i) + C_i^{\text{start}}\,u_i + C_i^{\text{noload}}\,u_i \Big)\]
subject to energy balance, transmission limits, ramp and min-up/min-down constraints. The LMP is the dual variable (shadow price) of the nodal energy-balance constraint — which is exactly why it decomposes into energy + congestion + loss (see Section 15):
\[\text{LMP}_n = \lambda_{\text{energy}} + \sum_{l}\mu_l\,\text{PTDF}_{l,n} + \nu_n\]
The critical design split is how capacity gets paid:
Why the design split matters for risk. Energy-only markets (ERCOT, SPP) concentrate a generator’s annual margin into a handful of scarcity hours, so revenue is fat-tailed and hedging shifts to the far right tail of the price distribution modelled in Section 11. Capacity markets (PJM, ISO-NE) pre-pay availability, damping energy spikes but adding a separate capacity-auction price to forecast. The 2021 Winter Storm Uri event — ERCOT pinned at the then-$9,000 cap for days — is the canonical illustration of energy-only tail risk, and drove the cap down to $5,000.
Generators are sorted by marginal cost and committed in ascending order until cumulative capacity meets demand. The price-setting plant sits at the intersection of the supply stack and the vertical demand curve.
\[\text{Supply}(Q) = \sum_{i} \mathbf{1}\left[MC_i \leq P^*\right] \cdot C_i\]
where \(C_i\) is installed capacity of plant \(i\) and \(P^*\) is the clearing price.
The merit order:
\[MC_{\text{wind}} \approx MC_{\text{solar}} < MC_{\text{nuclear}} < MC_{\text{hydro}} < MC_{\text{coal}} < MC_{\text{gas CCGT}} < MC_{\text{gas peaker}}\]
The duck curve describes the characteristic shape of net load (total load minus solar generation) in solar-heavy grids. As utility-scale and rooftop solar generation peaks around midday, net load plummets — creating a surplus that must be absorbed by storage, exports, or curtailment. By late afternoon, solar drops and net load surges, requiring extreme evening ramp capability.
\[\text{Net Load}_t = \text{Total Load}_t - \text{Solar}_{t} - \text{Wind}_t\]
Evening Ramp Rate (CAISO, 2023): ~14,000 MW over 3 hours (~4,700 MW/hr)
REPRESENTATIVE DATA
LIVE EIA DATA
# EIA API v2: Henry Hub Natural Gas Spot Price (monthly)
hh_raw <- if (api_ok) {
fetch_eia(
route = "natural-gas/pri/sum",
frequency = "monthly",
data_cols = "value",
facets = list(duoarea = "NUS", series = "RNGWHHD"),
start = "2015-01",
n = 500
)
} else NULL
if (!is.null(hh_raw) && nrow(hh_raw) > 0) {
hh_df <- hh_raw %>%
mutate(date = as.Date(paste0(period, "-01")),
price = as.numeric(value)) %>%
filter(!is.na(price)) %>%
arrange(date)
hh_source <- "EIA API v2 (live)"
} else {
# Representative monthly Henry Hub data
set.seed(42)
hh_df <- tibble(
date = seq(as.Date("2015-01-01"), as.Date("2023-12-01"), by = "month"),
price = c(
# 2015: low prices
rep(2.9,3), rep(2.6,3), rep(2.7,3), rep(2.3,3),
# 2016: trough
rep(2.2,3), rep(2.6,3), rep(2.9,3), rep(3.2,3),
# 2017
rep(3.0,3), rep(2.8,3), rep(2.9,3), rep(2.8,3),
# 2018: rise
rep(2.8,3), rep(2.9,3), rep(3.1,3), rep(4.5,3),
# 2019
rep(3.0,3), rep(2.5,3), rep(2.4,3), rep(2.4,3),
# 2020: COVID crash
rep(1.9,3), rep(1.7,3), rep(2.1,3), rep(2.6,3),
# 2021: recovery
rep(2.6,3), rep(2.7,3), rep(4.0,3), rep(5.9,3),
# 2022: energy crisis
rep(4.5,3), rep(7.2,3), rep(8.8,3), rep(5.5,3),
# 2023: collapse
rep(3.2,3), rep(2.2,3), rep(2.5,3), rep(2.7,3)
)
) %>% slice(1:108)
hh_source <- "EIA (representative historical data)"
}LIVE EIA DATA
storage_raw <- if (api_ok) {
fetch_eia(
route = "natural-gas/stor/wkly",
frequency = "weekly",
data_cols = "value",
facets = list(duoarea = "NUS", series = "NW2_EPG0_SWO_NUS_BCF"),
start = "2018-01-01",
n = 1000
)
} else NULL
if (!is.null(storage_raw) && nrow(storage_raw) > 0) {
stor_df <- storage_raw %>%
mutate(date = as.Date(period),
storage = as.numeric(value),
week = lubridate::week(date),
year = lubridate::year(date)) %>%
filter(!is.na(storage))
stor_source <- "EIA API v2 (live)"
} else {
# Synthetic seasonal storage pattern (Bcf)
weeks <- 1:52
base_seasonal <- 1700 + 1300 * sin((weeks - 12) * pi / 26)
set.seed(99)
stor_df <- bind_rows(lapply(2018:2023, function(yr) {
tibble(year = yr,
week = weeks,
date = as.Date(paste0(yr,"-01-01")) + (weeks - 1) * 7,
storage = pmax(100, base_seasonal + rnorm(52, 0, 80) +
ifelse(yr == 2022, -300, 0)))
}))
stor_source <- "EIA (representative seasonal pattern)"
}
# Compute 5-year average and bounds
stor_5yr <- stor_df %>%
filter(year %in% 2018:2022) %>%
group_by(week) %>%
summarise(avg = mean(storage, na.rm = TRUE),
lo = min(storage, na.rm = TRUE),
hi = max(storage, na.rm = TRUE))
stor_2023 <- stor_df %>% filter(year == max(stor_df$year))\[P_{\text{gas}} = f(\underbrace{S}_{\text{storage surplus/deficit}},\ \underbrace{W}_{\text{HDD/CDD}},\ \underbrace{L}_{\text{LNG exports}},\ \underbrace{Q}_{\text{production}})\]
Henry Hub embodies a balance of domestic supply, seasonal weather-driven demand, and the increasingly important LNG export floor (~12 Bcf/d in 2023). Storage injections/withdrawals versus the 5-year average norm are the primary short-term price signal used by trading desks.
LNG Arbitrage Profit:
\[\text{LNG}_{\text{arb}} = P_{\text{Asia (JKM)}} - \left(P_{\text{US (HH)}} + C_{\text{liquefaction}} + C_{\text{shipping}} + C_{\text{regasification}}\right)\]
Typical cost components (2023):
| Component | $/MMBtu |
|---|---|
| Liquefaction tolling | ~2.50 |
| LNG shipping (Gulf → NE Asia) | ~1.50–2.50 |
| Regasification | ~0.30–0.50 |
| Total supply chain cost | ~4.30–5.50 |
At HH = $2.57 and JKM = $13.00: arb ≈ $13.00 − ($2.57 + $5.00) = $5.43/MMBtu ✔ Profitable
LIVE EIA DATA
# WTI Cushing spot price
wti_raw <- if (api_ok) {
fetch_eia(
route = "petroleum/pri/spt",
frequency = "monthly",
data_cols = "value",
facets = list(series = "RWTC"),
start = "2015-01",
n = 500
)
} else NULL
# Brent spot price
brent_raw <- if (api_ok) {
fetch_eia(
route = "petroleum/pri/spt",
frequency = "monthly",
data_cols = "value",
facets = list(series = "RBRTE"),
start = "2015-01",
n = 500
)
} else NULL
# Parse or use fallback
parse_oil <- function(raw) {
if (!is.null(raw) && nrow(raw) > 0) {
raw %>%
mutate(date = as.Date(paste0(period, "-01")),
price = as.numeric(value)) %>%
filter(!is.na(price)) %>%
arrange(date) %>%
select(date, price)
} else NULL
}
wti_df <- parse_oil(wti_raw)
brent_df <- parse_oil(brent_raw)
if (is.null(wti_df) || is.null(brent_df)) {
# fall back for BOTH series if either fetch fails — otherwise the inner_join
# below joins a data frame against NULL and aborts the knit
dates <- seq(as.Date("2015-01-01"), as.Date("2023-12-01"), by = "month")
# Representative WTI & Brent monthly closing prices
wti_vals <- c(47,44,47,52,59,60,55,45,44,46,48,37,
31,26,38,45,47,48,43,42,45,44,52,53,
54,55,52,48,47,46,50,48,55,57,60,58,
57,62,57,55,62,68,70,68,70,71,60,55,
61,55,58,60,55,71,73,68,70,78,78,77,
80,91,100,95,88,102,112,95,88,95,92,80,
74,78,76,75,80,73,71,69,72,68,78,77,
77,72,78,81,75,73,72,77,80,78,77,79,
75,76,78,77,78,75,72,74,82,83,77,74)
brent_vals <- wti_vals + runif(length(wti_vals), 2, 6)
wti_df <- tibble(date = dates, price = head(wti_vals, length(dates)))
brent_df <- tibble(date = dates, price = head(brent_vals, length(dates)))
oil_source <- "EIA (representative historical data)"
} else {
oil_source <- "EIA API v2 (live)"
}
spread_df <- inner_join(wti_df, brent_df, by = "date", suffix = c("_wti","_brent")) %>%
mutate(spread = price_brent - price_wti)\[\text{Yield}_{\text{gasoline}} = f(\underbrace{\text{API Gravity}}_{\text{lighter = more gasoline}},\ \underbrace{\text{Sulfur \%}}_{\text{lower = cheaper to process}})\]
The Strategic Petroleum Reserve (SPR) is the world’s largest emergency crude oil stockpile, maintained by the U.S. Department of Energy in four salt cavern sites along the Gulf Coast (Bryan Mound, Big Hill, West Hackberry, Bayou Choctaw). Authorised capacity: ~714 million barrels. The SPR holds sour crude (~65%) and sweet crude (~35%) to match U.S. refinery feedstock needs.
LIVE EIA DATA
spr_raw <- if (api_ok) {
fetch_eia(
route = "petroleum/stoc/wstk",
frequency = "weekly",
data_cols = "value",
facets = list(duoarea = "NUS", product = "EPC0", series = "WCSSTUS1"),
start = "2000-01-01",
n = 2000
)
} else NULL
if (!is.null(spr_raw) && nrow(spr_raw) > 0) {
spr_df <- spr_raw %>%
mutate(date = as.Date(period),
spr_mb = as.numeric(value)) %>%
filter(!is.na(spr_mb)) %>%
arrange(date)
spr_source <- "EIA API v2 (live)"
} else {
# Representative SPR levels: ~600mb peak (2010) → 350mb (2023)
spr_df <- tibble(
date = seq(as.Date("2000-01-01"), as.Date("2023-12-01"), by = "month"),
spr_mb = c(
seq(542, 700, length.out = 120), # 2000–2010 fill
seq(700, 695, length.out = 60), # 2010–2015 plateau
seq(695, 640, length.out = 60), # 2015–2020 drawdown
seq(640, 595, length.out = 24), # 2020–2021
seq(595, 372, length.out = 24) # 2022–2023 Biden drawdown
) + rnorm(288, 0, 5)
) %>% slice(1:288)
spr_source <- "EIA (representative historical data)"
}Electricity demand is modelled as a linear function of weather, economic activity, and structural trends:
\[D_t = \alpha + \beta_1 \cdot HDD_t + \beta_2 \cdot CDD_t + \beta_3 \cdot GDP_t + \beta_4 \cdot EV_t + \varepsilon_t\]
where: - \(HDD_t\) = heating degree days (cold-weather demand) - \(CDD_t\) = cooling degree days (hot-weather AC demand) - \(GDP_t\) = real GDP index (industrial/commercial load) - \(EV_t\) = EV charging load term (structural electrification)
FRED API EIA STEO
# ── Live macro drivers ───────────────────────────────────────────────────────
# GDP proxy: FRED INDPRO (monthly industrial production index)
# Weather: EIA STEO ZWHDPUS / ZWCDPUS national degree days (FRED has none)
fetch_macro_drivers <- function(start = "2019-01-01", n_months = 60) {
end <- format(Sys.Date(), "%Y-%m-%d")
# ── GDP proxy: FRED industrial production index (INDPRO, monthly, live) ──────
indpro <- fetch_fred_series("INDPRO", start, end)
if (!is.null(indpro) && nrow(indpro) >= 12) {
indpro <- indpro %>%
dplyr::mutate(month = lubridate::floor_date(date, "month")) %>%
dplyr::group_by(month) %>%
dplyr::summarise(GDP_idx = mean(value, na.rm = TRUE), .groups = "drop") %>%
dplyr::arrange(month) %>%
dplyr::mutate(GDP_idx = 100 * GDP_idx / GDP_idx[1]) # rebase to 100
} else indpro <- NULL
# ── Weather: EIA STEO national degree days ──────────────────────────────────
# FRED does NOT publish HDD/CDD (verified: series 400 / empty search). EIA's
# Short-Term Energy Outlook does: ZWHDPUS (heating) / ZWCDPUS (cooling), U.S.
# population-weighted, monthly. STEO also carries a forward forecast, so we
# trim to observed months (<= today) before fitting.
start_ym <- format(as.Date(start), "%Y-%m")
hdd <- fetch_steo_series("ZWHDPUS", start = start_ym)
cdd <- fetch_steo_series("ZWCDPUS", start = start_ym)
weather <- NULL
if (!is.null(hdd) && !is.null(cdd)) {
weather <- dplyr::full_join(
hdd %>% dplyr::transmute(month = lubridate::floor_date(date, "month"), HDD = value),
cdd %>% dplyr::transmute(month = lubridate::floor_date(date, "month"), CDD = value),
by = "month"
) %>%
dplyr::filter(month <= Sys.Date()) %>% # drop STEO forward forecast
dplyr::arrange(month)
}
if (!is.null(indpro) && !is.null(weather)) {
df <- indpro %>%
dplyr::inner_join(weather, by = "month") %>%
dplyr::arrange(month) %>%
dplyr::slice_tail(n = n_months)
if (nrow(df) >= 12)
return(list(df = df, source = "FRED INDPRO (GDP) + EIA STEO HDD/CDD (live)"))
}
if (!is.null(indpro)) {
df <- indpro %>%
dplyr::arrange(month) %>%
dplyr::slice_tail(n = n_months) %>%
dplyr::mutate(
month_n = dplyr::row_number(),
HDD = pmax(0, 500 * (cos((month_n - 3) * pi / 6) + 1) / 2),
CDD = pmax(0, 300 * sin((month_n - 3) * pi / 6))
)
return(list(df = df %>% dplyr::select(month, HDD, CDD, GDP_idx),
source = "FRED INDPRO (GDP) + seasonal HDD/CDD proxy (STEO unavailable)"))
}
NULL
}
# Live GDP needs FRED; live weather needs the EIA key — run if either is set
macro_live <- if (fred_ok || api_ok) fetch_macro_drivers() else NULL
if (!is.null(macro_live)) {
demand_df <- macro_live$df %>%
dplyr::transmute(
month = month,
month_n = dplyr::row_number(),
HDD = HDD, CDD = CDD, GDP_idx = GDP_idx
)
macro_source <- macro_live$source
} else {
set.seed(42)
n <- 60
demand_df <- tibble(
month = seq(as.Date("2019-01-01"), by = "month", length.out = n),
month_n = 1:n,
CDD = pmax(0, 300 * sin((month_n - 3) * pi / 6) + rnorm(n, 0, 40)),
HDD = pmax(0, 500 * (cos((month_n - 3) * pi / 6) + 1) / 2 + rnorm(n, 0, 60)),
GDP_idx = 100 + cumsum(rnorm(n, 0.3, 0.5))
)
macro_source <- if (fred_ok) "FRED (partial — fallback synthetic weather/GDP)"
else "Simulated (set fred_key param for live FRED data)"
if (!fred_ok)
message("Note: Set params$fred_key for live GDP/HDD/CDD from FRED.")
}
# True demand model coefficients (representative regression estimates)
alpha <- 3500; b1 <- 0.45; b2 <- 0.60; b3 <- 8.5
if (!"demand_mw" %in% names(demand_df)) {
demand_df <- demand_df %>%
mutate(
demand_mw = alpha + b1 * HDD + b2 * CDD + b3 * (GDP_idx - 100) + rnorm(n(), 0, 80)
)
}
# Fit OLS
model <- lm(demand_mw ~ HDD + CDD + GDP_idx, data = demand_df)
demand_df$fitted <- predict(model)
cat("Macro data source:", macro_source, "\n")## Macro data source: FRED INDPRO (GDP) + EIA STEO HDD/CDD (live)cat("Model Coefficients:\n")## Model Coefficients:print(round(coef(model), 3))## (Intercept) HDD CDD GDP_idx
## 1788.011 0.433 0.642 17.521cat("\nR-squared:", round(summary(model)$r.squared, 4))##
## R-squared: 0.4706EIA RTO + ISO APIs
# Live ISO fuel mix — EIA RTO grid monitor (hourly, aggregated over a recent window)
iso_fuel_live <- fetch_iso_fuel_mix(days = 3)
if (!is.null(iso_fuel_live)) {
iso_fuel_summary <- iso_fuel_live %>% dplyr::arrange(ISO, dplyr::desc(gen_gwh))
cat("Live ISO fuel mix — trailing ~3 days (GWh), EIA RTO grid monitor:\n")
print(iso_fuel_summary %>% dplyr::group_by(ISO) %>%
dplyr::slice_head(n = 3) %>% dplyr::ungroup())
iso_data_source <- "EIA API v2 — electricity/rto/fuel-type-data (hourly)"
} else {
iso_fuel_summary <- NULL
iso_data_source <- "Representative (EIA RTO API unavailable)"
}## Live ISO fuel mix — trailing ~3 days (GWh), EIA RTO grid monitor:
## # A tibble: 21 × 3
## ISO fuel gen_gwh
## <chr> <chr> <dbl>
## 1 CAISO Solar 445.
## 2 CAISO Natural Gas 321.
## 3 CAISO Wind 145.
## 4 ERCOT Natural Gas 1672.
## 5 ERCOT Solar 676.
## 6 ERCOT Wind 633.
## 7 ISO-NE Natural Gas 474.
## 8 ISO-NE Nuclear 199.
## 9 ISO-NE Other 52.0
## 10 MISO Natural Gas 2081.
## # ℹ 11 more rows# Live nodal LMP is pulled directly from CAISO OASIS in Section 15.3.
# fetch_iso_lmp_snapshot() also offers a gridstatus/Python path for manual use,
# but it is NOT called here to avoid triggering a py_install during knit.
cat("\nNodal LMP: see Section 15.3 for a live CAISO OASIS decomposition.\n")##
## Nodal LMP: see Section 15.3 for a live CAISO OASIS decomposition.This section builds the merit-order dispatch in idiomatic R, adds a deterministic price–demand sensitivity curve, and then runs a multi-factor Monte Carlo that jointly shocks demand, fuel prices, and generator availability — yielding a realistic clearing-price distribution with tail-risk metrics (P95, CVaR, scarcity probability) rather than a handful of discrete price points.
# ── Plant Stack ──────────────────────────────────────────────────────────────
plants_sim <- tibble::tribble(
~plant, ~capacity_mw, ~mc_per_mwh,
"Wind", 1000, 5,
"Solar", 800, 5,
"Nuclear", 1200, 15,
"Coal", 1500, 30,
"Gas CCGT", 2000, 45,
"Gas Peaker", 500, 120
) %>% arrange(mc_per_mwh) %>%
mutate(cum_cap = cumsum(capacity_mw))
# ── Single dispatch ───────────────────────────────────────────────────────────
dispatch_price <- function(demand_mw, stack = plants_sim) {
marginal <- stack %>% filter(cum_cap >= demand_mw) %>% slice(1)
if (nrow(marginal) == 0) return(list(price = NA, plant = "Insufficient Capacity"))
list(price = marginal$mc_per_mwh, plant = marginal$plant)
}
# Base case
d0 <- 3500
result <- dispatch_price(d0)
cat("Base demand:", d0, "MW\n")## Base demand: 3500 MWcat("Market clearing price: $", result$price, "/MWh\n")## Market clearing price: $ 30 /MWhcat("Marginal plant:", result$plant, "\n\n")## Marginal plant: Coal# ── Sensitivity across demand range ──────────────────────────────────────────
sensitivity <- tibble(
demand = seq(500, 7000, by = 100)
) %>%
mutate(price = sapply(demand, function(d) {
p <- dispatch_price(d)
if (is.na(p$price)) 200 else p$price
}),
marginal_plant = sapply(demand, function(d) dispatch_price(d)$plant))
# ── Multi-factor Monte Carlo ─────────────────────────────────────────────────
# The old model varied demand ONLY, through a fixed stack, so the clearing price
# could take just 6 discrete values (the six plant MCs) — a near-degenerate
# histogram. A credible distribution needs the joint, correlated drivers a desk
# actually hedges:
# (1) demand shocks (2) fuel-price shocks that RE-PRICE the gas fleet
# (3) forced generator outages that pull capacity out of the merit order.
# A common latent factor z couples cold/hot snaps to BOTH high load and gas spikes.
mc_stack <- tibble::tribble(
~plant, ~capacity_mw, ~fuel, ~heat_rate, ~vom, ~for_rate,
"Wind", 1000, "wind", NA_real_, 1, 0.05,
"Solar", 800, "solar", NA_real_, 1, 0.03,
"Nuclear", 1200, "uranium", 10.4, 4, 0.02,
"Coal", 1500, "coal", 10.2, 5, 0.08,
"Gas CCGT", 2000, "gas", 7.0, 4, 0.05,
"Gas Peaker", 500, "gas", 11.0, 12, 0.10
)
fuel_base <- c(gas = 3.00, coal = 2.50, uranium = 0.70, wind = 0, solar = 0) # $/MMBtu
VOLL <- 200 # scarcity price cap ($/MWh) when demand exceeds available capacity
# One stochastic dispatch: joint demand + gas-price + outage draws re-price and
# re-sort the stack, then the marginal available unit sets the clearing price.
draw_clearing_price <- function(i, mean_d = 3500, sd_d = 400) {
z <- rnorm(1) # latent weather factor
demand <- mean_d + sd_d * z + rnorm(1, 0, 150) # load rises with z
gas_p <- fuel_base["gas"] * exp(0.28 * z + rnorm(1, 0, 0.18)) # gas spikes with z (lognormal)
avail <- runif(nrow(mc_stack)) > mc_stack$for_rate # forced-outage draw
s <- mc_stack[avail, , drop = FALSE]
fp <- fuel_base; fp["gas"] <- gas_p
mc <- ifelse(is.na(s$heat_rate), s$vom, s$heat_rate * fp[s$fuel] + s$vom)
ord <- order(mc); s <- s[ord, ]; mc <- mc[ord]
k <- which(cumsum(s$capacity_mw) >= max(100, demand))[1]
if (is.na(k)) VOLL else unname(mc[k])
}
set.seed(2024)
mc_n <- MC_DRAWS
if (USE_PAR && requireNamespace("furrr", quietly = TRUE) &&
requireNamespace("future", quietly = TRUE)) {
n_workers <- max(1L, parallel::detectCores() - 1L)
future::plan(future::multisession, workers = n_workers)
mc_prices <- furrr::future_map_dbl(seq_len(mc_n), draw_clearing_price,
.options = furrr::furrr_options(seed = TRUE))
future::plan(future::sequential)
mc_engine <- sprintf("furrr parallel (%d workers)", n_workers)
} else {
mc_prices <- purrr::map_dbl(seq_len(mc_n), draw_clearing_price)
mc_engine <- "purrr sequential"
}
# ── Risk metrics ─────────────────────────────────────────────────────────────
mc_q <- quantile(mc_prices, c(.05, .50, .95), names = FALSE)
cvar95 <- mean(mc_prices[mc_prices >= mc_q[3]]) # expected shortfall, worst 5%
mc_df <- tibble(price = mc_prices)
cat(sprintf("Monte Carlo engine : %s\n", mc_engine))## Monte Carlo engine : furrr parallel (7 workers)cat(sprintf("Draws : %s (demand x gas price x outages)\n", format(mc_n, big.mark = ",")))## Draws : 100,000 (demand x gas price x outages)cat(sprintf("Mean / median price : $%.2f / $%.2f per MWh\n", mean(mc_prices), mc_q[2]))## Mean / median price : $25.21 / $25.35 per MWhcat(sprintf("P5 - P95 range : $%.2f - $%.2f per MWh\n", mc_q[1], mc_q[3]))## P5 - P95 range : $11.28 - $30.53 per MWhcat(sprintf("CVaR(95%%) shortfall : $%.2f per MWh\n", cvar95))## CVaR(95%) shortfall : $50.60 per MWhcat(sprintf("P(price > $60 spike) : %.2f%%\n", 100 * mean(mc_prices > 60)))## P(price > $60 spike) : 0.53%cat(sprintf("P(scarcity = VOLL) : %.2f%%\n", 100 * mean(mc_prices >= VOLL)))## P(scarcity = VOLL) : 0.34%Battery storage profits by charging during low-price hours and discharging during high-price hours:
\[\text{Profit} = \sum_{t} \left(P_t^{\text{peak}} - P_t^{\text{off-peak}}\right) \times E \times \eta_{\text{RT}}\]
where \(\eta_{\text{RT}}\) is the round-trip efficiency and \(E\) is energy throughput (MWh).
Key Parameters (4-hr BESS, 2023 U.S. Average):
| Parameter | Value |
|---|---|
| Round-trip efficiency \(\eta_{RT}\) | 85–90% |
| Cycle degradation | ~20% over 3,000 cycles |
| Installed cost ($/kWh) | \(250–350/kWh | | LCOE (\)/MWh dispatched) |
| Target spread for profitability | >$40/MWh |
# Simulate a year of hourly prices with a peaky structure
set.seed(55)
hours <- 0:23
base_price <- 35
price_profile <- tibble(
hour = hours,
avg_price = base_price + 15 * sin((hours - 6) * pi / 12) +
20 * pmax(0, sin((hours - 16) * pi / 4)) +
rnorm(24, 0, 5)
)
# Battery strategy: charge cheapest 4 hrs, discharge most expensive 4 hrs
n_charge <- 4 # hours to charge
RTE <- 0.87
capacity_mw <- 100 # MW
energy_mwh <- 400 # MWh (4-hr battery)
charge_hrs <- price_profile %>% slice_min(avg_price, n = n_charge)
discharge_hrs <- price_profile %>% slice_max(avg_price, n = n_charge)
charge_cost <- sum(charge_hrs$avg_price) * capacity_mw / 1000 # $k
discharge_rev <- sum(discharge_hrs$avg_price) * capacity_mw * RTE / 1000
daily_profit <- discharge_rev - charge_cost
annual_profit <- daily_profit * 365
cat("Daily charge cost: $", round(charge_cost * 1000), "\n")## Daily charge cost: $ 9194cat("Daily discharge rev: $", round(discharge_rev * 1000), "\n")## Daily discharge rev: $ 21436cat("Daily gross profit: $", round(daily_profit * 1000), "\n")## Daily gross profit: $ 12242cat("Annual gross profit: $", round(annual_profit * 1000, -3), "\n")## Annual gross profit: $ 4468000\[\text{H}_2 = \frac{\text{Electricity Input (kWh)}}{\text{Electrolyzer Efficiency}} \quad \text{(Green)}\]
The U.S. Inflation Reduction Act (IRA, 2022) created the Clean Hydrogen PTC (§ 45V), offering up to $3/kg for green hydrogen with lifecycle emissions < 0.45 kg CO₂e/kg H₂. This subsidy effectively bridges the current grey/green cost gap, making green H₂ competitive at renewable electricity prices below ~$30/MWh.
Load following requires reactors to ramp output up and down to track demand, unlike traditional baseload operation. Key constraints: - Xenon poisoning: \(^{135}\)Xe buildup after power reduction temporarily limits ramping up - Fuel rod thermal cycling: fatigue reduces fuel lifetime - Economic incentive: negative power prices during high-renewable periods may justify curtailment
French EDF fleet successfully load-follows at 25–100% within 30 min. U.S. fleet operates at maximum output to maximise revenue given high fixed costs. New SMR designs (Natrium) incorporate molten salt thermal storage to decouple electricity output from reactor thermal power, enabling flexible grid services.
In an unconstrained market, the price \(P^*\) solves:
\[S(P^*) = D(P^*)\]
where: \[S(P) = \sum_i C_i \cdot \mathbf{1}[MC_i \leq P] \quad \text{(supply stack)}\] \[D(P) = D_0 - \epsilon \cdot P \quad \text{(price-elastic demand, } \epsilon \ll 1 \text{)}\]
# Supply function: cumulative capacity at each MC level
supply_pts <- plants_sim %>%
select(mc_per_mwh, cum_cap) %>%
mutate(cum_start = lag(cum_cap, default = 0)) %>%
rename(price = mc_per_mwh, q_supply = cum_cap)
# Demand function: D(P) = 5000 - 2*P (slightly elastic)
D0 <- 5000
d_elas <- 2
price_range <- seq(0, 150, by = 1)
demand_curve <- tibble(
price = price_range,
q_demand = pmax(0, D0 - d_elas * price_range)
)
# Find equilibrium numerically
eq_data <- inner_join(
tibble(price = price_range,
q_supply = sapply(price_range, function(p) {
plants_sim %>% filter(mc_per_mwh <= p) %>% pull(cum_cap) %>% max(., 0)
})),
demand_curve, by = "price"
)
# Equilibrium: where supply first meets or exceeds demand
eq_row <- eq_data %>%
filter(q_supply >= q_demand) %>%
slice(1)
cat("Equilibrium Price: $", eq_row$price, "/MWh\n")## Equilibrium Price: $ 45 /MWhcat("Equilibrium Quantity:", round(eq_row$q_demand), "MW\n")## Equilibrium Quantity: 4910 MWReal ISO/RTO markets do not clear as a simple unconstrained equilibrium. SCED adds:
\[\min_{g_i} \sum_i MC_i \cdot g_i\]
subject to: - \(\sum_i g_i = D_t\) (energy balance at each node) - \(0 \leq g_i \leq \bar{C}_i\) (capacity limits) - \(|F_l| \leq \bar{F}_l, \quad \forall l \in \mathcal{L}\) (transmission flow limits) - \(\Delta g_i / \Delta t \leq R_i\) (ramp rate limits)
This yields Locational Marginal Prices (LMPs) that vary by node, reflecting energy cost, congestion rent, and marginal loss:
\[\text{LMP}_n = \underbrace{\lambda_{\text{energy}}}_{\text{system balance}} + \underbrace{\sum_{l \in \mathcal{L}} \mu_l \cdot \text{PTDF}_{l,n}}_{\text{congestion rent}} + \underbrace{\nu_n}_{\text{marginal loss}}\]
Energy-only markets (ERCOT, SPP) clear on energy alone; capacity-structured markets (PJM, ISO-NE, NYISO, MISO, CAISO) add a reliability price adder.
# ── ISO-specific plant stacks with nodal assignments ─────────────────────────
build_iso_stack <- function(iso) {
switch(iso,
PJM = tibble::tribble(
~plant, ~node, ~capacity_mw, ~mc_per_mwh,
"Nuclear", "EAST", 3500, 15,
"Coal", "WEST", 2000, 32,
"Gas CCGT","EAST", 4500, 42,
"Wind", "WEST", 1200, 5,
"Solar", "EAST", 800, 5,
"Gas Peak","EAST", 1000, 110
),
ERCOT = tibble::tribble(
~plant, ~node, ~capacity_mw, ~mc_per_mwh,
"Wind", "WEST", 3500, 5,
"Solar", "SOUTH", 2000, 5,
"Gas CCGT","NORTH", 4000, 38,
"Coal", "NORTH", 800, 35,
"Gas Peak","SOUTH", 1500, 95
),
MISO = tibble::tribble(
~plant, ~node, ~capacity_mw, ~mc_per_mwh,
"Coal", "MISO_N", 3500, 30,
"Gas CCGT","MISO_N", 3000, 40,
"Wind", "MISO_N", 2500, 5,
"Nuclear", "MISO_S", 1500, 15,
"Solar", "MISO_S", 600, 5
),
CAISO = tibble::tribble(
~plant, ~node, ~capacity_mw, ~mc_per_mwh,
"Solar", "SP15", 4000, 5,
"Gas CCGT","NP15", 2500, 50,
"Imports", "ZP26", 1500, 55,
"Hydro", "NP15", 1200, 10,
"Gas Peak","SP15", 800, 120
),
NYISO = tibble::tribble(
~plant, ~node, ~capacity_mw, ~mc_per_mwh,
"Nuclear", "UPST", 2000, 15,
"Gas CCGT","NYC", 3000, 48,
"Hydro", "UPST", 1500, 8,
"Wind", "LONG", 800, 5,
"Gas Peak","NYC", 600, 115
),
`ISO-NE` = tibble::tribble(
~plant, ~node, ~capacity_mw, ~mc_per_mwh,
"Gas CCGT","MA", 2500, 52,
"Nuclear", "MA", 1500, 15,
"Wind", "ME", 900, 5,
"Hydro", "ME", 600, 8,
"Gas Peak","CT", 500, 125
),
SPP = tibble::tribble(
~plant, ~node, ~capacity_mw, ~mc_per_mwh,
"Wind", "SPP_N", 4000, 5,
"Gas CCGT","SPP_S", 2500, 35,
"Coal", "SPP_S", 2000, 28,
"Solar", "SPP_N", 800, 5,
"Gas Peak","SPP_S", 600, 90
)
)
}
iso_names <- c("PJM","ERCOT","MISO","CAISO","NYISO","ISO-NE","SPP")
iso_stacks <- setNames(lapply(iso_names, build_iso_stack), iso_names)
# Demand is set to ~78% of each (illustrative) stack's capacity so a mid/upper-
# merit unit is always marginal. The earlier hard-coded GW peak loads were ~10x
# these toy stacks, leaving no marginal unit and emptying the LMP table.
iso_demands <- vapply(iso_stacks, function(s) round(sum(s$capacity_mw) * 0.78), numeric(1))
sced_results <- purrr::imap(iso_demands, function(demand, iso) {
sced_lmp_iso(iso, demand, iso_stacks[[iso]], iso_network_params(iso))
})
# Compile nodal LMP table
lmp_table <- purrr::map_dfr(sced_results, function(r) {
if (is.null(r$nodes)) return(NULL)
r$nodes %>%
dplyr::mutate(
ISO = r$iso, lambda_energy = r$lambda_energy,
market_type = r$market_type, binding = r$binding_constraint
)
})
cat("SCED/LMP Decomposition — 7 ISOs/RTOs\n\n")## SCED/LMP Decomposition — 7 ISOs/RTOsif (nrow(lmp_table) > 0) {
lmp_table %>%
dplyr::select(ISO, node, lmp_energy, lmp_cong, lmp_loss, lmp_total,
capacity_adder, lmp_all_in, market_type) %>%
dplyr::arrange(ISO, dplyr::desc(lmp_all_in)) %>%
dplyr::mutate(dplyr::across(where(is.numeric), ~ round(.x, 2))) %>%
print(n = 30)
} else {
cat("No nodal LMPs produced (check ISO stacks vs demand).\n")
}## # A tibble: 18 × 9
## ISO node lmp_energy lmp_cong lmp_loss lmp_total capacity_adder lmp_all_in
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 CAISO SP15 55 2.7 0.3 58 18 76
## 2 CAISO NP15 55 0 0 55 18 73
## 3 CAISO ZP26 55 -1.2 0.13 53.9 18 71.9
## 4 ERCOT NORTH 38 2.4 0.19 40.6 0 40.6
## 5 ERCOT WEST 38 0.6 0.05 38.6 0 38.6
## 6 ERCOT SOUTH 38 -1.2 0.1 36.9 0 36.9
## 7 ISO-NE CT 52 1.08 0.15 53.2 20 73.2
## 8 ISO-NE MA 52 0 0 52 20 72
## 9 ISO-NE ME 52 -0.72 0.1 51.4 20 71.4
## 10 MISO MISO… 40 1.5 0.22 41.7 11 52.7
## 11 MISO MISO… 40 -0.72 0.11 39.4 11 50.4
## 12 NYISO NYC 48 3 0.4 51.4 22 73.4
## 13 NYISO LONG 48 1 0.13 49.1 22 71.1
## 14 NYISO UPST 48 0 0 48 22 70
## 15 PJM EAST 42 2.8 0.29 45.1 14 59.1
## 16 PJM WEST 42 -1.2 0.13 40.9 14 54.9
## 17 SPP SPP_N 35 0.75 0.11 35.8 0 35.8
## 18 SPP SPP_S 35 -0.5 0.07 34.6 0 34.6
## # ℹ 1 more variable: market_type <chr>cat("\nSystem energy component (λ) by ISO:\n")##
## System energy component (λ) by ISO:purrr::map_dfr(sced_results, ~ tibble(
ISO = .x$iso, lambda_energy = .x$lambda_energy,
marginal_unit = .x$marginal_unit, market_type = .x$market_type
)) %>% print()## # A tibble: 7 × 4
## ISO lambda_energy marginal_unit market_type
## <chr> <dbl> <chr> <chr>
## 1 PJM 42 Gas CCGT capacity
## 2 ERCOT 38 Gas CCGT energy_only
## 3 MISO 40 Gas CCGT capacity
## 4 CAISO 55 Imports capacity
## 5 NYISO 48 Gas CCGT capacity
## 6 ISO-NE 52 Gas CCGT capacity
## 7 SPP 35 Gas CCGT energy_onlyLIVE CAISO OASIS
The SCED model in 15.2 is illustrative. Here we pull actual day-ahead LMPs from CAISO OASIS (public, no key) and decompose each into its true components. CAISO’s LMP is distinctive: on top of energy + congestion + loss it carries a marginal GHG cost (MGHG) from California’s cap-and-trade program — a live carbon price embedded directly in the wholesale clearing price.
\[\text{LMP}_n^{\text{CAISO}} = \underbrace{\text{MCE}}_{\text{energy}} + \underbrace{\text{MCC}}_{\text{congestion}} + \underbrace{\text{MCL}}_{\text{loss}} + \underbrace{\text{MGHG}}_{\text{cap-and-trade carbon}}\]
# NP15 (north) and SP15 (south) trading hubs — the classic CAISO congestion pair.
caiso_hubs <- c(NP15 = "TH_NP15_GEN-APND", SP15 = "TH_SP15_GEN-APND")
caiso_lmp <- fetch_caiso_hubs(caiso_hubs, days_back = 7) # OASIS is public, no key
if (!is.null(caiso_lmp) && nrow(caiso_lmp) > 0) {
caiso_summary <- caiso_lmp %>%
dplyr::group_by(hub) %>%
dplyr::summarise(dplyr::across(c(energy, congestion, loss, ghg, lmp),
~ round(mean(.x, na.rm = TRUE), 2)), .groups = "drop")
cat("CAISO day-ahead LMP decomposition (24-hr mean, $/MWh):\n")
print(caiso_summary)
caiso_src <- "CAISO OASIS PRC_LMP (live day-ahead)"
} else {
caiso_summary <- NULL
caiso_src <- "CAISO OASIS unavailable (throttled/offline) — see modeled SCED in 15.2"
cat(caiso_src, "\n")
}## CAISO day-ahead LMP decomposition (24-hr mean, $/MWh):
## # A tibble: 2 × 6
## hub energy congestion loss ghg lmp
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 NP15 19.7 -0.17 -0.54 8.58 27.6
## 2 SP15 19.7 -1.6 -1.09 8.58 25.6Model vs reality. The illustrative SCED (15.2) can’t
be literally replaced by OASIS data — real PTDF/shift-factor matrices
aren’t recoverable from prices alone — but this live pull
validates it: the same energy + congestion + loss structure
appears in the observed CAISO decomposition, plus the
California-specific GHG adder. Congestion (MCC) and loss (MCL) are
typically small at the hubs on an unstressed day and swing hardest
during the evening ramp and heat events; the NP15–SP15 spread is the
tradeable north–south congestion signal. To extend this to other ISOs,
swap in ERCOT (ercot.com public API),
PJM (Data Miner 2), or the Python
gridstatus library via
fetch_iso_lmp_snapshot().
The spark spread measures the profitability of a gas-fired power plant:
\[\text{Spark Spread} = P_{\text{power}} - \left(\text{Heat Rate} \times P_{\text{gas}}\right)\]
A positive spark spread indicates gas generators are in-the-money; negative spreads force economic curtailment.
# Simulate spark spread using model price series
set.seed(33)
n_months <- 96 # 8 years
dates_ss <- seq(as.Date("2016-01-01"), by = "month", length.out = n_months)
gas_prices <- pmax(1.5, 3.0 + cumsum(rnorm(n_months, 0, 0.3)) +
1.5 * sin(1:n_months * pi / 6)) # seasonal
power_prices <- 25 + 7 * gas_prices + rnorm(n_months, 0, 8) +
5 * sin((1:n_months + 3) * pi / 6) # correlated + seasonal
heat_rate <- 7.0 # MMBtu / MWh (efficient CCGT)
spread_ss_df <- tibble(
date = dates_ss,
gas = gas_prices,
power = power_prices,
spark_spread = power_prices - heat_rate * gas_prices
)
# Dark spread (coal)
coal_prices <- 2.5 + rnorm(n_months, 0, 0.2)
heat_rate_coal <- 10.2
spread_ss_df <- spread_ss_df %>%
mutate(dark_spread = power - heat_rate_coal * coal_prices)
cat("Mean Spark Spread:", round(mean(spread_ss_df$spark_spread), 2), "$/MWh\n")## Mean Spark Spread: 25.65 $/MWhcat("% Months Negative Spark:", round(mean(spread_ss_df$spark_spread < 0) * 100, 1), "%\n")## % Months Negative Spark: 0 %Electricity Markets: Merit-order dispatch with rising renewable penetration is compressing average wholesale prices while simultaneously increasing intra-day volatility. The duck curve is intensifying in CAISO and spreading to ERCOT and PJM as solar capacity scales. Battery storage is the structural beneficiary, with 4-hour BESS increasingly co-located with solar to capture evening ramp spreads.
Natural Gas: The U.S. has transitioned from price taker to global price setter via LNG exports (~12 Bcf/d, 2023). Henry Hub is now loosely tethered to JKM and TTF via arb mechanics. Storage trajectory vs the 5-year norm remains the highest-frequency tradeable signal. Marcellus oversupply structurally caps HH below $4/MMBtu absent an LNG demand shock or severe winter.
Crude Oil: WTI–Brent spread dynamics are driven by U.S. export logistics capacity and Cushing inventory. API gravity premiums for light-sweet crude (Permian, Bakken) persist as U.S. refiners are configured for heavier feedstocks. SPR drawdown risk remains elevated with stocks at 40-year lows.
Emerging Markets: Green hydrogen economics hinge on IRA § 45V PTC + sub-$30/MWh renewable electricity. SMR deployment (NuScale, TerraPower) faces project finance and regulatory timeline risk but represents a long-dated optionality on 24/7 carbon-free firm capacity. Watch for advanced nuclear licensing activity at the NRC as a leading indicator.
## Report generated: 2026-07-13 17:18:10 UTC## R version: R version 4.4.0 (2024-04-24)## Data sources: EIA Live API + FRED Live API + Representative / Simulated## Monte Carlo draws: 100,000