Producing a Time Series model to accurately forecast electricity demand in Victoria

Time Series Analysis - Assignment 3

Setup

# Clear the environment and the console
rm(list = ls()); cat("\f")



#----------------------------#
# Install required libraries #
#----------------------------#
packages <- c("ggplot2", "ggpubr", "readxl", "readr", "dplyr", "tidyr", "psych", "stringr",
              "lubridate", "knitr", "outliers", "MVN", "TSA", "tseries", "lmtest", "FSAdata",
              "forecast", "matrixcalc", "car", "corpcor", "scales", "QuantPsyc",
              "urca", "rugarch", "fGarch", "tswge",
              "imputeTS", "shiny", "coda", "rjags", "runjags", "ks", "epiDisplay", "fastDummies")
# Install any packages not already installed
installed_packages <- packages %in% rownames(installed.packages())
if (any(installed_packages == FALSE)) { install.packages(packages[!installed_packages]) }
#-------------------------#
# Load libraries          #
#-------------------------#
invisible(lapply(packages, library, character.only = TRUE))
#-------------------------#
# Clear the environment and the console
rm(list = ls()); cat("\f")

Required packages

library(TSA)
library(readr)
library(stringr)
library(dplyr)
library(lubridate)

Introduction

Aims/Objectives

This report summarises the process followed to identify a suitable model for forecasting electricity demand in Victoria, Australia. The report covers analysis of the Victorian Electricity Demand series, and the process followed to identify and fit suitable models for the series, select the most appropriate model, perform diagnostic checking, and forecast future electricity demand using the selected model.

Methodology

Creating custom functions to reduce repetition

To reduce repetition of some R codes, custom functions for producing ACF and PACF plots, and for sorting AIC and BIC scores were created and loaded to the R environment.

#-------------------------#
# function for producing ACF and PACF plots
#-------------------------#
acf_pacf <- function(ts_object, plot_note){
 # Produce ACF and PACF plots for a time series object.
 # Optional: specify a plot note (as a string) to improve plot titles.
 # Example: acf_pacf(price_ts, "Box-Cox transformed, first differenced price series.")
 par(mar=c(5,6,4,1)+.1) #set plot window margins
 par(mfrow=c(1,2))
 if (missing(plot_note)) {
 acf(ts_object, main ="ACF plot", lag.max = min(c(length(ts_object), 50)))
 pacf(ts_object, main ="PACF plot", lag.max = min(c(length(ts_object), 50)))
 } else {
 acf(ts_object, main=paste("ACF plot of", plot_note), lag.max = min(c(length(ts_object),50)))
 pacf(ts_object, main=paste("PACF plot of", plot_note), lag.max = min(c(length(ts_object),50)))
 }
 par(mfrow=c(1,1))
}
#-------------------------#
# function for sorting AIC and BIC scores
#-------------------------#
sort.score <- function(x, score = c("bic", "aic")){
 if (score == "aic"){
 x[with(x, order(AIC)),]
 } else if (score == "bic") {
 x[with(x, order(BIC)),]
 } else {
 warning('score = "x" only accepts valid arguments ("aic","bic")')
 }
}

Data description

The Victorian Electricity Demand data contains 2,106 observations of daily electricity demand in Victoria, measured in Mega Watt hours (MWh) from January 1st, 2015, to October 6th, 2020. The data was compiled by Alex Kozlov and published to the Kaggle website in 2020 (Kozlov, 2020).

Retrieving the data

After importing the Victorian Electricity Demand data into the R environment, a data summary confirmed theimported data contained 2,106 observations of electricity demand, with corresponding dates from January 1st, 2015, to October 6,th 2020, which aligned with the data description provided on Kaggle (Kozlov, 2020).

ve <- read.csv("complete_dataset.csv", header = TRUE) %>% dplyr::select(c(date, demand))
summary(ve)

     date               demand      
 Length:2106        Min.   : 85094  
 Class :character   1st Qu.:109964  
 Mode  :character   Median :119586  
                    Mean   :120035  
                    3rd Qu.:130436  
                    Max.   :170654  

ve$date <- as_date(ve$date)
summary(ve)

      date                demand      
 Min.   :2015-01-01   Min.   : 85094  
 1st Qu.:2016-06-10   1st Qu.:109964  
 Median :2017-11-18   Median :119586  
 Mean   :2017-11-18   Mean   :120035  
 3rd Qu.:2019-04-28   3rd Qu.:130436  
 Max.   :2020-10-06   Max.   :170654  

Table 1: Result from a summary of the Daily Victorian Electricity Demand series.

The minimum and maximum values of electricity demand suggested daily electricity demand could vary widely over time. This was not a surprising finding when considering Victoria’s average maximum temperature in Winter is just 14 degrees Celsius, while average temperatures in Autumn and Spring are fairly moderate (Victorian Government, 2022), and cold temperatures could lead to high demand for electric heating and electricity, compared to the demand in the Autumn and Spring months.

Exploring the data

As the Victorian Electricity Demand data contained daily observations, the data was converted to a time series object with a frequency of 365 and plotted to observe the key time series elements.

ts.ve <- ts(ve$demand, start=c(2015, 1), frequency=365)
plot(ts.ve, type="o", xlab="Time", ylab="Electricity Demand (MWh)", 
 main = "Time series plot of the Victorian Daily Electricity Demand series.")

Figure 1: Time series plot daily electricity demand in Victoria, Australia (Jan 2015 to Oct 2020)

The time series plot of the Daily Electricity Demand in Victoria series showed the following time series traits;

1. Trend - A very slight downward trend

2. Seasonality – clear seasonality, with peaks in the middle of each year, coinciding with cold temperatures in Winter and high demand for electricity for heating.

3. Changing variance – Some evidence of changing variance was observed. For example, there was a clear difference in the amount of variance in January 2018 compared to July 2016.

4. Change point – There was no clear evidence of a change point or intervention observed in the series.

5. Behaviour – The behaviour of the series appeared to be moving average (MA) although the presence of seasonality in the series made it difficult to clearly identify the behaviour of the series.

The clear presence of seasonality in the series suggested a SARIMA model would be appropriate for modelling the Victorian Daily Electricity Demand series. However, SARIMA models cannot handle ultra-high frequency data, such as daily or hourly observations, so the raw Victorian Daily Electricity Demand data was aggregated to monthly averages and converted to a time series object that would be used for SARIMA modelling and forecasting electricity demand.

ve$month <- floor_date(ve$date, "month")
ve <- ve %>% group_by(month) %>% summarise(mean = mean(demand))
summary(ve)

     month                 mean       
 Min.   :2015-01-01   Min.   :104801  
 1st Qu.:2016-06-08   1st Qu.:113529  
 Median :2017-11-16   Median :119516  
 Mean   :2017-11-15   Mean   :119848  
 3rd Qu.:2019-04-23   3rd Qu.:126645  
 Max.   :2020-10-01   Max.   :137856  

Table 2: Summary of the Victorian Electricity Demand data, aggregated to monthly averages.

A summary of the monthly average electricity demand identified that the minimum value for electricity demand had risen from 85,904 in the daily data, to 104,801 MWh in the monthly averages. Meanwhile, the maximum had fallen from 170,654 in the daily data, to 137,856 MWh in the monthly averages. The reduced range for the monthly average electricity demand was not surprising considering a monthly average would “smooth-out” demand, meaning days with particularly high electricity demand (e.g., particularly cold days), would be aggregated with milder days that coincide with lower values for electricity demand.

