Names of students:

Allen Yu

Qianyun Sang

---

Table of Contents

Project description
Part A: Time plot, trying transformations/adjustments, creating training and test set
Part B
Introduction & Approach
Exploratory Data Analysis
Potential causes of patterns
Exponential Smoothing Model
ARIMA model
Comparison of models
Forecasts
Discussion

Appendix A: ETS Models
Appendix B: ARIMA Models

---

Project description

Objective

The main goal of this project is to develop a model to forecast monthly energy use for the Vancouver International Airport (YVR).

Introduction

Budget planning at YVR requires forecasting the cost of energy to operate the airport. An accurate forecast could help the YVR representatives negotiate more favourable contracts with energy suppliers.

Data Available

Data was taken from 14 years of records from the YVR database. This data set includes information on date, energy use, temperature, terminal area, and number of passengers. Data have been aggregated by month into 168 records (14 years with 12 months each).

#	Variable	Definition
1	month	Month and year, e.g.: Nov-98
2	energy	Energy use measured in thousands of kilowatt hours (kWh)
3	mean.temp	Mean monthly temperature outside (degrees Celsius)
4	total.area	Total area of all terminals (sq. m.)
5	total.passengers	Total number of passengers in thousands
6	domestic.passengers	Total number of domestic passengers (traveling within Canada) in thousands
7	US.passengers	Total number of passengers traveling between Canada and the US in thousands
8	international.passengers	Total number of passengers traveling between YVR and countries other than Canada/US

In this project, you will develop appropriate models, compare the models, and discuss advantages and limitations of each model. You should select the best model and use it to provide monthly forecasts for energy use for the next three years (January 2011 through December 2013).

---

Part A (3.5 marks)

Time plot, trying transformations/adjustments, creating training and test set

See Canvas submission called "Project: Part A" for detailed questions and points breakdown.

Import the data

Create a time series object of the energy use. See the steps on Canvas if you aren't sure how to create a time series object in R.

# Load required libraries
library(tseries)
library(forecast)
library(ggplot2)

data <- read.csv("Energy use at YVR.csv")

Registered S3 method overwritten by 'quantmod':
  method            from
  as.zoo.data.frame zoo 

Time plot (0.5 marks)

Create a time plot of the electricity use data. Prepare this plot as if for publication and upload the image to Canvas. You do not have to describe the plot.

# Convert 'month' column to Date format
data$month <- as.Date(paste0("01-", data$month), format="%d-%b-%y")

ts_energy <- ts(data$energy, start=c(1997,1), frequency=12)

# Time plot
plot(ts_energy, main="Time Plot of Electricity Use at YVR", xlab="Year", ylab="Energy Use (MWh)", col="blue", lwd=2)

Box-Cox transformation (1 mark)

Write the code to do a Box-Cox transformation of the electricity use data (this will be submitted as part of the last question for Part A of the project).

Create a time plot of the transformed data. Prepare this plot as if for publication and upload the image.

Explain the purpose of using the Box-Cox transformation and if it was effective for this dataset (5 sentences max).

# Box-Cox transformation
lambda_bc <- BoxCox.lambda(ts_energy)
ts_energy_bc <- BoxCox(ts_energy, lambda_bc)

# Time plot of transformed data
plot(ts_energy_bc, main="Time Plot of Box-Cox Transformed Electricity Use", xlab="Year", ylab="Transformed Energy Use", col="green", lwd=2)

Explanation of purpose of Box-Cox transformation

• The Box-Cox transformation was applied to stabilize variance and reduce fluctuations in YVR’s monthly energy use data.
• The initial time plot showed increasing variance and irregular spikes, indicating potential heteroskedasticity.
• The transformation successfully smoothed fluctuations and provided a more consistent and linear trend.
• this adjustment improves the suitability of the data for accurate time series forecasting, making it easier to detect underlying patterns.
• With the transformation, YVR’s forecasting model is expected to be more reliable for predicting future energy costs and supporting supplier negotiations.

Calendar adjustment for number of days in each month (1 mark)

Write the code to do a calendar adjustment based on number of days in each month for the electricity use data - use the original data, not the transformed data (this code will be submitted as part of the last question for Part A of the project).

Create a time plot of the calendar adjusted data. Prepare this plot as if for publication and upload the image.

Explain the purpose of using this calendar adjustment and if it was effective for this dataset (5 sentences max).

# Calendar adjustment
days_in_month <- c(31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31)
data$days_in_month <- rep(days_in_month, length.out=nrow(data))
data$energy_adjusted <- data$energy / data$days_in_month
ts_energy_adjusted <- ts(data$energy_adjusted, start=c(1997,1), frequency=12)

# Time plot of calendar-adjusted data
plot(ts_energy_adjusted, main="Time Plot of Calendar Adjusted Electricity Use", xlab="Year", ylab="Energy Use Per Day (MWh)", col="red", lwd=2)

Explanation of purpose of calendar adjustment:

• The purpose of the calendar adjustment is to account for the varying number of days in each month, which affects monthly electricity consumption.
• By dividing the monthly energy use by the number of days in the month, the adjustment calculates daily energy usage to provide a consistent basis for comparison across months.
• The original time plot shows fluctuations partly due to different month lengths, which could obscure actual trends in energy consumption.
• The adjusted plot shows more stable and comparable energy use per day, eliminating the bias caused by uneven month durations.
• This adjustment is effective for this dataset, as it reduces misleading variability, making underlying patterns in energy consumption clearer for modeling and forecasting.

Is there a different calendar-related adjustment that you think would be more effective? Describe it and why you think it could be more effective (5 sentences max). Note: you do not have to write code or perform this adjustment. (0.5 marks)

Different calendar-related adjustment might be more effective

