Introduction

The focus of this study is to explore the probability of extreme pricing events in the electricity spot price using a binomial Bayesian model, the data we will be using was collected in a previous study the focus of which was examining energy pricing dynamics in relation to weather conditions.

The previous study found that although there were often many contributing factors which result in an extreme electricity pricing event (an extreme pricing event is when the spot price exceeds $5,000 MW/h) the weather and particularly high temperature weather is a common factor.

For example a plot from the previous study which shows the mean RRP of NSW spot price where we have examined the report issued by the Australian Energy Regulator (AER) during the period of 2009 - 2019, the points marked out on this plot represent a report issued by the regulator.

This plot shows 20 points where the AER have issued a report because the spot price has exceeded reporting threshold of $5,000 MW/h, 13 out of 20 of these reports indicated that the high temperature was a contributing factor in the extreme pricing event.

However these extreme pricing events are rare, in this data set once all data sources have been aggregate up to a daily level the total number of days with extreme pricing events is 101 out of 19,325 which represents 0.523 percent of the entire data set, so these are quite rare events, for this reason the aim of this modeling exercise is for analytical purposes only.

The Data Set

These data sets were collected in a previous study in which an attempt to model energy pricing trends using a linear model was conducted.

  1. Bureau of Meteorology (BOM) - This data set contains daily readings of minimum/maximum temperature, rainfall in millimeters and solar exposure readings, these readings have been collected from over 700 weather stations across the country, however for this paper we are only interested in the readings from weather stations located within major cities, as this is where weather conditions are most likely to have a significant impact on the electricity network.
  2. Australian Energy Market Operator (AEMO) - This data set consists of half hourly Recommended Retail Price (RRP) and total demand figures for wholesale electricity per state.
  3. Australian Bureau of Statistics (ABS) - We are utilising a shape file made available by the ABS to filter down the BOM data to readings from weather stations located in major cities.

Although the AEMO data set provides half hourly figures for electricity pricing we are constrained by the granularity of the BOM data set which is provided as daily figures, additionally the AEMO pricing is per state therefore we need to group the BOM data by state, and aggregate the AEMO data to daily figures.

Additionally since weather conditions in the major cities are likely to have a greater impact on the electricity grid than elsewhere in the state, weather stations located outside of major metropolitan areas will be filtered out of the data set before performing aggregations.

Binomial Logistic Bayesian Model

Initially in our modeling we will start with the absolute_max_temp (the maximum temperature reading from a major city) variable as our predictor of the extreme_price_event variable, this will be either 1 if the spot price exceeded $5,000.00/MWh else it will be 0, we will then expand the model by adding the state, followed by the year and season.

The rationale for choosing these variables as predictors are firstly maximum temperature in previous work high temperature days have been shown to align with extreme pricing events in the AEMO data, secondly in analysing the data across states some states tend to suffer more frequent extreme pricing events, and thirdly year and season as these events tend to occur in the warmer seasons and in some states appear to be increasing over time.

Now before we start modeling let’s take a quick look at the data visually through a histogram density plot of the avarage_max_temp variable in our data set, in the below plot the blue bars represent 101 observations where there was an extreme pricing event, and the dark grey bars represent the 1,9224 observations where there was no extreme pricing events, from the above plot we can see that the extreme pricing events are predominantly concentrated around the higher temperature readings.

ggplot(combined_data, aes(x = avarage_max_temp, y = ..density.., fill = extreme_price_event == 1)) +
  geom_histogram(bins = 30) + 
  scale_fill_manual(values = c("gray30", "skyblue"))

To model the probability of an extreme price event using the BOM weather data we will use and Binomial Bayesian model implemented using the rstanarm package, there are a number of link functions which could be used in training a model however the two most commonly used link functions for a binomial model are logit and probit, logit is generally recommended by rstanarm documentation for a number of reasons such as performance and better numerical stability, therefore our link function for this model will be logit (Mc-stan.org. 2019).

Since we are using the rstanarm package we will start out with using the default settings for the prior distributions, out of the box rstanarm provides weakly informative priors which provides moderate regularisation and are designed to perform well in a wide variety of scenarios, and for the first iteration of our model these defaults will be used (Goodrich, J. 2019).

extreme_price_fit1 <- stan_glm(
  extreme_price_event ~ absolute_max_temp
  , data = combined_data
  , family = binomial(link = "logit")
  , cores = 4
  , seed = 12345
)

Since there appears to be only a weak relationship between our predictor variables we will try setting the priors in the next iteration of our model to a Student t distribution which is heavier tailed and better suited when the coefficients and closer to zero, which will likely suite this data set (Mc-stan.org. 2019).

We will now update our model and set the prior and prior_intercept to a Student t distribution with 7 degrees of freedom, location 0 and scale of 2.5 than compare the updated model with the first iteration with default priors.