An ACF plot of the Monthly Average Electricity demand data (Figure 2) identified recurring peaks at a frequency of 12, which made sense considering the data contained monthly averages. As such, a frequency of 12 was applied to the data when converting it to a time series object.

ts.ve <- ts(ve$mean)
acf(ts.ve, lag.max = 50, main="ACF plot of the Monthly Electricity Demand series.")

Figure 2: ACF plot of the monthly average electricity demand in Victoria data (Jan 2015 to Oct 2020).

ts.ve <- ts(ve$mean, start=c(2015, 1), frequency=12)
plot(ts.ve, type="o", xlab="Time", ylab="Monthly Average Electricity Demand (MWh)", 
 main = "Time series plot of the Monthly Electricity Demand series.")

Figure 3: Time series plot monthly average electricity demand in Victoria data (Jan 2015 to Oct 2020).

The time series plot of the Monthly Average Electricity Demand in Victoria series (Figure 3) showed the following time series traits;

1. Trend - A very downward trend that was clearer than the daily demand series

2. Seasonality – there was clear evidence of seasonality in the monthly average demand series, and Figure 4 supported this observation further. June, July, and August consistently had the highest average demand, which was not surprising given these are winter months, and demand for electricity for heating would be particularly high during these times.

3. Changing variance – Some evidence of changing variance was observed, similar to the daily demand series.

4. Change point – There was no clear evidence of a change point or intervention observed in the series.

5. Behaviour – The series appeared to exhibit a combination of autoregressive (AR) and moving average (MA) behaviour, although the presence of seasonality made it difficult to clearly identify the behaviour of the series.

plot(ts.ve, type="l", xlab="Time", ylab="Monthly Average Electricity Demand (MWh)", 
 main = "Time series plot of the Monthly Electricity Demand series.")
points(y=ts.ve, x=time(ts.ve), pch=as.vector(season(ts.ve)))

Figure 4: Time series plot monthly average electricity demand in Victoria with point markers for months (Jan 2015 to Oct 2020).

The points in Figure 3 also showed evidence of succeeding measurements being related to one another. This relationship was highlighted in the scatter plot of neighbouring pairs in Figure 5, where a moderate-strength, positive correlation was identified.

par(mfrow=c(1,1))
plot(y=ts.ve, x=zlag(ts.ve), 
 ylab="Monthly Average Electricity Demand (MWh)", 
 xlab="Average Electricity Demand (MWh) in previous month",
 main = "Scatter plot of Average Monthly Electricity Demand (MWh) in consecutive Months")

Figure 5: Scatter plot of each monthly electricity demand observation and the measurement taken in the previous month.

Lagging the series confirmed the strength of the relationship, with a correlation between monthly demand of 0.617, indicating that a month with low demand tended to be followed by another month with low demand, a month with moderate demand with another month of moderate demand, and a month with high demand with another month of high demand.

y <- ts.ve # Assign the demand data to y
x <- zlag(ts.ve) # Generate first lag of the demand series
index <- 2:length(x) # Create an index to get rid of the first NA value in x
cor(y[index],x[index]) # Calculate correlation between numerical values in x and y

[1] 0.6170813

Table 3: correlation of neighbouring pairs in the monthly average electricity demand in Victoria series.

par(mfrow=c(1,2))
hist(ts.ve, xlab="Monthly Average Electricity Demand (MWh)",
 main="Distribution of the Monthly \nAverage Electricity Demand series")
## Normality checks
qqnorm(ts.ve, main="Normal Q-Q plot of the Monthly \nAverage Electricity Demand series")

qqline(ts.ve, col = 2) # tails depart from normal, indicating non-normal
par(mfrow=c(1,1))

Figure 6: Histogram and Normal Q-Q plot of the monthly electricity demand in Victoria series.

Although a histogram of the monthly average demand series suggested the series was approximately symmetrical, a Q-Q plot of the series showed deviations from the normal in both tails, suggesting the series was not normally distributed. However, a Shapiro-Wilk test resulted in a p-value of 0.056, indicating the series was approximately normal at the α=0.05 level (Table 4).

shapiro.test(ts.ve)

    Shapiro-Wilk normality test

data:  ts.ve
W = 0.96623, p-value = 0.05575

Table 4: : Result of a Shapiro-Wilk test for normality on the monthly electricity demand in Victoria series

acf_pacf(ts.ve, "\nthe Monthly Electricity Demand series.")

Figure 7: ACF and PACF plots of the monthly electricity demand in Victoria series.

An ACF plot of the monthly average electricity demand series revealed a slowly decaying wave pattern with seasonal peaks, indicative of a seasonal trend, and a PACF plot of the series displayed a highly significant first PACF lag, indicative of an autoregressive process with a trend.

Finally, tests for stationarity confirmed the series was stationary, but only just (Table 5). An Augmented Dickey-Fuller test produced a p-value of 0.040 which was in the “doubt range” of 0.03 to 0.1. As a result, a Phillips-Perron test and a KPSS test were performed, resulting in p-values of 0.01 and 0.055, respectively, and indicating that there was sufficient evidence to conclude that the series was stationary at the α=0.05 level.

adf.test(ts.ve) #significant = stationary, doubt range: (0.03 - 0.1)

    Augmented Dickey-Fuller Test

data:  ts.ve
Dickey-Fuller = -3.5992, Lag order = 4, p-value = 0.03989
alternative hypothesis: stationary

pp.test(ts.ve) #significant = stationary

    Phillips-Perron Unit Root Test

data:  ts.ve
Dickey-Fuller Z(alpha) = -30.944, Truncation lag parameter = 3, p-value = 0.01
alternative hypothesis: stationary

kpss.test(ts.ve) #significant = non-stationary

    KPSS Test for Level Stationarity

data:  ts.ve
KPSS Level = 0.45204, Truncation lag parameter = 3, p-value = 0.05472

Table 5: Statistical tests for stationarity of the monthly electricity demand in Victoria series.

Exploring Transformation

Before model specification, a Box-Cox transformation was applied to the series to observe the impact of transformation on the series. A vector of candidate power transformation was required to obtain candidate lambda values for the Box-Cox transformation due to computational limitations. Through trial and error, limits of 4 and 5 were identified as producing a range that could capture the maximum log likelihood lambda value for the Box-Cox transformation. The plot of candidate lambda values for the Box-Cox transformation are shown in Figure 8, including the chosen lambda value for the Box-Cox transformation of 4.6 (the middle vertical line in Figure 8).

BC <- BoxCox.ar(ts.ve, lambda=seq(4, 5, 0.01))

BC$ci #Values of the first and third vertical lines

[1] 4 5

# To find the lambda value of the middle vertical line
lambda <- BC$lambda[which(max(BC$loglike) == BC$loglike)]
lambda

[1] 4.6

Figure 8: Log likelihood plot for the lambda values in the Box-Cox transformation of the monthly electricity demand in Victoria series.

With an appropriate lambda value identified, the Box-Cox transformation was applied to the electricity demand series and plotted (Figure 9). Visual inspection could not identify any clear improvement in the stationarity of the Box-Cox transformed series, compared with the un-transformed series.

veBC <- ((ts.ve^lambda) - 1) / lambda
plot(veBC, type="o", xlab="Time", ylab="Monthly Average Electricity Demand (MWh)",
 main = "Time series plot of Box-Cox transformed \nMonthly Average Electricity Demand series.")