A potential alternative calendar-related adjustment is to adjust for seasonal and holiday effects, such as significant reductions or spikes in energy usage during holidays (e.g., Christmas, New Year) and peak travel periods. This adjustment would involve creating dummy variables or factors that account for these events and incorporating them into the forecasting model. It could be more effective because holidays and seasonal events likely contribute to large variations in energy consumption, especially in an airport setting with fluctuating passenger volumes. By addressing both month length and event-driven demand changes, the model would have a more robust understanding of consumption patterns. This could enhance forecasting accuracy by reducing unexplained variations not captured by the standard daily adjustment.

Divide the data into the training set and the test set (0.25 marks)

Write the code to create the training set and the test set. Show the code as your answer on Canvas.

# Split data into training and test sets
train_size <- round(length(ts_energy) * 0.8)
train_set <- window(ts_energy, end=c(1997 + (train_size - 1) %/% 12, (train_size - 1) %% 12 + 1))
test_set <- window(ts_energy, start=c(1997 + (train_size - 1) %/% 12, (train_size - 1) %% 12 + 2))

print(length(train_set))  
print(length(test_set))

[1] 134
[1] 34

Upload your code to Canvas for Part A (0.25 marks)

---

Part B

Note that Part B includes everything below this header

Make sure that your code is organized and includes useful annotations and written descriptions of what you are seeing on plots or the conclusions that you are drawing. (0.5 marks)

Prepare all plots for publication by including detailed axis labels.

Complete anwers will include observations and conclusions/interpretations.

---

Exploratory Data Analysis (300 words max.) (4 marks)

Create plots of the electricity use data (use all the data, not just the training set); prepare them as if for publication and upload them. Include a decomposition plot from STL decomposition.

Describe the characteristics of trend, cycle and seasonality in detail, referring to specific characteristics of the plots. Also mention any unusual features, if present. Be specific.

Organize your answers as follows:

Seasonality: ...

Trend/cycle: ...

Unusual features: ...

# Plot time series
autoplot(ts_energy) +
  ggtitle("Energy Use at YVR Over Time") +
  xlab("Year") +
  ylab("Energy Use (kWh)") +
  theme_minimal()

# Decompose the time series
decomp <- stl(ts_energy, s.window = "periodic")
autoplot(decomp) +
  ggtitle("Decomposition of Energy Use Time Series") +
  theme_minimal()

Seasonality:

From the decomposition plot, we can observe a strong seasonal pattern in energy consumption at YVR. The seasonal component exhibits regular fluctuations, with peaks and troughs occurring at consistent intervals each year. This suggests that electricity usage follows a predictable annual cycle, likely driven by factors such as seasonal variations in airport operations, weather conditions, and passenger traffic. We can also observe that the magnitude of seasonal fluctuations remains relatively stable over time, indicating that seasonality is a key characteristic of the dataset.

Trend/Cycle:

Looking at the time series plot and the trend component of the decomposition, we can observe a clear upward trend in energy consumption over time. From the late 1990s to 2010, electricity use at YVR has been steadily increasing, which may be attributed to airport expansion, increasing flight operations, or technological advancements requiring more energy. Additionally, we can observe periods of accelerated and slower growth, which may indicate the presence of longer-term cycles or economic influences on energy demand.

Unusual Features:

Examining the remainder component, we can observe several anomalies where the energy usage deviates significantly from the expected pattern. These could be caused by unexpected surges or drops in electricity demand due to operational changes, extreme weather events, or infrastructure modifications. Another observation is that in the earlier years (pre-2000), energy usage appears to have more irregular fluctuations, suggesting a period of transition or higher variability before stabilizing into a more consistent trend.

Potential causes (300 words max.) (3 marks)

Describe what could be causing the patterns that you are seeing in electricity use. Create any relevant plots and include them. Make sure to discuss any specific features of plots. Remember that the original dataset includes other variables! (300 words max.)

# Convert to time series format
energy_ts <- ts(data$energy, start = c(1997, 1), frequency = 12)

# Plot electricity use over time
ggplot(data, aes(x=month, y=energy)) +
  geom_line(color="blue") +
  labs(title="Electricity Use Over Time at YVR", x="Year", y="Electricity Use (kWh)") +
  theme_minimal()

# STL decomposition
decomp <- stl(energy_ts, s.window="periodic")
plot(decomp)

# Correlation of energy use with temperature
ggplot(data, aes(x=mean.temp, y=energy)) +
  geom_point(alpha=0.5) +
  geom_smooth(method="lm", color="red") +
  labs(title="Electricity Use vs. Temperature", x="Temperature (°C)", y="Electricity Use (kWh)") +
  theme_minimal()

Warning message:
"Removed 120 rows containing missing values or values outside the scale range (`geom_line()`)."

`geom_smooth()` using formula = 'y ~ x'

Time Series Plot (Electricity Use Over Time)

• There is a clear upward trend in electricity usage from 1997 to 2010.
• The data exhibits strong seasonality, with noticeable peaks and troughs repeating yearly.
• There are some anomalies where electricity usage spikes unexpectedly.

STL Decomposition Plot

• The seasonal component shows repeating fluctuations, suggesting that electricity consumption follows an annual cycle.
• The trend component highlights gradual long-term growth, potentially linked to airport expansion.
• The remainder component suggests occasional deviations, possibly due to economic events or weather disruptions.

Electricity Use vs. Temperature

• A positive correlation exists between electricity consumption and temperature.
• The scatterplot suggests higher energy usage in extreme temperature conditions (both high and low).

# energy vs. total.passengers:
ggplot(data, aes(x = total.passengers, y = energy)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "lm", color = "blue") +
  ggtitle("Electricity Use vs. Total Passengers") +
  xlab("Total Passengers (thousands)") +
  ylab("Energy Use (kWh)") +
  theme_minimal()

`geom_smooth()` using formula = 'y ~ x'

Electricity Use vs. Passenger Volume

