Introduction

Gas prices impact our daily lives significantly. Whether it’s driving to work, transporting goods, or planning a road trip, the cost of gas can affect our decisions and budgets. With this in mind, I set out to create a simple, yet thorough model that could predict future gas prices. In this blog post, we’ll take you through the project, explaining the process in simple terms and sharing our findings.

The primary engine behind this project is Fable, an easy-to-use package for building and evaluating forecasting models. It includes versatile models like ARIMA and ETS, integrates well with popular packages like tidyverse, and offers great visualization tools for forecasts. It also allows you to compare different models to determine the most accurate one, making forecasting more accessible and effective for time series data.

The Goal

My mission was to forecast gas prices for the upcoming months. By predicting these prices, I hoped to help better project for changes based on the national price trends. To achieve this, I utilized various data analysis tools in R, a popular programming language for statistics and data science.

Gas prices are affected by a number of external factors, making a prediction using multiple variables more reliable for month to month accuracy. However, we can still look at trends and historical data to make a guess at general, long-term trends.

https://rpubs.com/bseko/gas_forecasting

While these trends may appear simplistic. With as resource as volatile as the gasoline index, which can spike with a single hurricane, our goal is to get a general idea of where prices are heading and not read too much into the month-to-month values.

So how did I create this forecast?

Gathering the Data

First, I needed data. We collected historical gas prices from my curated database using the EIA api to grab histrocal data. You can learn more about that process here:

Focusing on monthly averages, this data provided me with a solid foundation, capturing how gas prices have fluctuated over time. Imagine looking at a detailed record of gas prices at your local station for the past several years – that’s the kind of information I worked with.

You can see a more comprehensive review of this data here: https://rpubs.com/bseko/gas_price_tracker

Defining the Time Frame and Splitting the Data

To make accurate predictions, I established specific time frames for training and testing, then used the entire timeline on the best model. I defined separate data sets for training (to build the model) and testing (to evaluate the model). It’s like studying past trends to predict future outcomes, similar to how meteorologists use historical weather data to forecast tomorrow’s weather.

The training set included older data to teach our model about past trends, while the testing set included more recent data to see how well the model could predict prices it hadn’t seen before. It’s like using old exams to prepare for a final test – the better you know the material, the better you perform on new questions. For this project, I train the model on all data from 1994 through June 2023, saving the last year’s worth of data for testing.

Given the extreme changes in prices over the past few years I don’t expect the model to mimic actual data. Rather, can it following the right seasonal patterns, and the current trend (which has been downward). I also care if it generalizes well, which means it will interpret new data just as well as the training and current data it is given. Good data models predict the data you have, and adapt for new data without needing to be constantly retrained.

Exploring Different Models

I explored several forecasting models each with its own strengths. These included:

1. Mean Model: A simple approach that uses the average of past prices.
2. Year-Over-Year Model: This model assumes that prices follow a yearly pattern.
3. Exponential Smoothing (ETS) Model: This one adjusts for trends and seasonal patterns in the data.
4. ARIMA Model: A more complex model that can capture various trends and patterns.
5. VAR Model: This model looks at multiple variables and their interactions over time.

Tuning the Models

To find the best model, I performed a grid search – a systematic way of testing different combinations of parameters. Imagine trying out different recipes to see which one makes the best cake. I evaluated each model based on how accurately it predicted gas prices, using measures like Mean Absolute Error (MAE) and Mean Absolute Percentage Error (MAPE).

Grid searches throw everything at the testing data, and see’s what is best. That doesn’t mean this is done without though though. The models and hyperparameters (the parameters the model doesn’t see) chosen, are those that generally do well with additive trends like the gas prices.

Choosing the Best Model

After testing, I selected the model that performed the best. This model could predict gas prices with the least error, making it the most reliable option. While minimal error isn’t a cure-all for forecasting. In this case, the trends are believable so the model with lowest error is a good choice. Sometimes, the input data needs to be adjusted to get better results, but that’s beyond the scope of this project.

Selecting the best model was like finding the recipe that consistently produced the tastiest cake. In this case we have an ARIMA model, that overall, gave us the lowest MAPE and intuitive trends. I’ll use this to train data for the upcoming months, and track the results.

Making Predictions

With my chosen model, I forecasted gas prices for the next few months. This step involved feeding the model my training data and letting it generate future prices. The result was a series of predictions that could help inform about potential price changes.

Saving and Sharing the Results

Finally, I saved my model and the predictions. I also stored these forecasts in a database, making it easily accessible for further analysis or reporting. For this project, I’ll come back to this output over several months before retraining the model, or even re-running the grid search code.

It’s always a good idea to save your model as well, you can train on new data, and continue to use the grid search results until the model is no longer generalizing well. Grid searches can take a lot of time, in this case ~40 minutes. So it makes sense to make the most out of your results.