Figure 9: Time series plot of the Box-Cox transformed monthly electricity demand in Victoria series.

Although an Augmented Dickey-Fuller Test of the Box-Cox transformed series (Table 6) resulted in a p-value in the “doubt range” (0.03 to 0.1), results from a Phillips-Perron unit root test and a KPSS test indicated the series was stationary (p-values of 0.01 and 0.08, respectively), with the KPSS result being less significant than the result from the KPSS test on the un-transformed series (Table 5), suggesting the transformation may have improved stationarity slightly.

adf.test(veBC) # significant = stationary, doubt range: (0.03 - 0.1)

    Augmented Dickey-Fuller Test

data:  veBC
Dickey-Fuller = -3.59, Lag order = 4, p-value = 0.04068
alternative hypothesis: stationary

pp.test(veBC) #significant = stationary

    Phillips-Perron Unit Root Test

data:  veBC
Dickey-Fuller Z(alpha) = -30.534, Truncation lag parameter = 3, p-value = 0.01
alternative hypothesis: stationary

kpss.test(veBC)

    KPSS Test for Level Stationarity

data:  veBC
KPSS Level = 0.39826, Truncation lag parameter = 3, p-value = 0.07791

Table 6: Results from tests of stationarity of the Box-Cox transformed monthly electricity demand series.

A Q-Q plot of the Box-Cox transformed series was then produced (Figure 10), which still showed deviations from the diagonal in both tails of the series, indicating the deviations had not been addressed by the transformation. Finally, a Shapiro-Wilk test on the Box-Cox transformed series resulted in a p-value of 0.004, indicating there was sufficient evidence to conclude that the Box-Cox transformed series was not approximately normal (Table 7).

As the improvement in stationarity from the Box-Cox transformation was very small, and the Box-Cox transformation resulted in a series that was no-longer approximately normal, the Box-Cox transformation was deemed inappropriate, and the un-transformed series would be used for model specification.

qqnorm(veBC, main="Normal Q-Q Plot for the \nBox-Cox transformed Monthly Electricity Demand series.")
qqline(veBC, col = 2)

Figure 10: Normal Q-Q Plot for the Box-Cox transformed monthly electricity demand series

shapiro.test(veBC)

    Shapiro-Wilk normality test

data:  veBC
W = 0.94446, p-value = 0.003759

Table 7: Result from a Shapiro-Wilk test on the Box-Cox transformed monthly electricity demand series.

Model Specification – SARIMA (Residual Approach)

As the series exhibited clear evidence of seasonality, SARIMA model specification was explored to observe how well a SARIMA model could capture the information in the monthly average electricity demand series. The residual approach was used for specifying SARIMA models for the series.

Seasonal Differencing

As the tests for stationarity indicated the monthly average electricity demand series was stationary, seasonal differencing was applied as a first step in SARIMA model specification. Figure 11: Residuals from the first seasonal difference SARIMA model, fit to the monthly electricity demand series.

m1.elec <- arima(ts.ve, order=c(0,0,0), #(p,d,q)
 seasonal=list(order=c(0,1,0), #(P,D,Q)
 period=12)) # s=12