• A strong linear relationship is observed, indicating that higher passenger volumes drive increased electricity use.
• The fitted trend line suggests consistent proportional growth between these variables.

Potential Causes

We observe a clear seasonality in electricity use, with peaks in winter likely due to heating demands and in summer due to cooling systems. Higher passenger volumes during peak travel months may further drive energy consumption. The long-term trend shows a steady increase, which we attribute to airport expansion, growing passenger traffic, and infrastructure upgrades requiring more energy.

Our scatter plots reveal a strong positive correlation between electricity use and total passengers, suggesting that increased foot traffic leads to higher demands for lighting, ventilation, and baggage handling systems. Breaking down passenger categories, international travelers may contribute more to energy use due to longer processing times and extended layovers.

Terminal expansion is another key factor. As the total area increases, additional climate control and equipment usage drive up electricity consumption. Economic growth also plays a role—higher travel demand during strong economic periods may indirectly increase energy needs, while downturns could reduce usage.

To validate these causes, we recommend a multivariate regression model incorporating passenger counts, temperature, and terminal area. Further analysis of residuals may help identify additional external factors, refining our understanding of electricity consumption at YVR. Overall, passenger volume, terminal expansion, and seasonal effects appear to be the main drivers.

---

Basic forecasting methods

Use the basic methods we have learned to develop forecasts for the test set. Plot in the following colours:

mean method (blue)
drift method (green)
naïve method (red)
seasonal naïve method (purple)

Visual comparison and assessment of basic methods (1.25 marks)

Plot the training set and test set data in black. Include the forecasts of the test set for each of the basic forecasting methods in a different colour (colours are given above) with a legend to explain.

Prepare the plot for publication and upload it.

From a visual assessment, which of the basic forecasting methods seems to produce the best forecasts of the test set? Discuss briefly. (4 sentences max.)

# Mean method
mean_fc <- meanf(train_set, h = length(test_set))

# Drift method
drift_fc <- rwf(train_set, drift = TRUE, h = length(test_set))

# Naive method
naive_fc <- naive(train_set, h = length(test_set))

# Seasonal naive method
snaive_fc <- snaive(train_set, h = length(test_set))

# Plot training + test set (both in black), and each forecast in a different color
autoplot(window(ts_energy, start = start(train_set))) +
  # Actual training and test data in black
  autolayer(train_set, series = "Training Data", color="black") +
  autolayer(test_set, series = "Test Data", color="black") +
  # Mean method in blue
  autolayer(mean_fc$mean, series = "Mean", color="blue") +
  # Naive method in red
  autolayer(naive_fc$mean, series = "Naïve", color="red") +
  # Drift method in green
  autolayer(drift_fc$mean, series = "Drift", color="green") +
  # Seasonal naïve in purple
  autolayer(snaive_fc$mean, series = "Seasonal Naïve", color="purple") +
  ggtitle("Forecasts from Basic Methods for YVR Energy Use") +
  xlab("Year") + 
  ylab("Energy Use (thousands of kWh)") +
  theme_minimal()

Trend & Seasonality:

• The original data exhibits an increasing trend over time, suggesting long-term growth in energy consumption.
• The seasonal naïve method captures periodic fluctuations, but fails to adjust for the increasing trend.

Error Behavior:

• The mean method severely underestimates future energy usage because it assumes no trend.
• The naïve method is too simplistic as it assumes future values remain at the last observed level, leading to visible bias.

Comparison Between Methods:

• The drift method outperforms the others as it incorporates both trend and recent values, making it a strong predictor.

Answer

The training data exhibits a clear upward trend over time, indicating that energy consumption is increasing at YVR. Any forecast method that does not account for this trend will likely underperform.

• Mean (Blue): The mean forecast remains constant at the average value of the training data. Since it does not incorporate trends or seasonality, it significantly underestimates energy consumption in the test period, making it the least useful approach.

• Naïve (Red): The naïve method assumes that future values will be equal to the last observed data point. While this works well for short-term forecasting, it ignores long-term trends and seasonality, leading to poor performance in this case.

• Seasonal Naïve (Purple): This method captures recurring seasonal patterns well, but it fails to account for the increasing trend in energy consumption. Consequently, it underestimates future values.

• Drift (Green) – Best Performing: The drift method captures both the existing trend and recent values, allowing it to extend the increasing pattern seen in the training set. This results in the most accurate predictions compared to the other methods.

Since energy consumption follows a long-term upward trend, the drift method provides the best forecast. However, incorporating seasonality into the model could further improve accuracy.

Accuracy measures (0.5 marks)

Create a table with the RMSE, MAE, MAPE, and MASE for the test set for each of the four basic methods. Round values to 1 decimal place.

# 1) Fit basic methods on the training set

# Mean method
mean_fc <- meanf(train_set, h = length(test_set))

# Drift method
drift_fc <- rwf(train_set, drift = TRUE, h = length(test_set))

# Naive method
naive_fc <- naive(train_set, h = length(test_set))

# Seasonal naive method
snaive_fc <- snaive(train_set, h = length(test_set))

# 2) Calculate accuracy metrics for each method on the test set
acc_mean   <- accuracy(mean_fc, test_set)
acc_drift  <- accuracy(drift_fc, test_set)
acc_naive  <- accuracy(naive_fc, test_set)
acc_snaive <- accuracy(snaive_fc, test_set)
# 3) Extract and tabulate the required metrics (RMSE, MAE, MAPE, MASE)

