# Packages
library(data.table)
library(dplyr)
library(ggplot2)
library(ggdensity)
library(knitr)
library(tidyr)
library(MASS)
library(plotly)
library(patchwork)
library(lubridate)BUSN_33320_HW3
# Read in the Data
data_PJM <- fread('/Users/neilstein/Documents/Academic/Spring 25/Financing the Grid/HW3/Data/PJM Data.csv')
data_Chicago <- fread('/Users/neilstein/Documents/Academic/Spring 25/Financing the Grid/HW3/Data/Data Chicago.csv')
data_Milwaukee <- fread('/Users/neilstein/Documents/Academic/Spring 25/Financing the Grid/HW3/Data/Data Milwaukee.csv')
data_CorpusChristi<- fread('/Users/neilstein/Documents/Academic/Spring 25/Financing the Grid/HW3/Data/Data Corpus Christi.csv')
data_Austin <- fread('/Users/neilstein/Documents/Academic/Spring 25/Financing the Grid/HW3/Data/Data Austin.csv')VI. How Weather Interacts with Energy Pricing Continued
Question 25
The code below plots the joint density of average daily temperature and the logged standard deviation of the daily volatility of the RT-DA LMP spread. Is there an obvious pattern that emerges? Try identifying notable outliers in the lmp_temp_Chicago dataset. Conjecture as to why outliers emerged on those days.
lmp_temp_Chicago <- data_Chicago[,.(mean_lmp_rt = mean(total_lmp_rt), mean_lmp_da = mean(total_lmp_da), vol_spread = sd(spread), temp = mean(avg_temperature_degrees_fahrenheit), humidity = mean(relative_humidity_percent)), by = .(date, month, year)]
ggplot(lmp_temp_Chicago, aes(x=temp, y=log(vol_spread)) ) +
geom_hdr(xlim = c(-25, 100), ylim = c(0,10)) +
geom_point(shape = 21) +
labs(fill = '', alpha = 'Probability', x = 'Average Daily Temperature (Degrees
Fahrenheit)', y = 'Logged Standard Deviation of (Lagged DA LMP - RT LMP)') +
theme_classic() +
theme(legend.position = 'top') +
facet_grid(~year)
# Lets grab some of the most volatile days
thr <- quantile(lmp_temp_Chicago$vol_spread, .995)
outliers <- lmp_temp_Chicago[vol_spread > thr, .(date, temp, vol_spread)]
print(outliers) date temp vol_spread
<IDat> <num> <num>
1: 2021-03-22 57 100.3618
2: 2022-01-13 35 110.0645
3: 2022-06-13 73 458.2674
4: 2022-12-23 -4 883.7946
5: 2022-12-24 7 749.1124
6: 2023-01-10 42 105.5087Admittedly it is hard to conjecture over the entire span of the days, but after digging, it turns out that 2022-12-23 was one of the coldest December days in Chicago in the last 40 years. This, combined with the fact that it was peak holiday season, likely put excessive strain on the grind. The next biggest day is 2022-12-24 which is the following day, essentially extending the same assumption.
Furthermore, on June 13th, 2022, “an exceptionally high-topped and powerful supercell impacted the Chicago metropolitan area, with a height of 60,000–70,000 ft (18–21 km) as measured by multiple NEXRAD sites,” which would have likely strained the grid and disrupted services (e.g. trees toppling over and affecting power lines). This was the case, as NWS notes that “hundreds of trees were downed or damaged, which also resulted in widespread power outages.”
Citation: US Department of Commerce, NOAA. “June 13, 2022: Supercell Storm Brings a Swath of Severe Wind Damage and Two Tornadoes to the Chicago Metro.” National Weather Service, August 6, 2022. https://www.weather.gov/lot/2022jun13.
Question 26
The code below renders the plots above using a 3D density plot. Give it a try! What do th epeaks and valleys represent?
dens_Chicago <- kde2d(lmp_temp_Chicago[date>'2021-01-01']$temp,log(lmp_temp_Chicago[date>'2021-01-01']$vol_spread))
plot_ly(x = dens_Chicago$x, y=dens_Chicago$y, z = dens_Chicago$z) %>% add_surface() %>% layout(scene = list(xaxis = list(title = 'Avg. Temperature'), yaxis = list(title = 'Log.SD (DA LMP - RT LMP)'), zaxis = list(title = 'PDF')))Peaks on the plot would essentially be the mode of the joint distribution - basically, the pairs of Log.SD and Avg. Temp that occur most frequently. The valleys in comparison would be be pairs that never appear in our distribution.
We also consistently see higher uncertainty when averaging around the 20 degrees Fahrenheit, with a plateau between two peaks in the log of the SD between Day-Ahead and Real-Time pricing around log.sd = 3 and log.sd = 6, which falls well below the freezing point for water. As we determined in the previous assignment, there are other predictors for volatility in electricity pricing besides temperature, driven by supply-side and demand-side factors.
While temperature is an underlying factor in driving demand for energy (e.g. HVAC systems supplying hot/cold air), we aren’t seeing a distinct linear relationship between pricing and temperature evidenced in this plot. Perhaps with the two distinct peaks, we are seeing results of confounding factors that typically emerge with cold-snap temperatures, such as extreme weather events.
VII. Comparing Interconnected Transmission Regions
The states west of PJM belong to the Midcontinent Independent System Operator (MISO). Milwaukee and Chicago, two cities merely 90 minutes apart by car, belong to two different ISOs (Milwaukee in MISO and Chicago in PJM). In this section, we explore differences in LMPs and their responsiveness to weather patterns, in two separate (yet interconnected) parts of the grid.
Question 27
The code below plots the spread between average daily DA LMPs in Chicago and Oak Creek Power Plant (near Milwaukee). Consider the outlier days. Could they be explained by variations in weather? Plot the same figure for the spread in average daily temperatures to test your hypothesis.
lmp_cities <- data.table( date = data_Milwaukee$date, Milwaukee = data_Milwaukee$DA_LMP, Chicago = data_Chicago$total_lmp_da,
da_spread = data_Chicago$total_lmp_da - data_Milwaukee$DA_LMP )
ggplot(lmp_cities, aes(x=date,y=da_spread)) +
stat_summary(size=0.2) +
labs(x='',y='Chicago DA LMP - Milwaukee DA LMP ($/MWh)')+
theme_classic()Warning: Removed 48 rows containing non-finite outside the scale range
(`stat_summary()`).No summary function supplied, defaulting to `mean_se()`
lmp_cities <- data.table( date = data_Milwaukee$date, Milwaukee = data_Milwaukee$avg_temp_f, Chicago = data_Chicago$avg_temperature_degrees_fahrenheit,
da_spread = data_Chicago$avg_temperature_degrees_fahrenheit - data_Milwaukee$avg_temp_f )
ggplot(lmp_cities, aes(x=date,y=da_spread)) +
stat_summary(size=0.2) +
labs(x='',y='Chicago DA LMP - Milwaukee DA LMP ($/MWh)')+
theme_classic()No summary function supplied, defaulting to `mean_se()`
library(data.table)
# 1) Combine into one table and compute daily means
daily_spread <- data.table(
date = data_Milwaukee$date,
DA_LMP = data_Chicago$total_lmp_da,
MKE_LMP = data_Milwaukee$DA_LMP,
CHI_temp = data_Chicago$avg_temperature_degrees_fahrenheit,
MKE_temp = data_Milwaukee$avg_temp_f
)[, .(
da_spread = mean(DA_LMP - MKE_LMP),
temp_spread = mean(CHI_temp - MKE_temp)
), by = date]
# 2) Select the 10 dates with largest absolute DA‑LMP spread
top10 <- daily_spread[order(-abs(da_spread))][1:10]
# 3) Display
print(top10) date da_spread temp_spread
<IDat> <num> <num>
1: 2021-02-17 -128.33501 0.000000
2: 2022-12-25 117.29370 -2.166667
3: 2021-02-16 -111.00049 -1.000000
4: 2022-12-24 66.62397 -1.833333
5: 2022-06-21 -56.12767 -1.750000
6: 2021-02-18 -53.23197 1.000000
7: 2022-06-16 -52.41635 2.375000
8: 2021-02-15 -51.93056 4.000000
9: 2022-12-27 46.96838 -1.166667
10: 2022-06-17 -45.96440 3.125000It does not look like temperature difference was the cause in the substantial difference in pricing. To avoid looking at repeat entries for the same day since its hourly data, we took the spread mean for the date. Here we see that on the most extreme date, the temp spread was 0 yet the DA cost was very different.
Question 28
One might expect that weather in Milwaukee and Chicago to be similar. The code below generates scatter plots comparing average daily temperatures in Milwaukee and Chicago,showing temperatures in the two cities are closely aligned. Try creating a similar plot for average daily DA LMPs in the two cities. How closely aligned are the average daily DA LMPs in the two cities? Why might LMPs in the two cities differ more than temperatures?
temperature_cities <- data.table( date = data_Milwaukee$date, Milwaukee = data_Milwaukee$avg_temp_f,
Chicago = data_Chicago$avg_temperature_degrees_fahrenheit) %>% unique()
ggplot(temperature_cities, aes(x = Milwaukee, y = Chicago)) +
geom_point(alpha = 0.2) +
geom_abline(intercept = 0, slope = 1,color='red') +
labs(x='Milwaukee Average Daily Temperature (Degrees Fahrenheit)', y='Chicago Average Daily Temperature (Degrees Fahrenheit)', color = '') +
facet_wrap(~year(date)) +
theme_classic()
da_cities <- data.table( date = data_Milwaukee$date, Milwaukee = data_Milwaukee$DA_LMP,
Chicago = data_Chicago$total_lmp_da) %>% unique()
ggplot(da_cities, aes(x = Milwaukee, y = Chicago)) +
geom_point(alpha = 0.2) +
geom_abline(intercept = 0, slope = 1,color='red') +
labs(x='Milwaukee Average Daily DA', y='Chicago Average Daily DA', color = '') +
facet_wrap(~year(date)) +
theme_classic()Warning: Removed 48 rows containing missing values or values outside the scale range
(`geom_point()`).
This highlights what we saw in the spread of the DA and temperature which indicated that temperature was not the primary cause of DA differences. We see the two cities have clear relationship in terms of temperature which makes sense given proximity. We do not however find a linear relationship in the average daily DA cost between the two cities. That means the DA price in Chicago gives almost no information regarding the DA price in Milwaukee, especially when considering that they are covered by separate RTOs. There are many distinct factors, including the makeup of their urban economies, that would also play a role in driving energy demand as well.
Question 29
We can examine whether daily temperatures in Chicago and Milwaukee are drawn from a normal distribution using a Q-Q plot. Based on the Q-Q plot created from the code below, to what extent are these temperature series normal? Try creating a similar Q-Q plot comparing Austin and Corpus Christi temperatures to the normal distribution. As an optional step, try creating a Q-Q plot of Chicago average temperatures against Milwaukee temperatures. Are these series drawn from the same distribution?
ggplot(temperature_cities %>% pivot_longer(!date, names_to = 'series', values_to = 'temp'), aes(sample = temp, color = series)) +
geom_qq() +
geom_qq_line() +
labs(x='Theoretical Normal Distribution Quantiles', y='Sample Quantiles (Degrees Fahrenheit)', color = '') +
theme_classic() +
theme(legend.position = 'top')
Both Milwaukee and Chicago temperatures align closely with a normal distribution in the center but exhibit excess kurtosis in the tails. In our distribution we would expect fewer extremely hot days then a normal distribution but more extremely cold days then a normal distribution. Another thing to again note is the similarity between the two cities.
temperature_cities2 <- data.table( date = data_Austin$date, Austin = data_Milwaukee$avg_temp_f,
Corpus = data_CorpusChristi$avg_temp_f) %>% unique()
ggplot(temperature_cities2 %>% pivot_longer(!date, names_to = 'series', values_to = 'temp'), aes(sample = temp, color = series)) +
geom_qq() +
geom_qq_line() +
labs(x='Theoretical Normal Distribution Quantiles', y='Sample Quantiles (Degrees Fahrenheit)', color = '') +
theme_classic() +
theme(legend.position = 'top')
Again both cities align closely with a normal distribution in the center but are quite different in the tails. First we notice that in Corpus both extremely hot and cold days occur less frequently than a normal distribution would suggest. While Austin in comparison experiences more cold days then expected but less hot days than expected from a normal distribution. Interestingly the cities are only 200 miles away but their distributions are quite different.
If we take a look at the “Climate of Texas” document, we can see on page 146 a visual depiction of weather streams and air masses that run through Texas. Notably, Austin is along the pathway a very narrow polar jet stream from the Pacific and guided by the Rocky Mountains, experiences West winds coming through from Baja California and other areas Western areas, and from the Gulf of Mexico, which is “the predominant geographical feature affecting the state’s climate, moderating seasonal temperatures along the Gulf Coast and more importantly, providing the major source of precipitation for most of the state.” Corpus Christi lies directly along the Gulf, so it misses the narrow polar jet stream that Austin experiences and has a more intense relationship with the Gulf air.
Reference: https://www.twdb.texas.gov/publications/state_water_plan/2012/04.pdf
Question 30
In February 2021, Texas suffered a devastating winter storm that exacerbated the vulnerabilities of its power grid. The code below displays average daily temperatures and real-time LMPs in Chicago, Milwaukee, and Corpus Christi during this month. Note the scale for each plot. What do you observe from these plots? Conjecture as to the relationship between prices and weather. Do you think the relationship is linear? Describe the role and nuances of anomalous weather.
# 1) Define your date window
start <- as.Date("2021-02-01")
end <- as.Date("2021-02-28")
# 2) Compute daily summaries for each city
# — Chicago
chi <- data_Chicago[date >= start & date <= end,
.(avg_temp = mean(avg_temperature_degrees_fahrenheit),
avg_rt_lmp = mean(total_lmp_rt)),
by = date
][, city := "Chicago"]
# — Milwaukee
mke <- data_Milwaukee[date >= start & date <= end,
.(avg_temp = mean(avg_temp_f),
avg_rt_lmp = mean(RT_LMP)),
by = date
][, city := "Milwaukee"]
# — Corpus Christi (same structure as Milwaukee)
cc <- data_CorpusChristi[date >= start & date <= end,
.(avg_temp = mean(avg_temp_f),
avg_rt_lmp = mean(RT_LMP)),
by = date
][, city := "Corpus Christi"]
# 3) Combine into one long table
daily <- rbindlist(list(chi, mke, cc))
# 4) Top row: average temperature
p1 <- ggplot(daily, aes(x = date, y = avg_temp)) +
geom_line() +
facet_wrap(~city, nrow = 1) +
labs(x = NULL, y = "Average Temperature (°F)") +
theme_classic()
# 5) Bottom row: real‐time LMP
p2 <- ggplot(daily, aes(x = date, y = avg_rt_lmp)) +
geom_line() +
facet_wrap(~city, nrow = 1) +
labs(x = NULL, y = "Real‐Time LMP ($/MWh)") +
theme_classic()
# 6) Stack them
(p1 / p2) + plot_layout(heights = c(1, 1))
Question 30 - Interpretation of Our Results
The plots generated display a curious relationship with weather, rather than a straightforward ‘as X increases Y increases’ style of linear relationship, the plots here indicate a pronounced relationship between spikes in LMP and the the extreme ends of the weather relative to seasonal trends. Unseasonably warm and unseasonably cold temperatures cause discernible increases in the price, so the relationship is not strictly linear, however it appears correlated to extreme fluctuations. Crises in the power grid for Texas, in particular, seem to cause more pronounced swings in price likely due to the independent nature of their grid and its inability to import power from other, stable regions when crises occur, leading to massive spikes in price thousands of times greater than the baseline price (see the Corpus Christi Real Time LMP chart above).
Anomalous weather, in particular, draws excess power from the grid as customers engage their home climate controls simultaneously to maintain their homeostasis, additionally, in case of this particular 2021 winter storm in Texas damaged power grid infrastructure and lead to extremely hazardous road conditions, exacerbating the power needs from hospitals and other emergency services. Another quirk of anomalous weather such as this winter storm is that it disrupts the capacity of many energy generation sources, such as blocking sunlight for solar, wetting coal & reducing the energy generated when fired, and rendering wind turbines more volatile. Curiously, nuclear power actually functions somewhat better in cold weather such as this, so there is a strong case to be made for diversifying power sources; however Texas specifically is a generally warmer climate, where nuclear is somewhat less effective, so nuclear would not be a catch all solution to combat weather-related power concerns. Diversification of power sources and improved climate modeling would be the best tools to prepare for anomalous weather.
VIII. Contracts for Differences
Question 31
The code below plots day-over-day changes in the Chicago hourly DA LMP. What does each point on the plot represent? Why are the observations on for most days approximately symmetric around zero? Which days exhibited the largest day-over-day changes in DA LMP? How do day-to-day changes in DA LMP relate to contracts for differences (CfDs)?
ggplot(data_Chicago[order(date)] %>% .[,prev_lmp_da := shift(total_lmp_da, type = 'lag'), by = hour], aes(x = date, y = total_lmp_da - prev_lmp_da)) +
geom_point(shape = 21, size = 0.8, alpha = 0.5) +
labs(x = '',y = 'DA LMP (t+1) - DA LMP (t)') +
theme_classic()Warning: Removed 24 rows containing missing values or values outside the scale range
(`geom_point()`).
Each point represents the hour (h) day‑ahead price change from one day to the next: \[ \Delta \mathrm{DA\_LMP}_{t,h} \;=\; \mathrm{DA\_LMP}_{t,h} \;-\; \mathrm{DA\_LMP}_{t-1,h} \]
In an efficient forward market, upward forecast revisions (which push DA LMP up) and downward revisions occur with nearly equal frequency and magnitude. This assumption also tracks with the fact that Chicago is a well-developed market with sufficient alternative power sources and grid interconnections to mitigate extreme circumstances, such a chart in a smaller city that is more isolated (such as Honolulu, Hawaii) would see more aggressive fluctuation given the relative lack of grid interconnections to draw power and fewer standby power generation options. Hence most hourly changes cluster around zero, producing an approximately balanced “plus/minus” scatter.
# Compute the hourly deltas (day‑over‑day) by hour‑of‑day
setorder(data_Chicago, hour, date)
data_Chicago[, prev_lmp_da := shift(total_lmp_da), by = hour]
deltas <- data_Chicago[!is.na(prev_lmp_da),
.(date, hour, delta = total_lmp_da - prev_lmp_da)
]
# 2) Find the 3 largest single‑hour increases
top3_up <- deltas[order(-delta)][1:3, .(date, hour, biggest_up = delta)]
# 3) Find the 3 largest single‑hour decreases
top3_down <- deltas[order(delta)][1:3, .(date, hour, biggest_down = delta)]
# 4) Display
cat("Top 3 Hourly Increases:\n")Top 3 Hourly Increases:print(top3_up) date hour biggest_up
<IDat> <int> <num>
1: 2022-12-25 21 222.3772
2: 2022-12-25 18 214.0008
3: 2022-12-25 20 213.7910cat("\nTop 3 Hourly Decreases:\n")
Top 3 Hourly Decreases:print(top3_down) date hour biggest_down
<IDat> <int> <num>
1: 2023-07-29 16 -224.7552
2: 2023-07-29 17 -222.0931
3: 2023-07-29 15 -191.1384A Contract for Difference in power markets is a derivative whereby two parties agree to exchange, for each settlement period (hour/day), the difference between the reference price (usually the DA LMP) and a fixed strike price. In other words, your payoff in hour h and day t is:
\[ \Pi_{t,h} \;=\; \mathrm{Payoff}_{t,h} \;=\; \mathrm{DA\_LMP}_{t,h} \;-\; \mathrm{K}_{t,h} \]
Its mark‑to‑market change from day (t-1) to (t) is \[ \Delta \Pi_{t,h} \;=\; \Pi_{t,h} \;-\; \Pi_{t-1,h} \;=\; \bigl(\mathrm{DA\_LMP}_{t,h}-K\bigr) \;-\; \bigl(\mathrm{DA\_LMP}_{t-1,h}-K\bigr) \;=\; \Delta\mathrm{DA\_LMP}_{t,h}. \]
Thus, the hourly day‑ahead price changes you plotted are exactly the P&L drivers for a CfD position.
Question 32
Form your own hypotheses about how day-over-day changes in DA LMP may vary according to spreads between different dimensions of the Chicago data (e.g., day-of-week, temperature changes, etc.).Test them. Discuss how patterns in day-over-day changes in DA LMP may relate to energy storage needs.
Two main hypotheses come to mind when teasing out how DA LMP pricing data may vary in the Chicago data:
The day of the week: power consumption likely increases on weekends as people spend more time at home and their places of business do not fully shut down and still draw power. We predict that Thursday and Friday will be more volatile than other days of the week as people transition towards the weekend, but the lack of coordination as to when the weekend officially begins causes spikes in DA LMP.
Seasonal temperature instability: Summer and Winter seasons are generally predictable in terms of their power demands for cooling/heating respectively, however we theorize that the seasonal transitions in March & September are more volatile and therefore cause a level of unpredictability in the DA LMP, which may be borne out in the data.
# Hypothesis #1 -- Thurs. and Fri. are more volatile than the rest of the week
ggplot(data_Chicago %>%
arrange(date) %>%
mutate(prev_lmp_da = lag(total_lmp_da, default = NA), .by = hour) %>%
mutate(day_group = ifelse(day_of_week %in% c("Thursday", "Friday"), "Thursday & Friday", "Other Days"))
,aes(x = date, y = total_lmp_da - prev_lmp_da, color = day_group)) +
geom_point(shape = 21, size = 0.8, alpha = 0.5) +
labs(x = '', y = 'DA LMP (t+1) - DA LMP (t)', color = 'Day Grouping') +
scale_color_manual(values = c("Other Days" = "blue", "Thursday & Friday" = "green")) +
facet_wrap(~ day_group, ncol = 1) +
theme_classic()Warning: Removed 24 rows containing missing values or values outside the scale range
(`geom_point()`).
# Hypothesis #2 -- Sept. and March will be more volatile than than the rest of the year
ggplot(data_Chicago %>%
arrange(date) %>%
mutate(prev_lmp_da = lag(total_lmp_da, default = NA), .by = hour) %>%
# month data doesn't exist already in the data so it must be created
mutate(month = month(datetime_hb, label = TRUE)) %>%
mutate(month_group = ifelse(month %in% c("Mar", "Sep"), "March & September", "Other Months")),
aes(x = date, y = total_lmp_da - prev_lmp_da, color = month_group)) +
geom_point(shape = 21, size = 0.8, alpha = 0.5) +
labs(x = '', y = 'DA LMP (t+1) - DA LMP (t)', color = 'Seasonal Transition Months') +
scale_color_manual(values = c("Other Months" = "purple", "March & September" = "orange")) +
facet_wrap(~ month_group, ncol = 1) +
theme_classic()Warning: Removed 24 rows containing missing values or values outside the scale range
(`geom_point()`).
Looking through the results from the charting exercise above we don’t see a particularly strong response from testing hypothesis #2 (transition month weather impacting DA LMP), however for hypothesis #1 (Thursday/Friday vs. the rest of the week) we can see a meaningful split. For example, halfway through the 2023-2024 period there is a period of high volatility in DA LMP, and in our chart for hypothesis #1 you can see a massive spike upwards in price on Thursday & Friday and an equally significant DA LMP decrease on the other days of the week.
Earlier we discussed one possible explanation for this in the initial hypothesis, that being the lack of market coordination for when the weekend starts and power demand shifts among customers, which may highlight an opportunity for the market. Since we are unlikely to overcome the social/market coordination challenge of officially establishing when the weekend begins, improved energy storage from Tuesdays & Wednesdays prior to this consistent issue, where the DA LMP is negative, would help flatten the price towards zero on these days of difficult to predict demand.
Question 33
The code below creates a data set containing the ratio of day-over-day DA LMP spreads between Chicago and Milwaukee. That is, it computes the day-over-day change in hourly DALMP for Chicago and Milwaukee and then computes their ratio in the spread_ratio field. Use this data set to plot spread_ratio and interpret the magnitude of spread_ratio. What do larger values mean and how do they relate to incentives for energy transmission and storage given your interpretation, what has changed in 2023 relative to 2021 and 2022?
internodal_spread <- data.table(date = data_Chicago$date, hour = data_Chicago$hour,
day_of_week = data_Chicago$day_of_week, month = data_Chicago$month,
da_Chicago = data_Chicago$total_lmp_da, da_Milwaukee = data_Milwaukee$DA_LMP,
da_CorpusChristi = data_CorpusChristi$DA_LMP) %>% .[order(date)] %>% .[,':=' (prev_da_Chicago = shift(da_Chicago, fill = NA, type = 'lag'), prev_da_Milwaukee = shift(da_Milwaukee, fill = NA, type = 'lag'), prev_da_CorpusChristi = shift(da_CorpusChristi, fill = NA, type = 'lag')), by = hour] %>% .[, ':=' (Chicago_da_spread = da_Chicago - prev_da_Chicago, Milwaukee_da_spread = da_Milwaukee - prev_da_Milwaukee, CorpusChristi_da_spread = da_CorpusChristi - prev_da_CorpusChristi)] %>% .[, spread_ratio := Chicago_da_spread/Milwaukee_da_spread]# 1. Prepare a POSIX datetime so ggplot can draw a proper time series
internodal_spread[, datetime := as.POSIXct(paste(date, sprintf("%02d:00:00", hour)),
format = "%Y-%m-%d %H:%M:%S", tz = "UTC")]
ggplot(internodal_spread, aes(x = datetime, y = spread_ratio)) +
geom_line(alpha = 0.6) +
labs(
title = "Day over Day DA LMP Spread Ratio: Chicago vs Milwaukee",
x = "Date",
y = expression(frac(Delta~LMP[CHI], Delta~LMP[MKE]))
) +
theme_minimal()Warning: Removed 24 rows containing missing values or values outside the scale range
(`geom_line()`).
internodal_spread[, year := as.integer(format(date, "%Y"))]
summary_by_year <- internodal_spread[!is.na(spread_ratio), .(
median_ratio = median(spread_ratio, na.rm = TRUE),
p95_ratio = quantile(spread_ratio, 0.95, na.rm = TRUE)
), by = year]
print(summary_by_year) year median_ratio p95_ratio
<int> <num> <num>
1: 2021 0.010199560 14.61675
2: 2022 -0.027304295 29.62446
3: 2023 0.004694067 16.00362Question 33 - Interpretation
Our spread ratio charted above can be represented in simpler terms by the equation below:
\[ \text{spread_ratio} = \frac{\text{Chicago_da_spread } (\Delta LMP_{\text{CHI}})}{\text{Milwaukee_da_spread } (\Delta LMP_{\text{MKE}})} = \frac{LMP_{\text{CHI}}(t+1) - LMP_{\text{CHI}}(t)}{LMP_{\text{MKE}}(t+1) - LMP_{\text{MKE}}(t)} \]
Practically, what we are observing here is that the spread ratio details the comparative change in Chicago vs. neighboring Milwaukee. A value greater than 1, in this context specifically, means that Chicago has much more volatile pricing than Milwaukee, a value above 0 but below 1 means that both city’s respective volatility roughly aligned, and a negative value would mean that Chicago’s volatility moved in the opposite direction from Milwaukee’s. The larger the divergence (positive or negative) means that there prices swinging in comparatively different directions, meaning that we have big opportunities to introduce energy storage & arbitrage to settle the price differences.
From our output table, we can see that 2021 was a weakly positive year (median ratio = 0.01), so the prices moved in the same direction at the same time, so there was limited incentive for introducing energy transmission/storage measures since the prices approximately moved together. However, in 2022 we observe a weakly negative median ration (-0.02) as well as a much larger 95th percentile ratio (up to 29.6 from 14.6 in 2021), indicating that there is an increasing deviation between the two markets that introduces incentives for energy storage, arbitrage, and transmission. Further analysis would need to be conducted on the exact magnitude of the costs required to set up the energy storage and transmission systems, especially with 2023 introducing a median ratio of 0.004, meaning that the incentives are weakened and Chicago’s fluctuations in DA LMP are not responding much to the changes in Milwaukee’s DA LMP. There is still a weakly positive value, so there may be an opportunity here, but it is less likely to be significantly profitable given the 2013 numbers.
Question 34
In the previous question, we explored the ratio of day-over-day DA LMP spreads between Milwaukee and Chicago. The code below creates a similar dataset for the ratio of day-over-day DA LMP spreads between Chicago and Corpus Christi, TX; the spread_ratio_ercot variable holds this data (Corpus Christi belongs to the Electric Reliability Council of Texas,or “ERCOT”). What do larger values mean, and how do they relate to incentives for energy transmission and storage?
internodal_spread$spread_ratio_ercot <- internodal_spread$Chicago_da_spread / internodal_spread$CorpusChristi_da_spreadIn short, a larger spread_ratio_ercot flags hours when Chicago’s incremental value of energy (or your stored MWh) far exceeds what you’d get in ERCOT—driving both transmission bids from ERCOT→MISO and resulting in incentives for regional storage arbitrage strategies. It is worth noting that these ratios are based around simple pricing and do not include a discounting for the energy transfer loss from the point of origin to the point of destination, we are simply thinking through how the spread ratios could influence ERCOT’s incentives and decision making.
Question 35
The code below tells you the number of days that spread_ratio and spread_ratio_ercot were above 1, respectively, between 2021 and 2023. Interpret these results.
count1 <- sum(internodal_spread$spread_ratio > 1, na.rm = TRUE)
count2 <- sum(internodal_spread$spread_ratio_ercot > 1, na.rm = TRUE)
print(count1) [1] 8940print(count2) [1] 7361count3 <- sum(internodal_spread$spread_ratio > 0, na.rm = TRUE)
count4 <- sum(internodal_spread$spread_ratio < 0, na.rm = TRUE)
print(count3) [1] 13119print(count4)[1] 13088Interpretation:
34.1% (8940/26280) - Chicago’s day-over-day DA-LMP move exceeds Milwaukee’s. Which implies ~1/3rd of the time you would make money by shifting energy from Milwaukee to Chicago. This is a clear indication that in the sample time period (2021-2023) the improvements in transmission and storage could allow for premium generation for the individuals utilizing the technology. Additional research into the consistency of this trend and the cost of developing the expanded infrastructure to initiate these transfers would need to be explored further.
28.1% (7361/26280) - Chicago’s DA-LMP change exceeds Corpus Christi’s, which conversely means that 71.9% of the time swings in ERCOT exceed or are equal to that of Chicago. This implies that in most hours, moving energy from Chicago to ERCOT would capture more premium. As in part 34, this assumes feasible transfers with zero loss over the distance traveled, which we can essentially ignore with Milwaukee given how close it is to Chicago as compared to ERCOT (Texas).
50.1% (13119/26280) - the frequency at which the Chicago and Milwaukee nodes move in the same direction, while 49.9% (13088/26280) represents the frequency at which the Chicago and Milwaukee nodes move in the opposite direction. This essentially shows a 50/50 ratio, complicating the incentives for constructing energy transfer infrastructure. Further analysis would need to be done to see if this is common among neighboring cities in different ISOs, though more risk-tolerant business may decide to enter the arbitrage space based just on this coinflip-level of payoff odds.
Question 36
Construct a spread metric that considers the contemporaneous spread between Chicago and Milwaukee to those from the preceding day. Interpret this ratio. How does it differ from spread_ratio? Try comparing these two different series using a plot.
# 1) Compute the contemporaneous hourly spread CHI–MKE
internodal_spread[, cont_spread := da_Chicago - da_Milwaukee]
# 2) For each hour, grab the “same‐hour yesterday” spread
internodal_spread[, cont_spread_prev := shift(cont_spread, type = "lag"), by = hour]
# 3) Define the new ratio: today’s spread vs. yesterday’s spread
internodal_spread[, spread_ratio_contemp := cont_spread / cont_spread_prev]
# 4) Make a time series plot comparing the two ratios
# – spread_ratio = (ΔCHI)/(ΔMKE)
# – spread_ratio_contemp = (CHI−MKE)ₜ/(CHI−MKE)ₜ₋₂₄ₕ
melted <- melt(
internodal_spread,
id.vars = c("date","hour"),
measure.vars = c("spread_ratio","spread_ratio_contemp"),
variable.name = "metric",
value.name = "ratio"
)
ggplot(melted[!is.na(ratio)], aes(
x = as.POSIXct(paste(date, sprintf("%02d:00:00",hour)), tz="UTC"),
y = ratio, color = metric
)) +
geom_line(alpha = 0.6) +
facet_wrap(~metric, ncol = 1, scales = "free_y") +
labs(
x = "Date–Hour",
y = "Ratio",
title = "Comparing day‑over‑day‐spread‑ratio vs. contemporaneous‑spread‑ratio"
) +
theme_minimal()
The spread_ratio observe interday pricing, which is more of a macro approach of aggregating overall price ratio over a full day period, similar to how we observed the spread in temperature on a monthly basis and then narrowed our scope to a weekly average for a given month, resulting in a lower variance. In these cases, interday DA spread_ratio would have significant variability because of macro trends for expected supply and demand variables and possible externalities resulting from both.
By observing hourly for the contemporaneous spread, more often or not there is going to be consistency in how and why both cities price based on price and demand–we see spikes especially in cases of extreme weather events, but these likely last from a range of <1 hour to a 24-48 hour period and in target areas (i.e. not necessarily equally distributed between the two cities and moreso localized deviations). Hence why we see a positive contemporaneous ratio spike in 2022 for when Chicago experience high demand.
IX. Load and Price Duration Curves
Question 37
Use the code below to load a dataset, created by Bolarinwa Ajanaku,23, containing dailyobservations of PJM data from January 1, 2016 to December 31, 2019.24 This dataset includes a number of series, including load, RT LMP, as well as supply and prices for natural gas and wind.
daily_PJM <- fread('/Users/neilstein/Documents/Academic/Spring 25/Financing the Grid/HW3/Data/PJM Daily.csv')
colnames(daily_PJM) [1] "Date" "NG_Price" "Actual_load_GWh" "rtm_price"
[5] "Gas_GWh" "Gas_pct" "Wind_GWh" "Wind_pct"
[9] "Baseload_GWh" "Baseload_pct" "Rd_g_mln" "Rd_w_mln"
[13] "Rd_bl_mln" "Gas_UR" "Wind_UR" "Baseload_UR"
[17] "Gas_VF" "Wind_VF" "Baseload_VF" Question 38
The code below generates annual load duration curves for PJM, over the years 2016-2019. Based on the plot, which years featured higher demand for electricity? Try adapting the annual price duration curves for PJM, using the rtm_price variable in place of Actual_load_GWh. Which year featured the highest LMP prices? Is this consistent with your inference based on the load duration curve?
ldc <- daily_PJM[, .(Date, Actual_load_GWh)] %>% .[, Date := as.Date(Date, format = "%m/%d/%Y")] %>% .[, rank := frank(-Actual_load_GWh) / .N, by = year(Date)]
ggplot(ldc, aes(rank, Actual_load_GWh, color = factor(year(Date)))) +
geom_line() +
scale_x_continuous(expand = expand_scale(mult=c(0,0.1)), limits = c(0,NA))+
scale_y_continuous(expand = expand_scale(mult=c(0,0.1))) +
labs(x = 'Percent of Days', y = 'Actual Load (GWh)', color = '', title = 'PJM Annual Load Duration Curves') +
theme_classic() +
theme(plot.title = element_text(hjust=0.5), legend.position = 'top', plot.caption = element_text(hjust = 0))Warning: `expand_scale()` was deprecated in ggplot2 3.3.0.
ℹ Please use `expansion()` instead.
rtm <- daily_PJM[, .(Date, rtm_price)] %>% .[, Date := as.Date(Date, format = "%m/%d/%Y")] %>% .[, rank := frank(-rtm_price) / .N, by = year(Date)]
ggplot(rtm, aes(rank, rtm_price , color = factor(year(Date)))) +
geom_line() +
scale_x_continuous(expand = expand_scale(mult=c(0,0.1)), limits = c(0,NA))+
scale_y_continuous(expand = expand_scale(mult=c(0,0.1))) +
labs(x = 'Percent of Days', y = 'RTM Price)', color = '', title = 'PJM RTM Price Curves') +
theme_classic() +
theme(plot.title = element_text(hjust=0.5), legend.position = 'top', plot.caption = element_text(hjust = 0))
2018 for the load duration curve appears to dominate other years in actual load as we trend higher in percent of days, though there is a brief observation where 2016 does dominate 2018.
2018 looks to be consistently higher in overall PJM RTM pricing over all percent of days and is fairly consistent with our inference on the load duration curve. 2016 doesn’t appear to dominate, though, so we see that 2016 is crowded out.
Question 39
Try plotting load duration curves for Austin, Milwaukee, and Corpus Christi, using the“load_mw” field. Which cities typically have higher loads?
Question 39 - Interpretation
If we are observing in terms of overall GWh load, Milwaukee clearly dominates where even at ~87.5% of days, we are seeing 15,000 GWh load consistently across the three years whereas Corpus Christi tops out in 2023 around 6,000 GWh and Austin under 1500 GWh in 2022 at the same percent of days. It is also worth noting that this trend is consistent across all of the observed years (2021 - 2023)