res.m1 <- rstandard(m1.elec) # produce residuals from D=1 model
plot(res.m1, xlab="Time", ylab="Residuals", main="Time series plot of the residuals from 
SARIMA(0,0,0)x(0,1,0)_12 model.")

Figure 11: Residuals from the first seasonal difference SARIMA model, fit to the monthly electricity demand series.

Although there was still evidence of a trend in the residuals from the first seasonal difference model, the plot of residuals in Figure 11 indicated that seasonal differencing had brought the residuals closer to a white noise series.

acf_pacf(res.m1, "the residuals after \nfitting the first seasonal difference.")$acf

NULL

Figure 12: ACF and PACF plots of the first seasonal difference SARIMA model, applied to the monthly electricity demand series.

The ACF plot of the residuals from the first seasonal differenced model (Figure 12) showed no evidence of a slowly decaying seasonal pattern (significant lags at seasonal periods), indicating that first seasonal differencing had addressed the seasonality and further seasonal differencing was not necessary. The ACF plot also showed one significant seasonal lag (lag period 1), suggesting that the value of capital Q in the SARIMA model could be 1.

Meanwhile, the PACF plot showed one significant seasonal PACF lag (lag period 2), suggesting the value of capital P in the SARIMA model could be 1.

Evaluating Seasonal Parameters

With the seasonal differencing confirmed, and values of capital P and Q identified, a SARIMA(0,0,0)x(1,1,1)12 model was fit to the electricity demand series, with residuals plotted in Figure 13 and Figure 14.

# fit model with P=1, D=1, Q=1
m2.elec <- arima(ts.ve, order=c(0,0,0), #(p,d,q)
 seasonal=list(order=c(1,1,1), #(P,D,Q)
 period=12)) # s=12
res.m2 <- rstandard(m2.elec)
plot(res.m2, xlab="Time", ylab="Residuals", 
 main="Time series plot of the residuals from the second SARIMA model \n SARIMA(0,0,0)x(1,1,1)_12.")

Figure 13: Residuals from the second SARIMA model fit to the monthly electricity demand series.

The residuals from the second fitted SARIMA model (Figure 13) suggested that the trend observed in the residuals from the first fitted seasonal model had been reduced, but not entirely removed, indicating there was still work to be done.

acf_pacf(res.m2, "the residuals after \nfitting SARIMA(0,0,0)x(1,1,1)_12")

Figure 14: ACF and PACF plots of the residuals from the second SARIMA model applied to the monthly electricity demand series.

After fitting the SARIMA(0,0,0)x(1,1,1)12 model, there were no significant seasonal PACF lags in the PACF plot of the residuals (Figure 14), supporting a value of 1 for capital P.

Although there was still a significant seasonal ACF lag in the ACF plot of the residuals (Figure 14), there was also a slowly decaying pattern within the first period (between zero and lag period 1 of the ACF plot) which can be an indication of non-stationarity. As such, an Augmented Dickey Fuller test was conducted on the series (Table 8), resulting in a p-value of 0.216, indicating the series was not stationary.

As capital Q had already been specified, and the series was non-stationary, ordinary differencing would be applied to observe its impact before attempting to re-specify the value of capital Q.

adf.test(res.m2) #significant = stationary

    Augmented Dickey-Fuller Test

data:  res.m2
Dickey-Fuller = -2.8838, Lag order = 4, p-value = 0.2158
alternative hypothesis: stationary

Table 8: Result from an Augmented Dickey Fuller test of the residuals from the SARIMA(0,0,0)x(1,1,1)12 fitted model.

Ordinary Differencing

A SARIMA(0,1,0)x(1,1,1)12 model was fit to the electricity demand series, with residuals plotted in Figure 15.

m3.elec <- arima(ts.ve, order=c(0,1,0), #(p,d,q)
 seasonal=list(order=c(1,1,1), #(P,D,Q)
 period=12)) # s=12
res.m3 <- rstandard(m3.elec)
plot(res.m3, xlab="Time", ylab="Residuals", 
 main="Time series plot of the residuals from the third SARIMA model \n SARIMA(0,1,0)x(1,1,1)_12.")

Figure 15: Residuals from the third SARIMA model fit to the monthly electricity demand series. The residuals from the third SARIMA model fit to the electricity demand series showed no clear trend and was much closer to a white noise series.

ACF and PACF plots were produced separately to aid in visual inspection of the significant lags, and any significant seasonal lags that may be present.

acf(res.m3, lag.max = 50, main="ACF plot of the residuals from fitting the third SARIMA model \nSARIMA(0,1,0)x(1,1,1)_12.")$acf

, , 1

               [,1]
 [1,] -0.5308210642
 [2,]  0.1415684623
 [3,] -0.0108026966
 [4,] -0.1984253792
 [5,]  0.1529420880
 [6,] -0.0828465943
 [7,]  0.0318307741
 [8,]  0.1537205864
 [9,] -0.1289643114
[10,] -0.0105070881
[11,] -0.0569745873
[12,]  0.1058377207
[13,] -0.2361088351
[14,]  0.2394092349
[15,] -0.0907193598
[16,]  0.0359676039
[17,]  0.0288279059
[18,] -0.1219946303
[19,]  0.1028221519
[20,] -0.0328633481
[21,] -0.0912192990
[22,]  0.1897320009
[23,] -0.0575769929
[24,] -0.0206190057
[25,] -0.0426407124
[26,]  0.0692718078
[27,] -0.0727247285
[28,]  0.0320346636
[29,] -0.0173349686
[30,]  0.0686106853
[31,]  0.0688439704
[32,] -0.2094035858
[33,]  0.1375309536
[34,] -0.0965746566
[35,]  0.0397056746
[36,]  0.0315884651
[37,] -0.0297790426
[38,]  0.0594644266
[39,] -0.0398522627
[40,] -0.0090086898
[41,] -0.0179328898
[42,] -0.0193992948
[43,] -0.0237708406
[44,]  0.0364315564
[45,]  0.0587024317
[46,] -0.0360257610
[47,]  0.0453037920
[48,] -0.0909331566
[49,]  0.0463452859
[50,]  0.0007746394

Figure 16: ACF plot of the residuals from the third SARIMA model applied to the monthly electricity demand series.

The ACF plot of the residuals from the third SARIMA model (Figure 16) showed no evidence of a slowly decaying pattern in the first period, indicating that ordinary differencing had addressed the decaying pattern. Ordinary differencing had also lead to fewer significant lags around the first seasonal lag, making it clear that the value of capital Q did not need to be re-specified.

One significant lag was identified within the first period of the ACF plot (between zero and lag period 1), suggesting a value of 1 for “small q”.

pacf(res.m3, lag.max = 50, main="PACF plot of the residuals \nfrom fitting SARIMA(0,1,0)x(1,1,1)_12.")$acf

, , 1

              [,1]
 [1,] -0.530821064
 [2,] -0.195205903
 [3,] -0.035615490
 [4,] -0.280333531
 [5,] -0.153564767
 [6,] -0.111143284
 [7,] -0.098990683
 [8,]  0.128395797
 [9,]  0.082701810
[10,] -0.076244979
[11,] -0.155266152
[12,]  0.094910808
[13,] -0.276848656
[14,] -0.110452013
[15,] -0.046416114
[16,] -0.058438955
[17,] -0.069881576
[18,] -0.051688984
[19,] -0.023699229
[20,] -0.022148509
[21,] -0.100923854
[22,]  0.002191603
[23,]  0.084386934
[24,] -0.037338532
[25,] -0.116683733
[26,]  0.024164994
[27,] -0.006454089
[28,] -0.096504451
[29,] -0.061335276
[30,]  0.036987457
[31,]  0.140189597
[32,] -0.073156893
[33,]  0.017450711
[34,] -0.095488894
[35,]  0.031999754
[36,]  0.011156845
[37,] -0.026629571
[38,] -0.076179274
[39,]  0.039603414
[40,]  0.093518257
[41,] -0.032821560
[42,] -0.093024587
[43,] -0.121173196
[44,] -0.032201939
[45,] -0.087949176
[46,] -0.003318902
[47,]  0.032144011
[48,] -0.044253606
[49,] -0.018071927
[50,]  0.075539556

Figure 17: PACF plot of the residuals from the third SARIMA model applied to the monthly electricity demand series.

The PACF plot of the residuals from the third SARIMA model showed two significant lags within the first period (zero to lag period 1), so a value of 2 was proposed for “small p”.

Evaluating Ordinary Parameters

With proposed values of small p and small q identified, a SARIMA(2,1,1)x(1,1,1)12 model was fit to the electricity demand series, with residuals plotted in Figure 18.

m4.elec <- arima(ts.ve, order=c(2,1,1), #(p,d,q)
 seasonal=list(order=c(1,1,1), #(P,D,Q)
 period=12)) # s=12
res.m4 <- rstandard(m4.elec)
plot(res.m4, xlab="Time", ylab="Residuals", 
 main="Time series plot of the residuals from the fourth SARIMA model \n SARIMA(2,1,1)x(1,1,1)_12.")

Figure 18: Residuals from the fourth SARIMA model fit to the monthly electricity demand series. The residuals from the fourth SARIMA model fit to the electricity demand series appeared very similar to a white noise series, suggesting the fourth SARIMA model fitted to the electricity demand series had resulted in white noiseresiduals.

acf_pacf(res.m4, "the residuals after fitting \nSARIMA(2,1,1)x(1,1,1)_12.")

Figure 19: ACF and PACF plots of the residuals from the fourth SARIMA model applied to the monthly electricity demand series.

There were no significant lags in the ACF plot of the residuals from the fourth SARIMA model fit to the electricity demand series (Figure 19).

Although there was a significant PACF lag in the residuals of the fourth SARIMA model, it was not a seasonal lag, and occurred outside of the first period (after lag period 1), so it was ignored, and the fourth SARIMA model was accepted.

Further supporting this decision were tests of stationary (Table 9), that confirmed the residuals from the fourth SARIMA model were stationary at the α=0.05 level, with p-values of 0.05 from the ADF test (doubt range), 0.01 for the Phillips-Perron unit root test, and 0.1 for the KPSS test. Table 9: Results from tests for stationarity of the residuals from the fourth SARIMA model fit to the monthly electricity demand series.

adf.test(res.m4) #significant = stationary

    Augmented Dickey-Fuller Test

data:  res.m4
Dickey-Fuller = -3.4664, Lag order = 4, p-value = 0.05234
alternative hypothesis: stationary

pp.test(res.m4)

    Phillips-Perron Unit Root Test

data:  res.m4
Dickey-Fuller Z(alpha) = -70.485, Truncation lag parameter = 3, p-value = 0.01
alternative hypothesis: stationary

kpss.test(res.m4)

    KPSS Test for Level Stationarity

data:  res.m4
KPSS Level = 0.086709, Truncation lag parameter = 3, p-value = 0.1

Table 9: Results from tests for stationarity of the residuals from the fourth SARIMA model fit to the monthly electricity demand series.

Model Specification - EACF & BIC

When the fourth SARIMA model was fit to the electricity demand series, the residuals appeared similar to a white noise series, and there were no significant ACF or PACF lags in the residuals, indicating there was no information left in the residuals. Because there was no information left in the residuals from the fourth SARIAM model, the residuals from the third SARIMA model, SARIMA(0,1,0)x(1,1,1)12 were used for model specification with EACF and BIC.

EACF

To propose possible values of “small p” and “small q” within the SARIMA(p,d,q)x(1,1,1)12 model, the extended ACF function was applied to the residuals from the third SARIMA model identified in the residual approach with the output from the function shown in Table 10.

eacf(res.m3)

AR/MA
  0 1 2 3 4 5 6 7 8 9 10 11 12 13
0 x o o o o o o o o o o  o  o  o 
1 x o o x o o o o o o o  o  o  o 
2 o o x x o o o o o o o  o  o  o 
3 o x x x o o o o o o o  o  o  o 
4 x o o o o o o o o o o  o  x  o 
5 x o o x o o o o o o o  o  x  o 
6 x o o x o o o o o o o  o  x  o 
7 x x o o o o o o o o o  o  x  o 

Table 10: EACF plot of the residuals from the third SARIMA model fit to the electricity demand series.

In the EACF plot, possible values of “small p” and “small q” are identified by finding the top-left most “o” and taking the corresponding co-ordinates from the AR and MA axes. The co-ordinates of neighbouring “o”s are then used to propose additional models from the EACF plot.

In the EACF plot (Table 10), the most top-left point was identified as (0,1), and the proposed models were;

• { SARIMA(0,1,1)x(1,1,1)_12,
• SARIMA(1,1,1)x(1,1,1)_12,
• SARIMA(0,1,2)x(1,1,1)_12,
• SARIMA(1,1,2)x(1,1,1)_12 }

BIC

Values for “small p” and “small q were also proposed using the Bayesian Information Criterion, or “BIC”, with the output from the BIC plot of residuals from the third SARIMA model shown in Figure 20.

res <- armasubsets(y=res.m3, nar=10, nma=10, y.name='p', ar.method='ols')

Reordering variables and trying again:

plot(res)

Figure 20: : BIC plot of the residuals from the third SARIMA model fit to the electricity demand series.

In the BIC plot, values of “small p” and “small q” are identified using the shaded squares in the top rows of the plot, as the models in these rows have values of “small p” and “small q” that result in the lowest BIC score(s). In Figure 20, the values of “small p” correspond to columns labelled “p-lag n”, while values of “small q” correspond to columns labelled “error-lag n”.

Using the BIC plot in Figure 20, the lowest BIC score achieved was -12, corresponding to the top 3 rows. From these rows, the models proposed by the BIC plot were;

• { SARIMA(0,1,2)x(1,1,1)_12,
• SARIMA(0,1,3)x(1,1,1)_12,
• SARIMA(0,1,10)x(1,1,1)_12 }.

All 3 models were supported by multiple rows in the BIC plot, with the SARIMA(0,1,2)x(1,1,1)_12 and SARIMA(0,1,10)x(1,1,1)_12 models supported by the top 2 rows of the BIC plot.

Final set of models from EACF and BIC model specification

Consolidating the models proposed by the EACF and BIC tools, the final set of proposed models was;

• { SARIMA(0,1,1)x(1,1,1)_12,
• SARIMA(1,1,1)x(1,1,1)_12,
• SARIMA(0,1,2)x(1,1,1)_12,
• SARIMA(1,1,2)x(1,1,1)_12,
• SARIMA(0,1,3)x(1,1,1)_12,
• SARIMA(0,1,10)x(1,1,1)_12 }

Model Fitting

Each of the 6 models proposed using EACF and BIC tools was fit to the electricity demand series and a coefficient significant test was used to identify the coefficients that were significant at the α=0.05 level under Maximum Likelihood Estimation (ML) and Conditional Sum of Squares (CSS) methods. The significant coefficients identified under each method are summarised in Table 11.

SARIMA model

Source: data provided for assignment by RMIT.

Table 11: Summary table of coefficients identified as significant when fitting each SARIMA model proposed by EACF and BIC tools. * = coefficient was only significant at the α=0.1 level.

It was clear from the results in Table 11 that there were more significant coefficients under the “CSS” method compared to the “ML” method. Subsequent coefficient tests were also performed using the combined “CSS-ML” method to observe whether the combined method could bring the “ML” results closer to the results from the “CSS” method. In all coefficient tests under the combined “CSS-ML” method, the combined method made little or no difference in the significant coefficients compared to the results under the “ML” method.

With significant coefficients identified, AIC and BIC scores were calculated for each model and sorted to identify the model with the best (lowest) AIC and BIC scores. The results from the function created to sort AIC and BIC scores are shown in Table 12.

m011.elec <- arima(ts.ve, order=c(0,1,1), 
 seasonal=list(order=c(1,1,1), period=12), method="ML")
coeftest(m011.elec) #significant: ma1

z test of coefficients:

      Estimate Std. Error z value  Pr(>|z|)    
ma1  -0.809874   0.155859 -5.1962 2.034e-07 ***
sar1 -0.016535   0.308264 -0.0536    0.9572    
sma1 -0.831536   0.630801 -1.3182    0.1874    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

m011.elecCSS <- arima(ts.ve, order=c(0,1,1), 
 seasonal=list(order=c(1,1,1), period=12), method="CSS")
coeftest(m011.elecCSS) #significant: ma1, sar1, sma1 (ALL)

z test of coefficients:

     Estimate Std. Error z value Pr(>|z|)    
ma1  -0.69295    0.11182 -6.1971 5.75e-10 ***
sar1 -0.31489    0.13453 -2.3407  0.01925 *  
sma1 -0.32848    0.15480 -2.1220  0.03384 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

m011.elecCSSML <- arima(ts.ve, order=c(0,1,1), 
 seasonal=list(order=c(1,1,1), period=12), method="CSS-ML")
coeftest(m011.elecCSSML) #significant: ma1

z test of coefficients:

      Estimate Std. Error z value Pr(>|z|)    
ma1  -0.809923   0.155858 -5.1966 2.03e-07 ***
sar1 -0.016598   0.308228 -0.0539   0.9571    
sma1 -0.831455   0.630468 -1.3188   0.1872    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#..........................
# SARIMA(1,1,1)x(1,1,1)_12
m111.elec <- arima(ts.ve, order=c(1,1,1), 
 seasonal=list(order=c(1,1,1), period=12), method="ML")
coeftest(m111.elec) #significant: ma1

z test of coefficients:

      Estimate Std. Error z value  Pr(>|z|)    
ar1   0.173019   0.141183  1.2255    0.2204    
ma1  -0.999829   0.207679 -4.8143 1.477e-06 ***
sar1  0.028995   0.186409  0.1555    0.8764    
sma1 -0.996630   1.173327 -0.8494    0.3957    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

m111.elecCSS <- arima(ts.ve, order=c(1,1,1), 
 seasonal=list(order=c(1,1,1), period=12), method="CSS")
coeftest(m111.elecCSS) #significant: ma1, sma1. sar1 slightly significant

z test of coefficients:

      Estimate Std. Error z value  Pr(>|z|)    
ar1   0.002163   0.207080  0.0104 0.9916661    
ma1  -0.560411   0.166663 -3.3625 0.0007723 ***
sar1 -0.248362   0.134108 -1.8520 0.0640328 .  
sma1 -0.386302   0.156249 -2.4724 0.0134226 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

m111.elecCSSML <- arima(ts.ve, order=c(1,1,1), 
 seasonal=list(order=c(1,1,1), period=12), method="CSS-ML")
coeftest(m111.elecCSSML) #significant: ma1

z test of coefficients:

      Estimate Std. Error z value  Pr(>|z|)    
ar1   0.172988   0.141177  1.2253    0.2205    
ma1  -0.999915   0.207794 -4.8120 1.494e-06 ***
sar1  0.028972   0.186266  0.1555    0.8764    
sma1 -0.998376   1.178446 -0.8472    0.3969    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#..........................
# SARIMA(0,1,2)x(1,1,1)_12
m012.elec <- arima(ts.ve, order=c(0,1,2), 
 seasonal=list(order=c(1,1,1), period=12), method="ML")
coeftest(m012.elec) #significant: ma1

z test of coefficients:

       Estimate Std. Error z value  Pr(>|z|)    
ma1  -0.8068348  0.1682033 -4.7968 1.612e-06 ***
ma2  -0.0768209  0.1840409 -0.4174    0.6764    
sar1 -0.0054301  0.2963605 -0.0183    0.9854    
sma1 -0.8959873  1.0500212 -0.8533    0.3935    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

m012.elecCSS <- arima(ts.ve, order=c(0,1,2), 
 seasonal=list(order=c(1,1,1), period=12), method="CSS")
coeftest(m012.elecCSS) #significant: ma1, sar1, sma1 (3/4)

z test of coefficients:

      Estimate Std. Error z value  Pr(>|z|)    
ma1  -0.622411   0.148536 -4.1903 2.786e-05 ***
ma2  -0.086941   0.132888 -0.6542   0.51295    
sar1 -0.281615   0.141611 -1.9887   0.04674 *  
sma1 -0.362009   0.162519 -2.2275   0.02591 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

m012.elecCSSML <- arima(ts.ve, order=c(0,1,2), 
 seasonal=list(order=c(1,1,1), period=12), method="CSS-ML")
coeftest(m012.elecCSSML) #significant: ma1

z test of coefficients:

       Estimate Std. Error z value  Pr(>|z|)    
ma1  -0.8068450  0.1682116 -4.7966 1.614e-06 ***
ma2  -0.0768218  0.1840440 -0.4174    0.6764    
sar1 -0.0054236  0.2962780 -0.0183    0.9854    
sma1 -0.8959811  1.0494303 -0.8538    0.3932    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#..........................
# SARIMA(1,1,2)x(1,1,1)_12
m112.elec <- arima(ts.ve, order=c(1,1,2), 
 seasonal=list(order=c(1,1,1), period=12), method="ML")
coeftest(m112.elec) #significant: ar1, ma1 (2/5)

z test of coefficients:

      Estimate Std. Error z value  Pr(>|z|)    
ar1   0.683372   0.308947  2.2119   0.02697 *  
ma1  -1.494272   0.360186 -4.1486 3.345e-05 ***
ma2   0.494547   0.350503  1.4110   0.15826    
sar1  0.073609   0.232386  0.3168   0.75143    
sma1 -0.964360   1.657483 -0.5818   0.56069    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

m112.elecCSS <- arima(ts.ve, order=c(1,1,2), 
 seasonal=list(order=c(1,1,1), period=12), method="CSS")
coeftest(m112.elecCSS) #significant: ar1, ma1, ma2, sar1. (sma1 slightly)

z test of coefficients:

       Estimate Std. Error   z value  Pr(>|z|)    
ar1  -0.9389813  0.0017596 -533.6291 < 2.2e-16 ***
ma1   0.5081401  0.0966007    5.2602 1.439e-07 ***
ma2  -0.7288425  0.1082163   -6.7351 1.639e-11 ***
sar1 -0.3284227  0.1304098   -2.5184   0.01179 *  
sma1 -0.2837530  0.1575636   -1.8009   0.07172 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

m112.elecCSSML <- arima(ts.ve, order=c(1,1,2), 
 seasonal=list(order=c(1,1,1), period=12), method="CSS-ML")
coeftest(m112.elecCSSML) #significant: ar1, ma2

z test of coefficients:

      Estimate Std. Error z value  Pr(>|z|)    
ar1  -0.960448   0.232988 -4.1223 3.751e-05 ***
ma1   0.141024   0.263963  0.5343   0.59316    
ma2  -0.760185   0.236213 -3.2182   0.00129 ** 
sar1 -0.012274   0.317402 -0.0387   0.96915    
sma1 -0.851445   0.752609 -1.1313   0.25792    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#..........................
# SARIMA(0,1,3)x(1,1,1)_12
m013.elec <- arima(ts.ve, order=c(0,1,3), 
 seasonal=list(order=c(1,1,1), period=12), method="ML")
coeftest(m013.elec) #significant: ma1, ma3 (2/5)

z test of coefficients:

      Estimate Std. Error z value Pr(>|z|)    
ma1  -0.875865   0.172304 -5.0833 3.71e-07 ***
ma2   0.164355   0.189217  0.8686  0.38506    
ma3  -0.288444   0.146850 -1.9642  0.04951 *  
sar1  0.098258   0.182277  0.5391  0.58985    
sma1 -0.999449   0.816219 -1.2245  0.22077    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

m013.elecCSS <- arima(ts.ve, order=c(0,1,3), 
 seasonal=list(order=c(1,1,1), period=12), method="CSS")
coeftest(m013.elecCSS) #significant: ma1, sar1, sma1 (3/5)

z test of coefficients:

       Estimate Std. Error z value  Pr(>|z|)    
ma1  -0.6222375  0.1514560 -4.1084 3.985e-05 ***
ma2  -0.0049533  0.2101461 -0.0236   0.98119    
ma3  -0.0940025  0.1822346 -0.5158   0.60597    
sar1 -0.2846956  0.1406170 -2.0246   0.04291 *  
sma1 -0.3392088  0.1711043 -1.9825   0.04743 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

m013.elecCSSML <- arima(ts.ve, order=c(0,1,3), 
 seasonal=list(order=c(1,1,1), period=12), method="CSS-ML")
coeftest(m013.elecCSSML) #significant: ma1, ma3

z test of coefficients:

      Estimate Std. Error z value Pr(>|z|)    
ma1  -0.875760   0.172265 -5.0838  3.7e-07 ***
ma2   0.164216   0.189229  0.8678  0.38549    
ma3  -0.288303   0.146879 -1.9629  0.04966 *  
sar1  0.098114   0.182265  0.5383  0.59037    
sma1 -0.999639   0.818805 -1.2209  0.22214    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#..........................
# SARIMA(0,1,10)x(1,1,1)_12
m0110.elec <- arima(ts.ve, order=c(0,1,10), 
 seasonal=list(order=c(1,1,1), period=12), method="ML")
coeftest(m0110.elec) #significant: ma1, ma5, ma9 (ma6 slightly)

z test of coefficients:

       Estimate Std. Error z value  Pr(>|z|)    
ma1  -0.8088721  0.1639771 -4.9328 8.104e-07 ***
ma2   0.1711300  0.1577841  1.0846  0.278106    
ma3  -0.3170106  0.2058961 -1.5397  0.123643    
ma4  -0.0526888  0.1998735 -0.2636  0.792080    
ma5   0.3353477  0.1667891  2.0106  0.044367 *  
ma6  -0.3136305  0.1687111 -1.8590  0.063030 .  
ma7   0.2440080  0.1634153  1.4932  0.135391    
ma8   0.1535326  0.1776646  0.8642  0.387494    
ma9  -0.5809941  0.1819398 -3.1933  0.001406 ** 
ma10  0.1691848  0.1601178  1.0566  0.290682    
sar1 -0.0036632  0.3719782 -0.0098  0.992143    
sma1 -0.6717786  0.4118726 -1.6310  0.102883    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

m0110.elecCSS <- arima(ts.ve, order=c(0,1,10), 
 seasonal=list(order=c(1,1,1), period=12), method="CSS")
coeftest(m0110.elecCSS) #significant: ma1, ma2, ma3, ma5, ma9, sar1, sma1 (ma8 & ma10 slightly)

z test of coefficients:

      Estimate Std. Error z value Pr(>|z|)    
ma1  -0.422340   0.143916 -2.9346 0.003339 ** 
ma2   0.362223   0.137913  2.6265 0.008628 ** 
ma3  -0.567730   0.138874 -4.0881 4.35e-05 ***
ma4   0.029361   0.174006  0.1687 0.866005    
ma5   0.466174   0.181996  2.5615 0.010424 *  
ma6  -0.300337   0.197932 -1.5174 0.129172    
ma7   0.175793   0.206825  0.8500 0.395349    
ma8   0.322408   0.176227  1.8295 0.067324 .  
ma9  -0.413695   0.173925 -2.3786 0.017379 *  
ma10 -0.320318   0.165468 -1.9358 0.052889 .  
sar1 -0.083947   0.026625 -3.1529 0.001616 ** 
sma1 -0.345907   0.162045 -2.1346 0.032790 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

m0110.elecCSSML <- arima(ts.ve, order=c(0,1,10), 
 seasonal=list(order=c(1,1,1), period=12), method="CSS-ML")
coeftest(m0110.elecCSSML) #significant: ma1, ma5, ma9 (ma6 slightly)

z test of coefficients:

       Estimate Std. Error z value  Pr(>|z|)    
ma1  -0.8088367  0.1639624 -4.9331 8.095e-07 ***
ma2   0.1710508  0.1577668  1.0842  0.278276    
ma3  -0.3169169  0.2058768 -1.5394  0.123718    
ma4  -0.0527278  0.1998611 -0.2638  0.791917    
ma5   0.3352861  0.1667693  2.0105  0.044381 *  
ma6  -0.3135734  0.1687031 -1.8587  0.063066 .  
ma7   0.2439649  0.1633898  1.4931  0.135399    
ma8   0.1536371  0.1776375  0.8649  0.387099    
ma9  -0.5810232  0.1819189 -3.1939  0.001404 ** 
ma10  0.1691395  0.1601129  1.0564  0.290796    
sar1 -0.0037158  0.3719440 -0.0100  0.992029    
sma1 -0.6716938  0.4117835 -1.6312  0.102852    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# ========================= #
sort.score(AIC(m011.elec, m111.elec, m012.elec, m112.elec, m013.elec, m0110.elec), 
 score = "aic")

sort.score(BIC(m011.elec, m111.elec, m012.elec, m112.elec, m013.elec, m0110.elec), 
 score = "bic")

Table 12: Sorted AIC and BIC scores from the residuals of the models proposed by EACF and BIC tools, fit to the electricity demand series.

SARIMA(0,1,1)x(1,1,1)_12 was identified as the “best” model as it had the lowest AIC and BIC scores of all proposed models. This model was also supported by the coefficient tests, as it was the only model where all model coefficients were significant under the CSS method. All other models had at least one coefficient that was not significant when tested using the CSS method.

Overfitting

With the SARIMA(0,1,1)x(1,1,1)_12 model identified as the “best” model, overfitting was conducted to confirm that the coefficients of additional model parameters (“small p” and “small q”) were insignificant (overfitted).

The candidate overfitting models for SARIMA(0,1,1)x(1,1,1)_12 were SARIMA(1,1,1)x(1,1,1)_12 and SARIMA(0,1,2)x(1,1,1)_12. These models had been evaluated during the model fitting stage but were re-assessed to confirm that the overfitted coefficients were insignificant and that nothing had been overlooked.

The coefficient test results from SARIMA(1,1,1)x(1,1,1)_12 and SARIMA(0,1,2)x(1,1,1)_12 are shown in Table 13 and Table 14, respectively.

# Over fitting - SARIMA(1,1,1)x(1,1,1)_12
m111.elec <- arima(ts.ve, order=c(1,1,1), 
 seasonal=list(order=c(1,1,1), period=12), method="ML")
coeftest(m111.elec) #significant: ma1 (ar1 NOT significant)

z test of coefficients:

      Estimate Std. Error z value  Pr(>|z|)    
ar1   0.173019   0.141183  1.2255    0.2204    
ma1  -0.999829   0.207679 -4.8143 1.477e-06 ***
sar1  0.028995   0.186409  0.1555    0.8764    
sma1 -0.996630   1.173327 -0.8494    0.3957    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

m111.elecCSS <- arima(ts.ve, order=c(1,1,1), 
 seasonal=list(order=c(1,1,1), period=12), method="CSS")
coeftest(m111.elecCSS) #significant: ma1, sma1. sar1 slightly significant (ar1 NOT significant)

z test of coefficients:

      Estimate Std. Error z value  Pr(>|z|)    
ar1   0.002163   0.207080  0.0104 0.9916661    
ma1  -0.560411   0.166663 -3.3625 0.0007723 ***
sar1 -0.248362   0.134108 -1.8520 0.0640328 .  
sma1 -0.386302   0.156249 -2.4724 0.0134226 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#..........................
# Over fitting - SARIMA(0,1,2)x(1,1,1)_12
m012.elec <- arima(ts.ve, order=c(0,1,2), 
 seasonal=list(order=c(1,1,1), period=12), method="ML")
coeftest(m012.elec) #significant: ma1 (ma2 NOT significant)

z test of coefficients:

       Estimate Std. Error z value  Pr(>|z|)    
ma1  -0.8068348  0.1682033 -4.7968 1.612e-06 ***
ma2  -0.0768209  0.1840409 -0.4174    0.6764    
sar1 -0.0054301  0.2963605 -0.0183    0.9854    
sma1 -0.8959873  1.0500212 -0.8533    0.3935    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

m012.elecCSS <- arima(ts.ve, order=c(0,1,2), 
 seasonal=list(order=c(1,1,1), period=12), method="CSS")
coeftest(m012.elecCSS) #significant: ma1, sar1, sma1 (ma2 NOT significant)

z test of coefficients:

      Estimate Std. Error z value  Pr(>|z|)    
ma1  -0.622411   0.148536 -4.1903 2.786e-05 ***
ma2  -0.086941   0.132888 -0.6542   0.51295    
sar1 -0.281615   0.141611 -1.9887   0.04674 *  
sma1 -0.362009   0.162519 -2.2275   0.02591 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Table 13: Coefficient test results for the overfitted SARIMA(1,1,1)x(1,1,1)_12 model.

The coefficient tests in Table 13 indicated that the overfitted coefficient “ar1” was insignificant under both ML and CSS methods, confirming that a value of zero for “small p” was appropriate.

As the additional coefficients in the overfitted models were both insignificant, the SARIMA(0,1,1)x(1,1,1)_12 model was confirmed as an appropriate model for the electricity demand series.

sort.score(AIC(m011.elec, m111.elec, m012.elec, m112.elec, m013.elec, m0110.elec), 
 score = "aic")

sort.score(BIC(m011.elec, m111.elec, m012.elec, m112.elec, m013.elec, m0110.elec), 
 score = "bic")

Diagnostic Checking

To inspect the residuals from the SARIMA(0,1,1)x(1,1,1)12 model fit to the electricity demand series, the residuals were standardized and plotted in Figure 21.

st.res.m011 <- rstandard(m011.elec)
par(mar=c(5,6,4,1)+.1)
plot(st.res.m011, type="o", xlab="Time", ylab="Standardised Residuals",
 main="Standardised residuals from the SARIMA(0,1,1)x(1,1,1)_12 Model.")
abline(h=0)

Source: data provided for assignment by RMIT.

Figure 21: Time series plot of the standardized residuals from fitting the SARIMA(0,1,1)x(1,1,1)12 model to the electricity demand series. Markers added to aid in identifying possible outliers.

Although there were no major irregularities in the standardized residuals from the fitted SARIMA model, the standardized residuals for December 2016 and December 2017 (marked in red in Figure 21) looked like possible outliers, so further investigation was performed.

par(mar=c(5,6,4,2)+.1) #set plot window margins
acf(st.res.m011, lag.max=50, 
 main="ACF plot of the standardized residuals from the SARIMA(0,1,1)x(1,1,1)_12 Model.")

Figure 22: ACF plot of the standardized residuals from fitting the SARIMA(0,1,1)x(1,1,1)12 model to the electricity demand series.

There were no significant autocorrelation lags in the ACF plot of standardized residuals (Figure 22), indicating there was no evidence of autocorrelation in the standardized residuals from the fitted model.

#LBQPlot(st.res.m011, lag.max=30) ## function does not work 11/06/2023

Source: data provided for assignment by RMIT.

Figure 23: Ljung-Box test of the standardized residuals from fitting the SARIMA(0,1,1)x(1,1,1)12 model to the electricity demand series.

All points in the Ljung-Box test were above the 0.05 reference line, indicating that there was no issue with independence of errors in the standardized residuals (Figure 23).

par(mfrow=c(1,2))
qqnorm(st.res.m011, main="Q-Q plot of the standardized residuals: \nSARIMA(0,1,1)x(1,1,1)_12 Model.")
qqline(st.res.m011, col = 2)

hist(st.res.m011, xlab="Standardized residuals", 
 main="Standardized residuals from the \nSARIMA(0,1,1)x(1,1,1)_12 Model.",
 breaks = 24)
par(mfrow=c(1,1))

Figure 24: Q-Q plot and histogram of the standardized residuals from the fitted SARIMA(0,1,1)x(1,1,1)12 model.

shapiro.test(st.res.m011)

    Shapiro-Wilk normality test

data:  st.res.m011
W = 0.9518, p-value = 0.009016

Table 15: Result from a Shapiro-Wilk test on the standardized residuals the fitted SARIMA(0,1,1)x(1,1,1)12 model.

The Q-Q plot of the standardized residuals showed deviations from the normal in both tails of the residuals (Figure 24) and a Shapiro-Wilk test produced a p-value of 0.009 (Table 15) indicating that the standardized residuals may not be approximately normally distributed. However, the histogram of standardized residuals (Figure 24) appeared approximately symmetric, and inspection of the standardized residuals revealed that no standardized residuals were greater than +/- 3 (Table 16).

min(st.res.m011)

[1] -2.944563

max(st.res.m011)

[1] 2.831909

Table 16: Minimum and maximum standardized residual from the fitted SARIMA(0,1,1)x(1,1,1)12 model. As there were no significant points in the Ljung-Box test, no significant lags in the ACF plot, and the histogram of standardized residuals appeared approximately symmetric with no values greater than +/- 3, the suspected outliers were deemed non-problematic, and the model was accepted.

Results

Forecasting

With the SARIMA(0,1,1)x(1,1,1)_12 model identified as the most appropriate model for the Victorian electricity demand series, the model was used to forecast the average monthly electricity demand for the next 10 months. The historical data, forecasted points, and forecast limits are shown in Figure 25.

m011.elecCSSA <- Arima(ts.ve, order=c(0,1,1),
 seasonal=list(order=c(1,1,1), period=12),method = "CSS")
summary(m011.elecCSSA)

Series: ts.ve 
ARIMA(0,1,1)(1,1,1)[12] 

Coefficients:
          ma1     sar1     sma1
      -0.6929  -0.3149  -0.3285
s.e.   0.1118   0.1345   0.1548

sigma^2 = 12590667:  log likelihood = -552.01

Training set error measures:
                   ME     RMSE      MAE       MPE     MAPE      MASE       ACF1
Training set 300.0969 3116.537 1898.403 0.1955092 1.618582 0.4883199 0.04948418

frc <- forecast(m011.elecCSSA, h=10)
par(mfrow=c(1,1))
plot(frc, xlab="Month", ylab="Monthly Average Electricity Demand (MWh)",
main="Predictions for the monthly average Victorian Electricity Demand 
 for the next 10 months using the SARIMA(0,1,1)x(1,1,1)_12 Model.")

Figure 25: Time series plot of the historical monthly average electricity demand in Victoria, with forecasted values for the next 10 months, including forecast limits, based on the SARIMA(0,1,1)x(1,1,1)12 model.

frc

Table 17: Results from forecasting the next 10 months of average electricity demand in Victoria, using the SARIMA(0,1,1)x(1,1,1)12 model.

The forecasted values for electricity demand looked reasonable when compared to the historical data, including high demand in June, July, and August, which were identified as months that consistently had high demand for electricity in Victoria. The forecasts also appeared to do a good job of capturing the seasonality and the trend in the series.

The forecast limits, marked by light and dark grey in Figure 25, also looked reasonable, as the upper limits were not greater than historical highs, and the lower limits (although lower than any historical point) did not look unreasonable in the context of the downward trend in electricity demand.

Finally, the addition of the forecast limits appeared to improve the perceived precision and trustworthiness of the forecasts by providing a range of values the forecasts could fall within.

Summary & Conclusion

The objective of the project was to understand the “Daily Electricity Demand” dataset and identify the “best” model to accurately forecast electricity demand in Victoria.

Although the data had to be aggregated to average monthly demand, aggregating the data did aid in visual inspection of the series, and reduced the impact of days with unusually high demand resulting from unique circumstances. Daily forecasts could also prove less valuable than a monthly average, as many factors can influence the demand for electricity on a given day and a daily forecast could lead to a lack of confidence in forecasts when they cannot account for these many factors.

After model specification, the SARIMA(0,1,1)x(1,1,1)_12 model was identified as the “best” of the proposed models and appeared to do a good job of capturing both the seasonality and the trend in the monthly average electricity demand series. That all coefficients in the SARIMA(0,1,1)x(1,1,1)_12 model were significant was reassuring, as was diagnostic checking, which did not identify any major irregularities in the standardized residuals from the fitted model.

Forecasts produced using the SARIMA(0,1,1)x(1,1,1)_12 model appeared reasonable, and the forecast limits were reassuring.

Finally, further research and analysis is recommended to explore modelling of ultra-high frequency data to produce a model that can predict daily electricity demand, possibly through the use of Neural Networks (Crone, 2009).

References

Crone, S. F. (2009, June). Input-variable specification for neural networks-an analysis of forecasting low and high time series frequency. 2009 International Joint Conference on Neural Networks (pp. 619-626). Atlanta, GA: IEEE. doi:10.1109/IJCNN.2009.5179046

Kozlov, A. (2020). Daily Electricity Price and Demand Data. Retrieved from Kaggle: https://www.kaggle.com/datasets/aramacus/electricity-demand-in-victoria-australia

Victorian Government. (2022, 05 28). Practical Information - Weather. Retrieved from Visit Victoria: https://www.visitvictoria.com/practical-information/melbourne-weather