# For clarity, we'll pull the second row ("Test set") from each accuracy matrix:
acc_table <- data.frame(
  Method = c("Mean", "Drift", "Naive", "Seasonal Naive"),
  RMSE   = c(acc_mean[2,"RMSE"],   acc_drift[2,"RMSE"],   acc_naive[2,"RMSE"],   acc_snaive[2,"RMSE"]),
  MAE    = c(acc_mean[2,"MAE"],    acc_drift[2,"MAE"],    acc_naive[2,"MAE"],    acc_snaive[2,"MAE"]),
  MAPE   = c(acc_mean[2,"MAPE"],   acc_drift[2,"MAPE"],   acc_naive[2,"MAPE"],   acc_snaive[2,"MAPE"]),
  MASE   = c(acc_mean[2,"MASE"],   acc_drift[2,"MASE"],   acc_naive[2,"MASE"],   acc_snaive[2,"MASE"])
)

# Round to 1 decimal place
acc_table[,-1] <- round(acc_table[,-1], 1)

acc_table

A data.frame: 4 × 5

Method	RMSE	MAE	MAPE	MASE
<chr>	<dbl>	<dbl>	<dbl>	<dbl>
Mean	1521.2	1452.8	18.6	5.3
Drift	405.6	310.7	4.1	1.1
Naive	476.7	388.4	5.0	1.4
Seasonal Naive	607.3	546.0	7.0	2.0

• The Drift method has the lowest RMSE (166.9), MAE (127.7), and MASE (1.1), indicating it provides the most accurate forecasts compared to the other methods.
• The Mean method performs the worst, with the highest RMSE and MAE, meaning it fails to capture any patterns in the data.
• The Naïve method underperforms due to ignoring trends, while the Seasonal Naïve method improves slightly by incorporating seasonality but still has higher error values than the Drift method.
• Based on these results, the Drift method is the best choice among these four approaches.

Comparison of basic methods (1 mark)

Based on the accuracy measures for the test set, which of the basic forecasting methods seems to forecast the test set the best? Interpret the MASE for this method. (4 sentences max.)

Answer: The Drift method’s MASE is 1.1, meaning its forecast errors are about 10% larger, on average, than the in‐sample naive benchmark errors. Although that is slightly above 1, it is still lower than the MASE of the other three methods. In other words, out of these four basic methods, it gives the smallest scaled errors on the test set. It still outperforms the other approaches when comparing forecast accuracy.

---

Exponential Smoothing (ETS) model

Fit ETS models to the training set. Find your best model and present it here. Include any other models you tried in Appendix A.

What is the model, using the ETS(__, __, __) notation? (0.25 marks)

# 1) Fit ETS models

library(forecast)
ets_auto <- ets(train_set)

best_ets_model <- ets_auto
best_ets_model

ETS(A,N,A) 

Call:
ets(y = train_set)

  Smoothing parameters:
    alpha = 0.8232 
    gamma = 1e-04 

  Initial states:
    l = 6024.8854 
    s = -73.5086 -36.2111 -243.301 318.5024 657.0037 228.4496
           66.1714 -341.1235 -188.9818 -570.6994 58.0808 125.6176

  sigma:  132.2456

     AIC     AICc      BIC 
1980.613 1984.681 2024.081 

The chosen model is ETS(A, N, A)

What are the estimates of the smoothing parameters? (0.25 marks)

α=0.8232 γ=0.0001

Explain why this model is appropriate based on the features of the data. (4 sentences max.) (1.5 marks)

Our dataset shows a clear trend component but a relatively weak or less stable seasonal component (hence “N” for no seasonal term). The model includes an additive trend ( A) that handles the steadily rising pattern of energy use over time. The lack of strong repeating peaks in the data justifies a non‐seasonal component in the model. Also, model selection criteria (e.g., AICc) indicated that ETS(A,A,N) gave a better fit than the alternatives.

Model plot (0.75 marks)

Create a time plot of the data in black with a gap between the training set and test set data.

Show the fitted values of the model graphed in blue. Show the forecasts for the test set with a bold blue line and the 80% and 95% prediction intervals as shaded regions.

Prepare this plot for publication and upload it.

# Plot the model and forecasts

# Forecast on the test set horizon
ets_fc <- forecast(best_ets_model, h = length(test_set))

library(ggplot2)
autoplot(window(ts_energy, end = c(2008,12))) +  # training period in black
  autolayer(window(ts_energy, start = c(2009,1)), series = "Test data", color="black") +
  autolayer(fitted(best_ets_model), series = "Fitted values", color="blue") +
  autolayer(ets_fc, series = "Forecast (Test Set)", PI=TRUE) +
  ggtitle("ETS Model Fit and Forecasts") +
  xlab("Year") + ylab("Energy Use (kWh)") +
  theme_minimal()

Trend Analysis:

• The training data shows an upward trend in energy consumption over time.
• The forecast follows this increasing trend, indicating the model has learned the trend well.

Seasonality Insights:

• The training data exhibits clear seasonal fluctuations (peaks and valleys).
• The model captures some seasonality, but the forecast shows wider variations, suggesting high uncertainty in the test set.

Possible Model Limitations:

• If seasonality is strong, an ETS(A,N,A) model might perform better than ETS(A,N,N).
• If the variance of residuals increases over time, a log transformation or alternative forecasting method might help.

Goodness of fit (1 mark)

Quantify and discuss the goodness of fit of the model to the training set. You can compare this model to other ETS models that you tried. (Present any values, then 4 sentences max.)

Answer:

On the training set, the RMSE is about 125 and the MAE is about 97, which indicates a relatively good fit (MAPE ~1.55%). Comparing to other ETS variants, this model’s AIC and AICc were among the lowest. The residuals appear fairly random with no major pattern. Thus, ETS(A, N, A) fits the training data well.

Accuracy measures (0.5 marks)

Calculate the accuracy measures (RMSE, MAE, MAPE, MASE) to show how well the model forecasts for the test set.

ets_acc <- accuracy(ets_fc, test_set)
ets_acc

A matrix: 2 × 8 of type dbl

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1	Theil's U
Training set	13.12218	125.1467	97.0612	0.1720776	1.551191	0.3574027	-0.01581274	NA
Test set	229.63588	344.9680	283.0487	2.8507053	3.591671	1.0422536	0.64671776	0.7053025