t_prior <- student_t(df = 7, location = 0, scale = 2.5, autoscale = FALSE)

extreme_price_fit2 <- update(extreme_price_fit1
       , prior = t_prior
  , prior_intercept = t_prior
)
loo1 <- loo(extreme_price_fit1, save_psis = TRUE, cores = 4)
loo2 <- loo(extreme_price_fit2, save_psis = TRUE, cores = 4)
model_comparison <- loo_compare(loo1, loo2)
knitr::kable(model_comparison)
elpd_diff se_diff elpd_loo se_elpd_loo p_loo se_p_loo looic se_looic
extreme_price_fit2 0.0000000 0.0000000 -496.8988 46.94899 3.733907 0.5793626 993.7976 93.89797
extreme_price_fit1 -0.0144897 0.2555198 -496.9133 47.12607 3.781170 0.5927862 993.8266 94.25213

Based on the output of the comparison of these two models it appears the the second model with the Student t prior distribution is the superior model, the key values to take note of in the above table are the elpd_loo (expected log predictive density) and looic (leave one out cross validation information criterion), in the preferred model the elpd_loo value is higher and the looic value is lower indicating that this is the prevailing model.

Now that we have a model we can start to visualise the predicted probability of an extreme pricing event as a function of the maximum daily temperature by plotting the coefficient estimates together with the observed outcomes, in the below plot the estimated probability is represented by the grey line, and we can see in this plot that the probability of an extreme pricing event is around 25% to 30% when the temperatures reach the high forties.

# Predicted probability as a function of x
pr_extreme_price <- function(x, ests) plogis(ests[1] + ests[2] * x)

# A function to slightly jitter the binary data
jitt <- function(...) {
  geom_point(aes_string(...), position = position_jitter(height = 0.05, width = 0.1), 
             size = 2, shape = 21, stroke = 0.2)
}

plot_probability <- function(fit, plot_title){
  ggplot(combined_data, aes(x = absolute_max_temp, y = extreme_price_event, color = extreme_price_event)) + 
  scale_y_continuous(breaks = c(0, 0.5, 1)) +
  jitt(x="absolute_max_temp") + 
  stat_function(fun = pr_extreme_price, args = list(ests = coef(fit)), 
                size = 2, color = "gray35") +
  labs(
    title = plot_title
    ) +
  xlab("BOM Maximum Temperature") +
  ylab("AEMO Extreme Pricing Event")
}

event_prob_plot1 <- plot_probability(extreme_price_fit2, "Probability of Extreme Pricing Event (Initial Model)")
event_prob_plot1

To evolve this model further we will now perform another two updates on the model, firstly adding the state predictor variable, and secondly with the year and season predictor variables, in previous work we have seen that there are differences in the susceptibility to extreme pricing events between the states so adding this to the model should prove informative.

Additionally in some instances the number of extreme pricing events appears to be increasing over time so adding the year should provide useful information to the model, as well as the season since we see more extreme pricing events in the warmer times of the year.

We will then perform a comparison of the models to determine which combination of predictor variables is performing better for predicting the probability of extreme pricing events.

extreme_price_fit3 <- update(extreme_price_fit2, formula = extreme_price_event ~ absolute_max_temp + state) 

combined_data$year <- factor(combined_data$year)
extreme_price_fit4 <- update(extreme_price_fit2, formula = extreme_price_event ~ absolute_max_temp + state + year + season) 

loo3 <- loo(extreme_price_fit3, save_psis = TRUE, cores = 4)
loo4 <- loo(extreme_price_fit4, save_psis = TRUE, cores = 4)

model_comparison <- loo_compare(loo1, loo2, loo3, loo4)
knitr::kable(model_comparison)
elpd_diff se_diff elpd_loo se_elpd_loo p_loo se_p_loo looic se_looic
extreme_price_fit4 0.00000 0.00000 -452.2927 41.62686 18.921678 2.2280441 904.5853 83.25371
extreme_price_fit3 -41.41465 11.31417 -493.7073 46.03650 7.225645 0.8251606 987.4146 92.07301
extreme_price_fit2 -44.60613 13.34784 -496.8988 46.94899 3.733907 0.5793626 993.7976 93.89797
extreme_price_fit1 -44.62062 13.48349 -496.9133 47.12607 3.781170 0.5927862 993.8266 94.25213

The output of this model comparison is indicating that the fourth model extreme_price_fit4 with the largest number of predictor variable is the superior model, and this is despite the fact that the loo diagnostic applies a penalty for additional predictors to help prevent overfitting in this instance the additional parameters of year and season have proven to be informative enough to overcome loo penalty.

The initial diagnostics are looking pretty good especially considering this is a pretty simple model, so let’s take a closer look at this model by visualising the coefficients for out parameters.