The ARIMA model forecast output showcases both historical data and future predictions. The chart includes three key elements: the actual gas prices (in red), the forecasted values (in green), and the training forecast (in blue). The historical data shows fluctuating prices from 2021 to 2024, with notable peaks and troughs. The model’s forecast, starting in mid-2024, suggests a slight dip followed by a modest upward trend in gas prices through 2026. The training forecast closely follows the actual prices, indicating a well-fitted model that captures the underlying patterns in the historical data. This forecast provides valuable insights for stakeholders looking to anticipate future gas price trends and make informed decisions based on these predictions.

My journey to forecast gas prices was rewarding and intersting. By leveraging advanced data analysis tools in R, I was able to create a reliable model that provides valuable insights into future gas prices. This project highlights the power of data science and its potential to make a real-world impact.

To get a full breakdown of the code, keep reading, you can also see all the details in the gas_forecasting git repo.

Code Explaination

I established an engine, or a dedicated script of functions to run my code. This approach streamlines the forecasting process. For most projects, it is beneficial to have dedicated functions. These functions are portable to new projects and act like small, lightweight packages.

1) Creating Global Parameters

The create_global_start_end function sets up the global parameters for the time series, including the start and end dates for forecasting and testing periods. This setup is crucial for ensuring consistency across different functions.

create_global_start_end <- function(fcst_start, fcst_end, fcst_test_start, fcst_test_end) {
  global_fcst_start <<- as.Date(fcst_start)
  global_fcst_end <<- as.Date(fcst_end)
  global_fcst_test_start <<- as.Date(fcst_test_start)
  global_fcst_test_end <<- as.Date(fcst_test_end)
}

2) Calculating Interval in Months

The create_interval_months function calculates the interval in months between the provided dates. This is useful for understanding the time span of the training and forecasting periods and setting defaults so you don’t need to mess with it in multiple files.

create_interval_months <- function(){
  train_interval_months <<- floor((as.numeric(as.yearmon(global_fcst_test_end) - as.yearmon(global_fcst_test_start)) * 12) + 1)
  forecast_interval_months <<- 6
}

3) Splitting Data into Training and Testing Sets

The test_train_split function splits the input object, which in this case is a tsibble, into training, testing, and forecasting sets based on the global parameters. This ensures that the models are trained and tested on the correct portions of the data.

test_train_split <- function(input_tsibble) {
  train_start <- format(global_fcst_start, "%Y %b")
  train_end <- format(global_fcst_test_start %m-% months(1), "%Y %b")
  test_start <- format(global_fcst_test_start, "%Y %b")
  test_end <- format(global_fcst_test_end, "%Y %b")
  train <- input_tsibble %>% filter_index(train_start ~ train_end)
  test <- input_tsibble %>% filter_index(test_start ~ test_end)
  fcst_train <- input_tsibble %>% filter_index(train_start ~ test_end)
  return(list(train = train, test = test, fcst_train = fcst_train))
}

4) Hyperparameter Tuning with Grid Search

The grid_search function performs a grid search for hyperparameter tuning on various forecasting models, including mean, ETS, and ARIMA models. It identifies the optimal parameters that minimize error rates. This code will use a pre-defined set of models and hyperparameters to check for the lowest error. Each model and error will be put into a single dataframe for evaluation post processing. Then I have a function to grab the best configuration.

grid_search <- function(forecast_label, train, test) {
  model_tuning <- function(train, test, models, model_params) {
    results <- data.frame()
    for (i in seq(length(names(models)))) {
      model_name <- names(models[i])
      model_formula <- models[[model_name]]
      params <- model_params[[model_name]]
      if (length(params) > 0) {
        if (model_name == "mean") {
          for (param_value in params[["window_size"]]) {
            model_formula_updated <- MEAN(value ~ window(size = param_value))
            model <- train %>% model(!!model_name := model_formula_updated)
            fcsts <- model %>% forecast(h = glue::glue("{train_interval_months} months"))
            accuracy_result <- accuracy(fcsts, test) %>% select(.model, MAE, MAPE)
            accuracy_result[["window_size"]] <- param_value
            results <- bind_rows(results, accuracy_result)
          }
        }
      }
    }
    return(results)
  }
  models <- list(mean = MEAN(value ~ window(size = 6)), ets = ETS(value ~ error("A") + trend("A") + season("A")), arima = ARIMA(log(value) ~ pdq(0, 0, 2)))
  model_params <- list(mean = list(window_size = 1:12), ets = list(error = c("A", "M"), trend = c("A", "M", "Z"), seasonal = c("A", "M", "Z")))
  grid_search_results <- model_tuning(train, test, models, model_params) %>% mutate(run_date = format(Sys.time(), "%Y-%m-%d %H:%M:%S"), label = forecast_label)
  return(grid_search_results)
}