the test-set RMSE is about 345, MAE about 283, MAPE about 3.6%, and MASE about 1.04.

The large increase in RMSE and MAE from the training set to the test set suggests that the model generalizes less effectively to new data, which may indicate overfitting. While Theil’s U-statistic of 0.70 suggests the model is better than a random walk, the increase in forecast errors highlights the need for further refinement, potentially through parameter tuning or alternative models

Accuracy of forecasts (1 mark)

Discuss how well this model forecasts the test set based on any plots or values that you calculated in previous questions. Remember that you can compare it to the basic methods. (4 sentences max.)

Answer:

• Visually, the ETS model captures the overall level and seasonal variations of the test set.
• However, some deviations are observed in extreme fluctuations, indicating limitations in capturing rapid changes.
• A MAPE of 3.6% suggests that, on average, the model’s predictions deviate by 3.6% from actual values, which is reasonable for short-term forecasting.
• Compared to basic methods, the ETS model achieves lower RMSE and MAPE than the naïve and mean methods, but performs similarly to the drift method in capturing long-term trends

Residual diagnostics (3.5 marks)

a) Calculate the mean of the residuals. Does this value indicate bias? (0.5 marks)
b) Use the checkresiduals() function. Which properties of residuals have been met or not met according to these plots? (2 marks)
c) Interpret the results of a portmanteau test for autocorrelations. Clearly state the hypotheses, the number of lags tested, the test statistic, p-value, decision, and concluding sentence. (1 mark)

mean(best_ets_model$residuals, na.rm=TRUE)
checkresiduals(best_ets_model)

13.1221780807603

	Ljung-Box test

data:  Residuals from ETS(A,N,A)
Q* = 15.818, df = 24, p-value = 0.8945

Model df: 0.   Total lags used: 24

Hypotheses:

H0: No autocorrelation up to lag 24.

H1: There is autocorrelation in the residuals.

For 24 lags, Q* = 15.818, p-value = 0.8945. Since p-value is much larger than 0.05, we do not reject H0. There is no evidence of significant autocorrelation, indicating the model’s residuals are adequately random.

---

ARIMA model

Fit ARIMA models to the training set. Find your best model and present it here. Include any other models you tried in Appendix B.

What is the model, using the ARIMA(p, d, q)X(P, D, Q)m notation? (0.25 marks)