plot(extreme_price_fit4)

An important diagnostic to check when developing a Bayesian model is the is the trace plot for the Markov Chain Monte Carlo (MCMC) to check that the chains are converging well, the mcmc_trace function is a good option for performing this check as it has additional functionality which will help us identify where there may be problems with divergence, we will pass the nut parameters to the np argument which will highlight any divergence within the chains by adding a rug plot underneath, in our case there is no divergence to plot which is a good sign (Mc-stan.org. 2019).

np <- nuts_params(extreme_price_fit4)
posterior <- as.array(extreme_price_fit4)

color_scheme_set("mix-brightblue-gray")
mcmc_trace(posterior, np = np)
## No divergences to plot.

The diagnostics for our preferred model all looking pretty good, lets now revisit the plot of the probability estimates along with the observations we looked at earlier, and compare the two models side by side.

event_prob_plot2 <- plot_probability(extreme_price_fit4, "Probability of Extreme Pricing Event (Final Model)")
grid.arrange(event_prob_plot1 + theme(legend.position="bottom"), event_prob_plot2 + theme(legend.position="bottom"), ncol = 2)

We can see from these plots that the final model is showing a much higher probability of an extreme pricing event when the temperatures reach the high forties with a probability of around 45% as opposed to 35% in the previous model. The probability for the lower ranges remains flat at 0% if we were able to include information on planned maintenance of various network resources such as generators, transmission lines, electrical substations and state interconnectors this would likely improve our ability to gauge the probability of these events in the lower temperature ranges.

Conclusion

In this study we have shown that it is possible to model the probability of extreme electricity spot pricing events fairly well using only a very simple set of predictor variables temperature, state, year and season, and this model performs quite well in determining probability in the higher temperature ranges, however these pricing events are not limited to high temperature days and there are latent variables which could help in determining a more accurate probability especially for days of more normal weather conditions.

These latent variables would include data on the schedule of outages for maintenance on various electricity grid resources such as generators, transmission lines, electrical substations and state interconnectors, this information would prove quite informative as in previous work outages have been shown to be a common contributing factor to extreme electricity pricing events.

Additionally the role of interconnectors and how they create dependencies between various states has not been accounted for in this model, and these dependencies can have flow on effects causing high spot price events in one state to affect the spot price in a neighbouring state or event flowing onto states that are not directly connected to each other, there are already research efforts which are modeling these effects using copula models (Ignatieva, K. et al).

References

  1. Commonwealth of Australia, Bureau of Meteorology 2019, Solar Exposure Information, viewed 30 Aug 2019, http://www.bom.gov.au/climate/austmaps/about-solar-maps.shtml.

  2. Commonwealth of Australia, Bureau of Meteorology 2019, Daily Maximum Temperature, viewed 20 Aug 2019, http://www.bom.gov.au/jsp/ncc/cdio/weatherData/av?p_nccObsCode=122&p_display_type=dailyDataFile&p_startYear=&p_c=&p_stn_num=008137.

  3. Commonwealth of Australia, Bureau of Meteorology 2019, Daily Minimum Temperature, viewed 20 Aug 2019, http://www.bom.gov.au/jsp/ncc/cdio/weatherData/av?p_nccObsCode=123&p_display_type=dailyDataFile&p_startYear=&p_c=&p_stn_num=008137.

  4. Commonwealth of Australia, Bureau of Meteorology 2019, Daily rainfall, viewed 20 Aug 2019, http://www.bom.gov.au/jsp/ncc/cdio/weatherData/av?p_nccObsCode=136&p_display_type=dailyDataFile&p_startYear=&p_c=&p_stn_num=008137.

  5. Commonwealth of Australia, Bureau of Meteorology 2019, Daily global solar exposure, viewed 20 Aug 2019, http://www.bom.gov.au/jsp/ncc/cdio/weatherData/av?p_nccObsCode=193&p_display_type=dailyDataFile&p_startYear=&p_c=&p_stn_num=008137.

  6. Ignatieva, K. and Trück, S. (2016). Modeling spot price dependence in Australian electricity markets with applications to risk management. Computers & Operations Research, 66, pp.415-433.

  7. Mc-stan.org. Visual MCMC diagnostics using the bayesplot package. viewed 2 Nov 2019, https://mc-stan.org/bayesplot/articles/visual-mcmc-diagnostics.html

  8. Goodrich, J. Prior Distributions for rstanarm Models. viewed 2 Nov. 2019, https://cran.r-project.org/web/packages/rstanarm/vignettes/priors.html

  9. Mc-stan.org. (2019). Estimating Generalized Linear Models for Binary and Binomial Data with rstanarm. viewed 27 Oct. 2019, http://mc-stan.org/rstanarm/articles/binomial.html