5) Getting Results and Generating Forecasts

Once the grid search runs, we need to pull back the best model. It’s always good to check different model outputs and see how the best of each model type works. This code will either return the very lowest model, or a specific one choosen by the user.

It begins by checking if the results dataframe is empty, issuing a warning if no results are found. The function can handle specific labels to extract the nth best model, validating the nth value to ensure it is within a valid range. It addresses infinite MAPE values by replacing them with Mean Absolute Error (MAE) values for comparison, ensuring models with infinite MAPE are still considered based on their MAE. The function filters results to retain only those with the lowest error metrics for each model, ensuring that the best model is selected based on the latest run date. It further refines the selection by checking if the top model for each label is a ‘mean’ model and if better alternatives are within 2% of its error rate, avoiding over-reliance on simpler models. Finally, it combines specific label results with general results for a comprehensive and accurate set of best models.

get_best_tuning_results <- function(gr_results, specific_label_nths = NULL) {
  if (nrow(gr_results) == 0) {
    warning("No results found for the given search term.")
    return(NULL)
  }

  if (!is.null(specific_label_nths)) {
    gr_results_specific <- list()
    for(label in names(specific_label_nths)){
      nth_best = specific_label_nths[[label]]
      temp_data <- gr_results %>% 
        dplyr::filter(label == !!label) %>%
        dplyr::group_by(label) %>% 
        dplyr::filter(run_date == max(run_date)) %>%
        dplyr::arrange(MAPE)
      if(nth_best < 1 | nth_best > nrow(temp_data)){
        warning(paste("Invalid nth_best value for label", label))
      } else {
        gr_results_specific[[label]] <- temp_data %>% slice(nth_best)
      }
    }
    gr_results_specific <- bind_rows(gr_results_specific)
  }
  
  if (any(is.infinite(gr_results$MAPE))) {
    gr_results <- gr_results %>%
      dplyr::mutate(error_metric = ifelse(is.infinite(MAPE), MAE, MAPE))
  } else {
    gr_results <- gr_results %>%
      dplyr::mutate(error_metric = MAPE)
  }

  gr_results <- gr_results %>%
    dplyr::filter(!is.na(error_metric)) %>%
    dplyr::group_by(label) %>% 
    dplyr::filter(run_date == max(run_date)) %>%
    dplyr::ungroup() %>% 
    dplyr::group_by(label, .model) %>%
    dplyr::filter(error_metric == min(error_metric)) %>%
    dplyr::arrange(error_metric) %>%
    dplyr::distinct(label, .keep_all = TRUE)

  gr_results <- gr_results %>%
    dplyr::group_by(label) %>%
    dplyr::arrange(error_metric) %>%
    dplyr::mutate(replace_mean = (.model == "mean") & (lead(error_metric) <= 1.02 * error_metric)) %>%
    dplyr::filter(!replace_mean | is.na(replace_mean)) %>%
    dplyr::distinct(label, .keep_all = TRUE)

  if (!is.null(specific_label_nths)) {
    gr_results <- gr_results %>% filter(!label %in% names(specific_label_nths))
    gr_results <- rbind(gr_results, gr_results_specific)
  }

  return(gr_results)
}

The create_forecasts function generates forecasts based on the best-tuned model identified earlier. It starts by constructing the appropriate model formula, selecting between MEAN, ETS, or ARIMA models based on the best model’s specifications. For MEAN models, it uses a window size parameter; for ETS models, it incorporates error, trend, and seasonal components; for ARIMA models, it includes parameters for differencing, autoregression, and moving average, with an option for a Fourier series term. The function then trains the model using the specified training data and generates forecasts for the defined forecast interval. These forecasts are transformed into a tidy dataframe, which includes the date and predicted value, ensuring the output is ready for analysis or visualization.

create_forecasts <- function(best_model, fcst_training_data) {
  model_formula <- if (best_model$.model == "mean") {
    MEAN(value ~ window(size = best_model$window_size))
  } else if (best_model$.model == "ets") {
    ETS(value ~ error(as.character(best_model$error)) + trend(as.character(best_model$trend)) + season(as.character(best_model$seasonal)))
  } else {
    ARIMA(log(value) ~ pdq(best_model$p, best_model$d, best_model$q) + fourier(K = best_model$K))
  }
  trained_model <- fcst_training_data %>% model(!!best_model$.model := model_formula)
  fcsts <- trained_model %>% forecast(h = glue::glue("{forecast_interval_months} months"))
  fcsts_out <- fcsts %>% as_tibble() %>% transmute(date = as.Date(year_month), value = .mean)
  return(list(fcst = fcsts_out, model = trained_model))
}