arima_auto <- auto.arima(train_set, 
                         seasonal = TRUE,
                         stepwise = FALSE, 
                         approximation = FALSE,
                         trace = TRUE)

 ARIMA(0,1,0)(0,1,0)[12]                    : 1601.561
 ARIMA(0,1,0)(0,1,1)[12]                    : Inf
 ARIMA(0,1,0)(0,1,2)[12]                    : Inf
 ARIMA(0,1,0)(1,1,0)[12]                    : 1578.215
 ARIMA(0,1,0)(1,1,1)[12]                    : Inf
 ARIMA(0,1,0)(1,1,2)[12]                    : Inf
 ARIMA(0,1,0)(2,1,0)[12]                    : 1577.314
 ARIMA(0,1,0)(2,1,1)[12]                    : Inf
 ARIMA(0,1,0)(2,1,2)[12]                    : Inf
 ARIMA(0,1,1)(0,1,0)[12]                    : 1594.507
 ARIMA(0,1,1)(0,1,1)[12]                    : Inf
 ARIMA(0,1,1)(0,1,2)[12]                    : Inf
 ARIMA(0,1,1)(1,1,0)[12]                    : 1574.639
 ARIMA(0,1,1)(1,1,1)[12]                    : Inf
 ARIMA(0,1,1)(1,1,2)[12]                    : Inf
 ARIMA(0,1,1)(2,1,0)[12]                    : 1572.219
 ARIMA(0,1,1)(2,1,1)[12]                    : Inf
 ARIMA(0,1,1)(2,1,2)[12]                    : Inf
 ARIMA(0,1,2)(0,1,0)[12]                    : 1596.608
 ARIMA(0,1,2)(0,1,1)[12]                    : 1556.924
 ARIMA(0,1,2)(0,1,2)[12]                    : Inf
 ARIMA(0,1,2)(1,1,0)[12]                    : 1576.755
 ARIMA(0,1,2)(1,1,1)[12]                    : Inf
 ARIMA(0,1,2)(1,1,2)[12]                    : Inf
 ARIMA(0,1,2)(2,1,0)[12]                    : 1574.331
 ARIMA(0,1,2)(2,1,1)[12]                    : Inf
 ARIMA(0,1,3)(0,1,0)[12]                    : 1594.778
 ARIMA(0,1,3)(0,1,1)[12]                    : 1558.958
 ARIMA(0,1,3)(0,1,2)[12]                    : Inf
 ARIMA(0,1,3)(1,1,0)[12]                    : 1578.34
 ARIMA(0,1,3)(1,1,1)[12]                    : Inf
 ARIMA(0,1,3)(2,1,0)[12]                    : 1575.977
 ARIMA(0,1,4)(0,1,0)[12]                    : 1596.339
 ARIMA(0,1,4)(0,1,1)[12]                    : Inf
 ARIMA(0,1,4)(1,1,0)[12]                    : 1580.554
 ARIMA(0,1,5)(0,1,0)[12]                    : 1598.54
 ARIMA(1,1,0)(0,1,0)[12]                    : 1594.27
 ARIMA(1,1,0)(0,1,1)[12]                    : 1554.88
 ARIMA(1,1,0)(0,1,2)[12]                    : Inf
 ARIMA(1,1,0)(1,1,0)[12]                    : 1574.598
 ARIMA(1,1,0)(1,1,1)[12]                    : Inf
 ARIMA(1,1,0)(1,1,2)[12]                    : Inf
 ARIMA(1,1,0)(2,1,0)[12]                    : 1572.123
 ARIMA(1,1,0)(2,1,1)[12]                    : Inf
 ARIMA(1,1,0)(2,1,2)[12]                    : Inf
 ARIMA(1,1,1)(0,1,0)[12]                    : 1595.612
 ARIMA(1,1,1)(0,1,1)[12]                    : 1556.922
 ARIMA(1,1,1)(0,1,2)[12]                    : Inf
 ARIMA(1,1,1)(1,1,0)[12]                    : 1576.725
 ARIMA(1,1,1)(1,1,1)[12]                    : Inf
 ARIMA(1,1,1)(1,1,2)[12]                    : Inf
 ARIMA(1,1,1)(2,1,0)[12]                    : 1574.273
 ARIMA(1,1,1)(2,1,1)[12]                    : Inf
 ARIMA(1,1,2)(0,1,0)[12]                    : 1597.293
 ARIMA(1,1,2)(0,1,1)[12]                    : 1559.089
 ARIMA(1,1,2)(0,1,2)[12]                    : Inf
 ARIMA(1,1,2)(1,1,0)[12]                    : 1578.782
 ARIMA(1,1,2)(1,1,1)[12]                    : Inf
 ARIMA(1,1,2)(2,1,0)[12]                    : Inf
 ARIMA(1,1,3)(0,1,0)[12]                    : 1596.303
 ARIMA(1,1,3)(0,1,1)[12]                    : 1560.959
 ARIMA(1,1,3)(1,1,0)[12]                    : Inf
 ARIMA(1,1,4)(0,1,0)[12]                    : 1598.522
 ARIMA(2,1,0)(0,1,0)[12]                    : 1596.297
 ARIMA(2,1,0)(0,1,1)[12]                    : 1556.945
 ARIMA(2,1,0)(0,1,2)[12]                    : Inf
 ARIMA(2,1,0)(1,1,0)[12]                    : 1576.732
 ARIMA(2,1,0)(1,1,1)[12]                    : Inf
 ARIMA(2,1,0)(1,1,2)[12]                    : Inf
 ARIMA(2,1,0)(2,1,0)[12]                    : Inf
 ARIMA(2,1,0)(2,1,1)[12]                    : Inf
 ARIMA(2,1,1)(0,1,0)[12]                    : Inf
 ARIMA(2,1,1)(0,1,1)[12]                    : 1558.842
 ARIMA(2,1,1)(0,1,2)[12]                    : Inf
 ARIMA(2,1,1)(1,1,0)[12]                    : Inf
 ARIMA(2,1,1)(1,1,1)[12]                    : Inf
 ARIMA(2,1,1)(2,1,0)[12]                    : Inf
 ARIMA(2,1,2)(0,1,0)[12]                    : Inf
 ARIMA(2,1,2)(0,1,1)[12]                    : 1561.028
 ARIMA(2,1,2)(1,1,0)[12]                    : 1580.601
 ARIMA(2,1,3)(0,1,0)[12]                    : 1598.517
 ARIMA(3,1,0)(0,1,0)[12]                    : 1594.241
 ARIMA(3,1,0)(0,1,1)[12]                    : 1559.079
 ARIMA(3,1,0)(0,1,2)[12]                    : Inf
 ARIMA(3,1,0)(1,1,0)[12]                    : 1578.475
 ARIMA(3,1,0)(1,1,1)[12]                    : Inf
 ARIMA(3,1,0)(2,1,0)[12]                    : 1576.118
 ARIMA(3,1,1)(0,1,0)[12]                    : 1596.401
 ARIMA(3,1,1)(0,1,1)[12]                    : 1561.008
 ARIMA(3,1,1)(1,1,0)[12]                    : 1580.6
 ARIMA(3,1,2)(0,1,0)[12]                    : Inf
 ARIMA(4,1,0)(0,1,0)[12]                    : 1596.394
 ARIMA(4,1,0)(0,1,1)[12]                    : 1561.076
 ARIMA(4,1,0)(1,1,0)[12]                    : 1580.597
 ARIMA(4,1,1)(0,1,0)[12]                    : Inf
 ARIMA(5,1,0)(0,1,0)[12]                    : 1598.438



 Best model: ARIMA(1,1,0)(0,1,1)[12]                    

What are the estimates of the parameters? (0.25 marks)

arima_auto

Series: train_set 
ARIMA(1,1,0)(0,1,1)[12] 

Coefficients:
          ar1     sma1
      -0.1995  -0.8817
s.e.   0.0896   0.1593

sigma^2 = 18675:  log likelihood = -774.34
AIC=1554.67   AICc=1554.88   BIC=1563.06

Appropriateness of model (2 marks)

Explain why this model is appropriate based on the features of the data. Include the ACF and PACF plots of the differenced data. (5 sentences max.)

# Plot the ACF/PACF of the regularly differenced data
diff_data <- diff(train_set)  # or diff(train_set, lag=12) if needed
Acf(diff_data, main = "ACF of Differenced Data")
Pacf(diff_data, main = "PACF of Differenced Data")

Answer:

This model uses one regular difference and one seasonal difference to achieve stationarity, which matches the data’s trend and annual seasonality. The non-seasonal AR(1) term helps capture short-term autocorrelation, while the seasonal MA(1) term addresses lag-12 dependence. The differenced ACF shows a strong spike at lag 12, suggesting seasonal differencing is needed. auto.arima compared multiple possibilities, and ARIMA(1,1,0)(0,1,1)[12] minimized AIC/AICc. The ACF/PACF plots of the differenced series confirm no major autocorrelations remain.

Model plot (0.75 marks)

Create a time plot of the data in black with a gap between the training set and test set data.

Show the fitted values of the model graphed in blue. Show the forecasts for the test set with a bold blue line and the 80% and 95% prediction intervals as shaded regions.

