Fitting a regression with ARIMA errors

Today let's understand how to fit a regression whose leftover error isn't independent noise, using a real electricity demand series to see exactly why that matters.

Here's the running example. Melbourne's electricity grid has to match demand every half hour, and demand climbs on hot days because every air conditioner in the city runs harder. We'll use 90 days of the grid's actual daily electricity demand, in gigawatt-hours (GWh), alongside that day's peak temperature, in degrees Celsius, from November 1, 2013 to January 29, 2014.

Plot one against the other and the relationship is easy to see, one point per day.

The points climb fairly steadily from left to right, and the correlation shown above the chart, r = 0.66, says the same thing in one number: hotter days tend to bring higher demand. That's the relationship a regression is built to capture, and it's exactly where we'll start.

The plain regression: demand on peak temperature

The simplest way to capture that relationship is an ordinary regression: predict demand from peak temperature, and let a single straight line summarize the two.

Build the daily series from vic_elec's half-hourly readings, calling the day's peak temperature MaxTemp and its total demand Demand, the same column names you will see in every fit from here on, then split it into a 90-day training window and the 7 days right after it.

RInteractive R
# Turn vic_elec's half-hourly demand and temperature into one daily row (GWh and degrees C), then split it into a 90-day training window and a 7-day test window library(tsibble) library(tsibbledata) library(fable) library(fabletools) library(dplyr) daily <- tsibbledata::vic_elec |> index_by(Day = Date) |> summarise( Demand = sum(Demand) / 1000 / 2, # half-hourly MW averaged into a daily GWh total MaxTemp = max(Temperature) ) train <- daily |> filter(Day >= as.Date("2013-11-01"), Day <= as.Date("2014-01-29")) test <- daily |> filter(Day >= as.Date("2014-01-30"), Day <= as.Date("2014-02-05")) nrow(train) #> [1] 90 nrow(test) #> [1] 7

  

That's 90 days to fit the model on, and 7 held out to test it later. Now fit the regression itself, and read its report.

RInteractive R
# Fit an ordinary regression of demand on peak temperature, and report its coefficients fit_tslm <- train |> model(tslm = TSLM(Demand ~ MaxTemp)) report(fit_tslm) #> Series: Demand #> Model: TSLM #> #> Residuals: #> Min 1Q Median 3Q Max #> -31.4420 -9.2630 0.6055 9.0603 33.8528 #> #> Coefficients: #> Estimate Std. Error t value Pr(>|t|) #> (Intercept) 62.4632 5.6694 11.018 < 2e-16 *** #> MaxTemp 1.8334 0.2226 8.236 1.52e-12 *** #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 #> #> Residual standard error: 13.84 on 88 degrees of freedom #> Multiple R-squared: 0.4353, Adjusted R-squared: 0.4289 #> F-statistic: 67.84 on 1 and 88 DF, p-value: 1.522e-12

  

TSLM stands for time series linear model, fable's version of an ordinary regression that also understands a tsibble's time index. Read the coefficients the way you would for any regression: the intercept, 62.4632, is the fitted line's predicted demand at a peak temperature of 0 degrees C, and the slope, 1.8334, says every extra degree of peak temperature adds 1.8334 GWh of demand. Both are far from zero (p = 1.52e-12 for the slope), and the model explains 43.53% of the variation in daily demand (R-squared = 0.4353).

The Residuals block above summarizes the residual, the gap between what the model predicts for a day and what demand that day actually was, across all 90 training days: at the extremes, the model was 31.4420 GWh too high on one day and 33.8528 GWh too low on another.

That looks like a solid regression. But a regression's p-values and standard errors are only trustworthy if its assumptions hold, and one of those assumptions is easy to miss on a time series. That's next.