PA Homework 1
Cameron Turner
2024-03-24
Massachusetts Gas Prices
Setup
# Clear the workspace
rm(list = ls()) # Clear environment
gc() # Clear unused memory## used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
## Ncells 521515 27.9 1155397 61.8 NA 669291 35.8
## Vcells 961152 7.4 8388608 64.0 16384 1840035 14.1 cat("\f") # Clear the console dev.off # Clear the charts## function (which = dev.cur())
## {
## if (which == 1)
## stop("cannot shut down device 1 (the null device)")
## .External(C_devoff, as.integer(which))
## dev.cur()
## }
## <bytecode: 0x7f9f0cbcff20>
## <environment: namespace:grDevices>#make vector containing the packages
packages <- c("ggplot2",
"knitr",
"fpp3",
"forecast",
"forecastHybrid",
"dplyr")
#Install packages
for (i in 1:length(packages)) {# Iterate through each package in the 'packages' vector.
if (!packages[i] %in% rownames(installed.packages())) {#check if package is installed
install.packages(packages[i], #install package if it isn't already
repos = "http://cran.rstudio.com/",
dependencies = TRUE)}
library(packages[i], character.only = TRUE)} # Load the package into the R session.## ── Attaching packages ────────────────────────────────────────────── fpp3 0.5 ──## ✔ tibble 3.2.1 ✔ tsibbledata 0.4.1
## ✔ dplyr 1.1.1 ✔ feasts 0.3.2
## ✔ tidyr 1.3.0 ✔ fable 0.3.3
## ✔ lubridate 1.9.2 ✔ fabletools 0.4.1
## ✔ tsibble 1.1.4## ── Conflicts ───────────────────────────────────────────────── fpp3_conflicts ──
## ✖ lubridate::date() masks base::date()
## ✖ dplyr::filter() masks stats::filter()
## ✖ tsibble::intersect() masks base::intersect()
## ✖ tsibble::interval() masks lubridate::interval()
## ✖ dplyr::lag() masks stats::lag()
## ✖ tsibble::setdiff() masks base::setdiff()
## ✖ tsibble::union() masks base::union()## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo## Loading required package: thief# Remove unnecessary values from environment
rm(i,packages)
### Importing data
gas <- read.csv("/Users/cam/Desktop/PREDICTIVE ANALYTICS/Homework/MA_Gas.csv")
# Convert to date format
gas$Date <- as.Date(gas$Date)
# Rename column of focus to be simpler
names(gas)[2] <- "Price"
# Filter gas to last 5 years
gas <- subset(gas, Date >= max(gas$Date) - years(5))
# Exclude year 5 for the sake of testing data
gas_ts <- subset(gas, Date <= max(gas$Date) - years(1))
# Convert to time series
gas_ts <- ts(gas_ts$Price,
start = c(2018,12), #start year, start month
frequency = 12)
# Make separate time series for testing data
gas_ts_test <- subset(gas, Date > max(gas$Date) - years(1))
gas_ts_test <- ts(gas_ts_test$Price,
start = c(2023,1), #start year, start month
frequency = 12)Observing the data
Source: eia.gov
The specific metric of the price data is “Massachusetts Price of Natural Gas Delivered to Residential Consumers (Dollars per Thousand Cubic Feet).”
# Show summary stats
summary(gas)## Date Price
## Min. :2018-12-01 Min. :13.24
## 1st Qu.:2020-03-01 1st Qu.:14.97
## Median :2021-06-01 Median :16.44
## Mean :2021-06-01 Mean :17.58
## 3rd Qu.:2022-09-01 3rd Qu.:20.04
## Max. :2023-12-01 Max. :25.26#Show head of table
head(gas)## Date Price
## 360 2018-12-01 15.21
## 361 2019-01-01 14.97
## 362 2019-02-01 14.58
## 363 2019-03-01 14.68
## 364 2019-04-01 15.16
## 365 2019-05-01 14.70# Show plot of data
ggplot(gas, aes(x = Date, y = Price)) +
geom_line() +
labs(x = "Date",
y = "$ Per Thousand Cubic Feet",
title = "MA Gas Prices (Past 5 Years)") +
scale_x_date(date_labels = "%Y", date_breaks = "1 years")
Modeling
Model 1: Naive
\[ \hat{Price_{T+h}}=Price_{T} \]
This model simply forecasts that the future value will simply equal the last recorded value. In this case, forecasted values equal gas price for December 2022 ($22.64 per thousand cubic feet). This proved to make the model the most accurate relatively speaking. However, I would not prefer it to other models. The reason for its accuracy in terms of MSE is that the residuals are lower after the unexpectedly steep drop in early 2023. It does not appear to actually fit the data well despite its accuracy in predicting the price in December 2024 (end of the forecasted period).
# Create naive model
naive_mod <- naive(gas_ts)
autoplot(naive_mod) +
autolayer(gas_ts_test,
series = "True Values") +
labs(x = "Date",
y = "$ Per Thousand Cubic Feet",
title = "NAIVE Forecast")
# Evaluate accuracy
# Make a df
naive_data <- data.frame(gas_ts_test)
naive_forecast<- naive_mod$mean
naive_data$naive_forecast <- naive_forecast[1]
# residuals
naive_data$residuals <- naive_data$gas_ts_test - naive_data$naive_forecast
# MSE
MSE_naive <- mean(naive_data$residuals^2)
print(paste("MSE: ",MSE_naive))## [1] "MSE: 10.8997166666667"Model 2: Seasonal Naive (SNAIVE)
\[ \hat{y}_{T+h|T}=Y_{T+h-m(k+1)} \]
Seasonal naive simply predicts that the data will replicate its movement from the last period. So in this case, the predicted values for each month are the same as they were a year before. While there appears to be some seasonality to the data, the price dropped steeply in the first half of 2024 and significantly threw off the model’s accuracy
# Create SNAIVE model
snaive_mod <- snaive(gas_ts)
# forecast SNAIVE model
snaive_forecast <- forecast(snaive_mod, h = 12)$mean
# Display forecast
autoplot(snaive_mod) +
autolayer(gas_ts_test,
series = "true_values") +
labs(x = "Date",
y = "$ Per Thousand Cubic Feet",
title = "SNAIVE Forecast")
# Evaluate accuracy
# Make combined df
snaive_data <- data.frame(cbind(gas_ts_test,snaive_forecast))
snaive_data <- na.omit(snaive_data)
#Residuals
snaive_data$residuals <- snaive_data$gas_ts_test - snaive_data$snaive_forecast
#MSE
MSE_snaive <- mean(snaive_data$residuals^2)
print(paste("MSE: ",MSE_snaive))## [1] "MSE: 16.0562083333333"Model 3: ETS
I allowed R to determine which ETS model would be ideal for the data. It returned an MNA model which was not especially accurate. This flaw could mostly be attributed to the steep drop in price in early 2023. As was seen earlier in the SNAIVE model, the seasonal component did appear to be fairly accurate in predicting the direction of changes in prices. Adding a trend component (not seen in this document) hurt the accuracy of predictions for 2023. This is due to a noticeable upward trend in the data and the data, as mentioned earlier, made a steep drop. By omitting the trend component and focusing on seasonality, the data does appear to fit better. After the price fall, the data seems to almost perfectly outline the lower bound of the forecasted values. Accounting for trend in an MAA model may be good for longer term predictions, but the drop makes this model better for 2023.
# Create ets model
ets_mod <- ets(gas_ts)
# forecast ets model
ets_forecast <- forecast(ets_mod, h = 12)
autoplot(ets_mod)
# Display ETS model
autoplot(gas_ts)+
autolayer(ets_forecast)+
labs(x = "Date",
y = "$ Per Thousand Cubic Feet",
title = "ETS Forecast") +
autolayer(gas_ts_test,
series = "True Values")
# Evaluate accuracy
# Make combined df
ets_data <- ets_forecast$mean
ets_data <- data.frame(gas_ts_test,ets_data)
# Residuals
ets_data$residuals <- ets_data$gas_ts_test - ets_data$ets_data
#MSE
MSE_ets <- mean(ets_data$residuals^2)
print(paste("MSE: ",MSE_ets))## [1] "MSE: 11.38790221951"Model 4: Drift
An unexpected challenge arose where the “drift” operator was not working. Thus, I performed the function manually using a coefficient from a linear regression. This appears to be more accurate overall than seasonal data. While this model may not be the best with predicting monthly data, it seems to be a much better long term predictor as we can see prices spiking towards the end of 2023. Considering that prices increase with inflation as well, the prices should surely make their way upwards further on. Additionally, it appears to be more accurate than SNAIVE, which focuses solely on the seasonal component.
# Create drift model
#drift_mod <- drift(gas_ts) ## Error in drift(gas_ts) : could not find function "drift"
# Do it manually
# Create a time index
time_index <- 1:length(gas_ts)
# Fit a linear regression model
drift_reg <- lm(gas_ts ~ time_index)
# Extract the drift component from the model
drift_mod <- coef(drift_reg)[2]
# Generate forecasts
start_year <- 2023
start_month <- 1
forecast_length <- 12
# Create a sequence of months starting from January 2023
forecast_dates <- seq(as.Date(paste(start_year, start_month, "01", sep = "-")), by = "month", length.out = forecast_length)
# Generate time index for forecasting
time_index_forecast <- (length(gas_ts) + 1):(length(gas_ts) + forecast_length)
# Predict the monthly values for the next 12 months
forecast_values <- predict(drift_reg, newdata = data.frame(time_index = time_index_forecast))
# Create a time series object for the forecast
drift_forecast <- ts(forecast_values, start = c(start_year, start_month), frequency = 12)
# Display forecast
autoplot(gas_ts) +
autolayer(gas_ts_test,
series = "True Values")+
autolayer(drift_forecast,
series = "Forecast Values")
## Evaluate accuracy
drift_data <- data.frame(gas_ts_test,drift_forecast)
#residuals
drift_data$residuals <- drift_data$gas_ts_test - drift_data$drift_forecast
# MSE
MSE_drift <- mean(drift_data$residuals^2)
print(paste("MSE: ",MSE_drift))## [1] "MSE: 11.4389157421436"