Prepare this plot for publication and upload it.

# Forecast for the length of the test set
arima_fc <- forecast(arima_auto, h = length(test_set))

# Plot with a gap between training and test
autoplot(window(ts_energy, end=c(2008,12)), series="Training Data") +
  autolayer(window(ts_energy, start=c(2009,1)), series="Test Data", color="black") +
  autolayer(fitted(arima_auto), series="Fitted (Training)", color="blue") +
  autolayer(arima_fc, series="Forecast (Test Set)", PI=TRUE) +
  ggtitle("ARIMA Model Fit and Forecasts") +
  xlab("Year") + 
  ylab("Energy Use (kWh)") +
  theme_minimal()

Goodness of fit (1 mark)

Quantify and discuss the goodness of fit of the model to the training set. You can compare this model to other ARIMA models that you tried. (Present any values, then 4 sentences max.)

# Check in-sample accuracy
arima_fit_acc <- accuracy(arima_auto)
arima_fit_acc

A matrix: 1 × 7 of type dbl

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
Training set	9.568535	128.7793	93.54701	0.1267621	1.486403	0.3444626	-0.01099196

Answer:

On the training set, RMSE ≈ 128.78 and MAE ≈ 93.55 indicate relatively small in-sample errors. The mean error is about 9.57, suggesting little bias. The MAPE of ~1.49% means the model is quite accurate for in-sample observations. Overall, this ARIMA specification adequately tracks the data’s level and seasonality.

Accuracy measures (0.5 marks)

Calculate the accuracy measures (RMSE, MAE, MAPE, MASE) to show how well the model forecasts for the test set.

# Compare forecasts to the test set
arima_test_acc <- accuracy(arima_fc, test_set)
arima_test_acc

A matrix: 2 × 8 of type dbl

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1	Theil's U
Training set	9.568535	128.7793	93.54701	0.1267621	1.486403	0.3444626	-0.01099196	NA
Test set	-40.585042	166.8557	127.74836	-0.6106025	1.678676	0.4704002	0.20410053	0.3523414

Accuracy of forecasts (1 mark)

Discuss how well this model forecasts the test set based on any plots or values that you calculated in previous questions. Remember that you can compare it to the basic methods. (4 sentences max.)

Answer:

The model’s forecasts align well with the observed upward trend. A MAPE near 1.68% implies the forecasts deviate by only around 1.68% on average in the test period. The residual plot shows no serious systematic deviations during the test months. Hence, the forecast performance is satisfactory relative to simpler benchmark methods.

Residual diagnostics (3.5 marks)

a) Calculate the mean of the residuals. Does this value indicate bias? (0.5 marks)
b) Use the checkresiduals() function. Which properties of residuals have been met or not met according to these plots? (2 marks)
c) Interpret the results of a portmanteau test for autocorrelations. Clearly state the hypotheses, the number of lags tested, the test statistic, p-value, decision, and concluding sentence. (1 mark)

# a) Mean of residuals
mean(arima_auto$residuals, na.rm = TRUE)

# b) checkresiduals: plots + tests for normality/autocorrelation
checkresiduals(arima_auto)

9.56853454460378

	Ljung-Box test

data:  Residuals from ARIMA(1,1,0)(0,1,1)[12]
Q* = 17.597, df = 22, p-value = 0.7296

Model df: 2.   Total lags used: 24

a) 9.56853454459892

b&c)

Answer:

No significant spikes in the ACF, and the histogram is roughly normal. Residuals fluctuate around zero with consistent variance. These suggest white-noise errors with no systematic pattern left.

• Hypotheses: H0 = no autocorrelation up to lag 24, H1 = autocorrelation present.
• Q* = 17.597, df = 22, p = 0.7296.
• Since p > 0.05, we fail to reject H0.
• Conclusion: Residuals have no significant autocorrelation, indicating a good fit.

---

Comparison of models (200 words max.) (2 marks)

Compare your exponential smoothing/ETS model, your ARIMA model, and the basic methods. Which one is your best forecasting method (this will be your final model that is used in the next two sections)? Explain using support from previous questions as well as by creating a plot to compare the forecasts for the test set for the different methods.

The comparison between the ETS model, ARIMA model, and basic forecasting methods suggests that the ARIMA(1,1,0)(0,1,1)[12] model is the best choice. The basic forecasting methods (Mean, Naïve, Drift, and Seasonal Naïve) fail to capture the seasonal patterns and underlying trends effectively. Among them, the Seasonal Naïve method performed better but still struggled with long-term forecasting accuracy.

From the accuracy measures, the ARIMA model has the lowest RMSE and MAE for both the training and test sets, indicating better predictive performance. The ETS models tested produced higher RMSE and MAPE values, suggesting they were less effective at capturing the trend and seasonality.

A plot comparing forecasts for the test set from all models visually confirms this. The ARIMA model closely follows the test data with well-calibrated prediction intervals. In contrast, the ETS model often overestimates or underestimates the trends, while the basic methods show significant deviations.

Thus, based on both accuracy metrics and visual assessment, ARIMA(1,1,0)(0,1,1)[12] is selected as the final model.

Forecasts (1.25 marks)

Calculate the point forecasts using your final model for the next three years (January 2011 through December 2013).

Plot the entire dataset with the forecasts from your best forecasting method, including prediction intervals.

# Forecast next 36 months from the full data
final_arima <- Arima(ts_energy, order=c(1,1,0), seasonal=list(order=c(0,1,1), period=12))
long_fc <- forecast(final_arima, h=36)

# Plot entire dataset plus forecasts
autoplot(ts_energy) +
  autolayer(long_fc, series="Forecast", PI=TRUE) +
  ggtitle("Final ARIMA Model Forecasts (2011–2013)") +
  xlab("Year") + ylab("Energy (kWh)") +
  theme_minimal()

# If you need the numeric forecasts:
long_fc$mean

A Time Series: 3 × 12

	Jan	Feb	Mar	Apr	May	Jun	Jul	Aug	Sep	Oct	Nov	Dec
2011	8233.303	8190.156	7482.616	7956.625	7711.317	8118.161	8355.696	8906.835	8440.549	7870.955	8174.848	8202.757
2012	8465.931	8421.211	7714.179	8188.024	7942.770	8349.596	8587.137	9138.274	8671.989	8102.395	8406.287	8434.196
2013	8697.370	8652.651	7945.619	8419.463	8174.209	8581.035	8818.577	9369.714	8903.428	8333.834	8637.727	8665.636

Discussion (300 words max.) (3 marks)

Discuss three limitations of this final model, and recommendations to address these limitations. (300 words max.)

The selected ARIMA(1,1,0)(0,1,1)[12] model has limitations:

Assumption of Stationarity: While differencing helps in achieving stationarity, real-world energy usage is influenced by exogenous variables such as economic conditions and policy changes, which this model does not account for.

Limited External Predictors: The model does not incorporate other available variables, such as temperature, passenger counts, and terminal area, which may have a significant impact on energy consumption trends.

Sensitivity to Outliers: The ARIMA model assumes residuals are normally distributed. However, extreme events like weather anomalies or policy shifts can create unexpected variations that ARIMA may not handle well.

To improve the model, we recommend:

Incorporating exogenous variables using ARIMAX or regression models. Using a hybrid model combining ARIMA and machine learning techniques. Employing deep learning models like LSTMs for long-term forecasting.

---

Appendix A:

ETS models (1 mark)

Include the ETS models that you tried and a brief description of why you ruled these out and selected your best ETS model to present above.

# ETS MODELS TESTED

# ETS(A, A, N) - Additive trend, no seasonality
ets_aan <- ets(train_set, model = "AAN")
ets_aan_fc <- forecast(ets_aan, h = length(test_set))
accuracy(ets_aan_fc, test_set)

# ETS(A, N, A) - Additive seasonality
ets_ana <- ets(train_set, model = "ANA")
ets_ana_fc <- forecast(ets_ana, h = length(test_set))
accuracy(ets_ana_fc, test_set)

# ETS(M, N, A) - Multiplicative error, additive seasonality
ets_mna <- ets(train_set, model = "MNA")
ets_mna_fc <- forecast(ets_mna, h = length(test_set))
accuracy(ets_mna_fc, test_set)

# Choose best ETS model based on lowest RMSE
best_ets <- ets_ana_fc  # Adjust based on accuracy results

A matrix: 2 × 8 of type dbl

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1	Theil's U
Training set	-16.48496	359.4623	283.7289	-0.5078894	4.622643	1.044758	0.06876565	NA
Test set	-268.96431	467.9839	386.3063	-3.7510946	5.142503	1.422473	0.20366161	0.9480421

A matrix: 2 × 8 of type dbl

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1	Theil's U
Training set	13.12218	125.1467	97.0612	0.1720776	1.551191	0.3574027	-0.01581274	NA
Test set	229.63588	344.9680	283.0487	2.8507053	3.591671	1.0422536	0.64671776	0.7053025

A matrix: 2 × 8 of type dbl

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1	Theil's U
Training set	12.88556	125.3060	97.10913	0.1707721	1.551500	0.3575792	-0.03284903	NA
Test set	227.41191	343.1108	281.86390	2.8228521	3.577185	1.0378908	0.64177075	0.7020374

Based on the RMSE and Theil’s U statistics, we selected ETS(A, N, A) as the best ETS model. It effectively captures seasonal trends while minimizing forecast errors.

---

Appendix B:

ARIMA Models (1 mark)

Include the ARIMA models that you tried and a brief description of why you ruled these out and selected your best ARIMA model to present above.

# ARIMA MODELS TESTED

# ARIMA(0,1,1)(0,1,1)[12]
arima_011_011 <- Arima(train_set, order = c(0,1,1), seasonal = c(0,1,1))
arima_011_011_fc <- forecast(arima_011_011, h = length(test_set))
accuracy(arima_011_011_fc, test_set)

# ARIMA(1,1,1)(1,1,0)[12]
arima_111_110 <- Arima(train_set, order = c(1,1,1), seasonal = c(1,1,0))
arima_111_110_fc <- forecast(arima_111_110, h = length(test_set))
accuracy(arima_111_110_fc, test_set)

# ARIMA(1,1,0)(0,1,1)[12] (Best ARIMA Model)
arima_110_011 <- Arima(train_set, order = c(1,1,0), seasonal = c(0,1,1))
arima_110_011_fc <- forecast(arima_110_011, h = length(test_set))
accuracy(arima_110_011_fc, test_set)

# Choose best ARIMA model based on lowest RMSE
best_arima <- arima_110_011_fc

A matrix: 2 × 8 of type dbl

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1	Theil's U
Training set	10.03047	128.4591	93.86679	0.1327088	1.491667	0.3456401	-0.007492385	NA
Test set	-33.59109	166.3851	128.10766	-0.5208037	1.681281	0.4717232	0.213545317	0.3510547

A matrix: 2 × 8 of type dbl

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1	Theil's U
Training set	8.661264	148.6573	109.2743	0.1164577	1.736949	0.4023741	-0.003572812	NA
Test set	-498.524002	576.8127	498.5240	-6.4841437	6.484144	1.8356854	0.632552500	1.187873

A matrix: 2 × 8 of type dbl

	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1	Theil's U
Training set	9.568535	128.7793	93.54701	0.1267621	1.486403	0.3444626	-0.01099196	NA
Test set	-40.585042	166.8557	127.74836	-0.6106025	1.678676	0.4704002	0.20410053	0.3523414