Time Series Analysis of Average Fuel Price (2000-2021)

The Data

The data used for this time series analysis is the U.S city average gasoline price data set sourced from the Federal Reserve Economic Data (FRED) ,St. Louis, and the United States Bureau of Labor Statistics. It contains monthly average regular unleaded gasoline prices from January 1976 to December 2021.The two fields in the data are date of observation and the average price.

When I first came across this data set on the FRED website, I was impressed with the amount of historical data recorded. I recently moved here from India, where I’ve driven for nearly 10 years, and I plan to start driving in the US as soon as possible, so analyzing the trend of gasoline prices seemed relevant and interesting to pursue.

Fuel price fluctuations are common across the world.The concern is when the rates increase considerably and consistently over a short time period. The most common reasons listed for price fluctuations are market competitors, tax and transportation fees, inflation, supply and demand, etc.

The goal of this analysis, as of now, is to understand how gasoline prices have fluctuated over the years and to determine whether: (i) prices follow patterns based on time of the year, (ii) the fluctuations can be captured through the data generation process, and (iii) a robust time series model can be fit to the data to forecast prices in a set time horizon.

The data set is continuous and does not have any gaps in monthly data points. It may not be easy to forecast the prices for very far into the future as there could be potential unforeseen changes in economic or social factors that affect the prices. However, predicting fluctuations for the near future should yield price values closer to actual prices.

For the purpose of this study, we will be considering data from January 2000 to December 2021.

The data set is loaded and filtered for the period of January 2000 till December 2021. The columns have been renamed in the final table for readability.

Exploratory Data Analysis

Structure of the Dataset

str(fueldata)

## tibble [264 × 3] (S3: tbl_df/tbl/data.frame)
##  $ obs_date  : Date[1:264], format: "2000-01-01" "2000-02-01" ...
##  $ avg_price : num [1:264] 1.3 1.37 1.54 1.51 1.5 ...
##  $ month_year: chr [1:264] "2000-01" "2000-02" "2000-03" "2000-04" ...

The data set considered for the study from the period of January, 2000 to December, 2021, contains 264 records. The average price is numeric and the original observed date field is of date type. The derived column month_year is a string column.

Summary of the Dataset

summary <- sumtable(fueldata,add.median = TRUE,out = 'kable')
summary %>% kable_styling()

Summary Statistics

Variable	N	Mean	Std. Dev.	Min	Pctl. 25	Pctl. 50	Pctl. 75	Max
avg_price	264	2.525	0.73	1.13	1.965	2.541	3.023	4.09

The mean average price is 2.525 which is almost comparable with the median of 2.541. This indicates that there is a possibility for the price to follow a normal distribution. This needs to be further examined through boxplots and histograms. No significant outliers are evident at this stage.

Trend of Price over Time

# Trend chart
monthly_price_plot = ggplot(fueldata)+
  geom_line(aes(obs_date,avg_price))+
  theme_bw()+
  xlab("Date")+
  ylab("Average Price")+
  labs(
    title = 'Average Price of Unleaded Regular Gasoline, FRED',
    subtitle = 'January 2000 - December 2021'
  )

#monthly_price_plot

monthly_price_plot + 
  geom_smooth(aes(obs_date,avg_price),method='lm',color='orange')

The trend chart of average price over time reveals a gradual positive change in fuel price over the years.The peak in average fuel price of $4.09, during the period of study, is observed in 2008 which maybe attributed to the impact of the global recession. There was a sharp fall in prices immediately after, in 2009, and then again in 2015 and 2016. The next notable decrease in prices was in early 2020 followed again by a steep rise.

#Average Price by Month of Year
fueldata %>% mutate(month = lubridate::month(obs_date,label = TRUE)) %>% group_by(month) %>% 
  summarize(avg_monthly_price = mean(avg_price)) %>% ungroup -> monthdata

avg_month_price_plot = ggplot(monthdata)+
  geom_col(aes(month,avg_monthly_price,fill=month))+
  theme_bw()+
  xlab("Month")+
  ylab("Average Price")+
  labs(
    title = 'Average Price of Unleaded Regular Gasoline by Month of the Year',
    subtitle = 'January 2000 - December 2021'
  )

avg_month_price_plot

The average fuel prices seem to be highest in the second and third quarters. The dip in prices during the winter months may be due to demand and supply fluctuations associated with extreme weather.

Histogram of Average Price

#Histogram of Average Fuel Price

p <- fueldata  %>%
  ggplot( aes(x=avg_price)) +
    geom_histogram(binwidth = 0.5, fill="#69b3a2", color="#e9ecef", alpha=0.9) +
    labs(
    title = 'Histogram: Average Price of Unleaded Regular Gasoline, FRED',
    subtitle = 'January 2000 - December 2021'
  )+
  xlab("Average Price")+
  ylab ("Frequency")
p

The distribution does not seem to follow a proper normal distribution. as there is a slight dip in frequency when the average price is between 1.75 and 2.25 before it peaks at around 2.5. Next step is to attempt a logarithmic transformation on the fuel price to see if that reveals a clearer distribution.

Histogram of Transformed Average Price

#Logarithmic Conversion and histogram plot 
p2 <- fueldata  %>%
  ggplot( aes(x=log(avg_price))) +
    geom_histogram(binwidth = 0.5, fill="#69b3a2", color="#e9ecef", alpha=0.9) +
    labs(
    title = 'Histogram: Logarithmic Average Price of Unleaded Regular Gasoline, FRED',
    subtitle = 'January 2000 - December 2021'
  )+
  xlab("Log(Average Price)")+
  ylab ("Frequency")
    
p2

The histogram for the logarithmic avergae fuel price reveals a slightly left skewed distribution.

Boxplot of Average Price

ggplot(fueldata, aes(y= avg_price)) +
  geom_boxplot(fill='#e9ecef', color="#69b3a2")+
  ylab("Average Price")+
  labs(
    title = 'Boxplot: Average Price of Unleaded Regular Gasoline, FRED',
    subtitle = 'January 2000 - December 2021'
  )+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())

The boxplot further corroborates the earlier observation that there are no significant outliers in the dataset.

Time Series Modelling - ARIMA

Behaviour of the Series

# Trend chart
monthly_price_plot = ggplot(fueldata)+
  geom_line(aes(obs_date,avg_price))+
  theme_bw()+
  xlab("Date")+
  ylab("Average Price")+
  labs(
    title = 'Average Price of Unleaded Regular Gasoline, FRED',
    subtitle = 'January 2000 - December 2021'
  )

#monthly_price_plot

monthly_price_plot + 
  geom_smooth(aes(obs_date,avg_price),method='lm',color='orange')

The trend chart of average price over time reveals a gradual positive change in fuel price over the years.The peak in average fuel price of $4.09, during the period of study, is observed in 2008 which maybe attributed to the impact of the global recession. There was a sharp fall in prices immediately after, in 2009, and then again in 2015 and 2016. The next notable decrease in prices was in early 2020 followed again by a steep rise.

The time series of the average gasoline price seems to be mean non-stationary for the period between 2000 and 2020 as there is an evident upward movement in the price over the years. On the other hand, the variance looks to be slightly fluctuating unevenly at some points in the period, which makes it difficult to ascertain stationarity. If the plot is variance non-stationary, it can be fixed using the box cox transformation.

Transformation for Variance Stationarity

fueldata_boxcox = fueldata %>%
  mutate(avg_price_boxcox = forecast::BoxCox(avg_price,lambda='auto'))
fueldata_boxcox %>%
  ggplot()+
      geom_line(aes(obs_date,avg_price_boxcox))+
      theme_bw()+
      ggtitle("Average Fuel Price over Time (Box-Cox Transformation)")+
      ylab("Average Price (Box-Cox Tranformation)")+
      xlab("Date")

After the box cox transformation to resolve any variance non-stationary behavior, the graph reveals no change visually. Hence, it can be concluded that the time series was originally variance stationary and does not require any transformation for variance stabilization.

Transformation for Mean Stationarity

We take the first lag difference of average price to attempt to stabilize the mean of the series.

fueldata_diff = fueldata%>%
  mutate(avg_price_boxcox = forecast::BoxCox(avg_price,lambda='auto')) %>%
  mutate(price_diff_box = avg_price_boxcox - lag(avg_price_boxcox)) %>%
  mutate(price_diff = avg_price - lag(avg_price))
fueldata_diff %>%
  ggplot()+
      geom_line(aes(obs_date,price_diff))+
      theme_bw()+
      ggtitle("Fuel Price (First Difference)")+
      ylab("Avg Price (Difference))")+
       xlab("Date")

The transformation of average price using the first difference seems to have brought the mean of the time series to 0. The next step is to test original as well as transformed variable to ascertain stationarity/non-stationarity. For this we will be using the Augmented Dicky-Fuller test and the KPSS test.

Tests for Non-Stationarity

ADF and KPSS Test

fueldata_diff = fueldata_diff %>%
  drop_na()

adf.test(fueldata_diff$avg_price) 

## 
##  Augmented Dickey-Fuller Test
## 
## data:  fueldata_diff$avg_price
## Dickey-Fuller = -2.1983, Lag order = 6, p-value = 0.4924
## alternative hypothesis: stationary

kpss.test(fueldata_diff$avg_price) 

## 
##  KPSS Test for Level Stationarity
## 
## data:  fueldata_diff$avg_price
## KPSS Level = 1.592, Truncation lag parameter = 5, p-value = 0.01

For the original price variable, the ADF test has showed a p value greater than 0.05 and the KPSS test value has a p value less than 0.05. This indicates that the average price is not mean stationary.

adf.test(fueldata_diff$price_diff) 

## 
##  Augmented Dickey-Fuller Test
## 
## data:  fueldata_diff$price_diff
## Dickey-Fuller = -7.8846, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationary

kpss.test(fueldata_diff$price_diff)

## 
##  KPSS Test for Level Stationarity
## 
## data:  fueldata_diff$price_diff
## KPSS Level = 0.04479, Truncation lag parameter = 5, p-value = 0.1

On testing for the differenced average price, the ADF test showed a p value less than 0.05 and KPSS test showed a p value above 0.05, indicating that the null hypothesis can be rejected and that the differenced variable is mean stationary.

adf.test(fueldata_diff$price_diff_box)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  fueldata_diff$price_diff_box
## Dickey-Fuller = -7.7629, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationary

kpss.test(fueldata_diff$price_diff_box)

## 
##  KPSS Test for Level Stationarity
## 
## data:  fueldata_diff$price_diff_box
## KPSS Level = 0.051774, Truncation lag parameter = 5, p-value = 0.1

On testing for the difference of box cox transformed variable, the same results were observed as for the raw differenced average price.
Since the original variable with difference is stationary as per the two tests, it can be concluded that, average price series is variance stationary and the first difference of average price is both mean and variance stationary.

Trend and Seasonality

The average price plotted over time earlier revealed an increasing trend over time. We take a closer look at the differenced average price from 2000 - 2010 and 2010-2020 to determine existence of seasonality.

fueldata_subset = fueldata_diff %>%
  filter(obs_date < ymd("2010=01-01") & obs_date >= ymd("2000=01-01")) 

fueldata_subset_2 = fueldata_diff %>%
  filter(obs_date < ymd("2021-01-01") & obs_date >= ymd("2010-01-01")) 

# Period : 2000-2010
fueldata_subset %>%
  ggplot()+
      geom_line(aes(obs_date,avg_price))+
      theme_bw()+
      ggtitle("Fuel Price (First Difference), 2000-2010")+
      ylab("Avg Price (Difference))") +
  xlab("Date")+
  scale_x_date(date_breaks = "3 months", date_labels =  "%b %Y") +theme(axis.text.x=element_text(angle=60, hjust=1))

# Period 2010-2020
fueldata_subset_2 %>%
  ggplot()+
      geom_line(aes(obs_date,avg_price))+
      theme_bw()+
      ggtitle("Fuel Price (First Difference), 2010-2020")+
      ylab("Avg Price (Difference))") +
   xlab("Date")+
  scale_x_date(date_breaks = "3 months", date_labels =  "%b %Y") + theme(axis.text.x=element_text(angle=60, hjust=1))

Yearly seasonality patterns in some form can be observed in the average fuel price. The price tends to be higher in the second quarter of the year and decreases considerably in the fourth quarter. Historically gasoline prices do tend to be higher during the spring-summer periods when demand is highest. In winter months, the demand for gasoline goes down due to bad weather and poor driving conditions. This could be the driving for the seasonal patterns observed.

However, there are also variations observed in the pattern which could be attributed to other economic factors such as recession, global events and other macro-economic factors. Examples of this is the drastic drops in fuel prices around the second quarter of 2020, when the first Covid-19 wave was in its climb and in 2008 during the global recession. Another example is the dip in the second quarter of 2016 which is attributed to the global drop in crude oil prices by the end of 2015.

ACF and PACF Plots of the Transformed Series

par(mfrow = c(1,2))
acf(fueldata_diff$price_diff_box,lag.max=30)
pacf(fueldata_diff$price_diff_box,lag.max=20)

The ACF plot of the first difference of average price shows resemblance of dampening over time in the first few periods, but it does not seem as evident as a pure AR model. It is possible that it is a mixed AR/MA model. If considered as an auto-regressive model, the PACF plot suggests two significant lags. Hence, the series can be considered to follow a second order auto-regressive model.

The analysis so far points towards the series following a ARIMA model of order (2,1,0) -a second order auto-regression element,with first order integration and no discernible behavior of a moving average model.

ARIMA

To determine whether it is a mixed model and the order of the auto-regression and moving average elements that best fits the time series, we try a few different variants as below. Since we were able to observe stationarity with first order differencing, we put the integration parameter as 1, and try different orders of AR and MA models.

arima(fueldata_diff$avg_price,order=c(0,1,2))

## 
## Call:
## arima(x = fueldata_diff$avg_price, order = c(0, 1, 2))
## 
## Coefficients:
##          ma1     ma2
##       0.5843  0.1532
## s.e.  0.0627  0.0651
## 
## sigma^2 estimated as 0.02013:  log likelihood = 139.69,  aic = -273.38

BIC(arima(fueldata_diff$avg_price,order=c(0,1,2)))

## [1] -262.679

arima(fueldata_diff$avg_price,order=c(1,1,2))

## 
## Call:
## arima(x = fueldata_diff$avg_price, order = c(1, 1, 2))
## 
## Coefficients:
##           ar1     ma1     ma2
##       -0.1839  0.7618  0.2431
## s.e.   0.2950  0.2841  0.1433
## 
## sigma^2 estimated as 0.0201:  log likelihood = 139.89,  aic = -271.78

BIC(arima(fueldata_diff$avg_price,order=c(1,1,2)))

## [1] -257.506

arima(fueldata_diff$avg_price,order=c(1,1,1))

## 
## Call:
## arima(x = fueldata_diff$avg_price, order = c(1, 1, 1))
## 
## Coefficients:
##          ar1     ma1
##       0.2147  0.3516
## s.e.  0.1026  0.0937
## 
## sigma^2 estimated as 0.02024:  log likelihood = 138.97,  aic = -271.93

BIC(arima(fueldata_diff$avg_price,order=c(1,1,1)))

## [1] -261.2297

arima(fueldata_diff$avg_price,order=c(2,1,2))

## 
## Call:
## arima(x = fueldata_diff$avg_price, order = c(2, 1, 2))
## 
## Coefficients:
##          ar1      ar2      ma1      ma2
##       1.1727  -0.3775  -0.6306  -0.1894
## s.e.  0.1123   0.1038   0.1186   0.1127
## 
## sigma^2 estimated as 0.01947:  log likelihood = 143.97,  aic = -277.94

BIC(arima(fueldata_diff$avg_price,order=c(2,1,2)))

## [1] -260.1002

arima(fueldata_diff$avg_price,order=c(2,1,1))

## 
## Call:
## arima(x = fueldata_diff$avg_price, order = c(2, 1, 1))
## 
## Coefficients:
##          ar1      ar2      ma1
##       1.2662  -0.5103  -0.7691
## s.e.  0.1010   0.0559   0.1071
## 
## sigma^2 estimated as 0.01967:  log likelihood = 142.64,  aic = -277.27

BIC(arima(fueldata_diff$avg_price,order=c(2,1,1)))

## [1] -262.9995

arima(fueldata_diff$avg_price,order=c(2,1,0))

## 
## Call:
## arima(x = fueldata_diff$avg_price, order = c(2, 1, 0))
## 
## Coefficients:
##          ar1      ar2
##       0.5766  -0.2313
## s.e.  0.0602   0.0602
## 
## sigma^2 estimated as 0.01999:  log likelihood = 140.57,  aic = -275.14

BIC(arima(fueldata_diff$avg_price,order=c(2,1,0)))

## [1] -264.4327

After trying different variations of AR and MA orders, we see that the AIC is lowest for the second order AR, first order integration, second order moving average model. However on calculating the BIC, the lowest value is observed for an ARIMA model of order (2,1,0). We will run an auto.arima to help determine best model.

auto.arima(fueldata_diff$avg_price,stationary=FALSE,allowdrift=FALSE,
           seasonal=FALSE,stepwise=FALSE,approximation=FALSE,max.order = 6,max.p = 7,max.q = 5)

## Series: fueldata_diff$avg_price 
## ARIMA(2,1,2) 
## 
## Coefficients:
##          ar1      ar2      ma1      ma2
##       1.1727  -0.3775  -0.6306  -0.1894
## s.e.  0.1123   0.1038   0.1186   0.1127
## 
## sigma^2 = 0.01977:  log likelihood = 143.97
## AIC=-277.94   AICc=-277.71   BIC=-260.1

Auto.arima above also reveals that the ARIMA(2,1,2) model is the best fit based on lowest AIC and BIC value. Therefore, we proceed with the model ARIMA(2,1,2) for residual and in-sample performance testing.

Residuals and In-Sample Performance Test

best_mod = arima(fueldata_diff$avg_price,order=c(2,1,2))
resid = best_mod$residuals
pred = resid+fueldata_diff$avg_price
ggplot()+
  geom_line(aes(fueldata_diff$obs_date,fueldata_diff$avg_price))+
  geom_line(aes(fueldata_diff$obs_date,pred),color='blue',alpha=0.4)+
  theme_bw()+
  xlab("Date")+
  ylab("Avg Price")

The in sample performance of the model is visualized above. The predicted values of the model for the same time period follow the same pattern of the observed values. Almost all fluctuations in observed values of the series are captured by the model predictions.

RMSE_insample = sqrt(mean((expm1(pred) - expm1(fueldata_diff$avg_price))^2,na.rm=T))
RMSE_insample

## [1] 3.098569

The RMSE for the in sample predictions is also favorable to the model selected. The RMSE reveals that the model predicts fuel prices within approximately 3 points of actual price on average in-sample.

forecast::checkresiduals(best_mod)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,1,2)
## Q* = 6.6341, df = 6, p-value = 0.356
## 
## Model df: 4.   Total lags used: 10

par(mfrow=c(1,2))
resid = best_mod$resid
acf(resid,lag.max=20)
pacf(resid,lag.max=20)

The residuals of the selected ARIMA model visually resemble a white noise process. However, from the histogram, the residuals do not exactly follow a normal distribution on account of a slight left skew caused by the large dip in the average price in 2008.

The ACF and PACF plots show that the autocorrelation function for the series is generally very close to 0.

Ljung - Box Test for Residual Independence

We further test for auto-correlation using the Ljung- Box test, to determine if the residuals are independent as expected.

Box.test(resid,type='Ljung-Box',lag=1)

## 
##  Box-Ljung test
## 
## data:  resid
## X-squared = 0.0050491, df = 1, p-value = 0.9434

Box.test(resid,type='Ljung-Box',lag=5)

## 
##  Box-Ljung test
## 
## data:  resid
## X-squared = 0.80015, df = 5, p-value = 0.977

Box.test(resid,type='Ljung-Box',lag=10)

## 
##  Box-Ljung test
## 
## data:  resid
## X-squared = 6.6341, df = 10, p-value = 0.7595

Box.test(resid,type='Ljung-Box',lag=20)

## 
##  Box-Ljung test
## 
## data:  resid
## X-squared = 17.039, df = 20, p-value = 0.6505

The test is performed for different lags. The p value for lags 1,5,10 and 20 are all greater than 0.05. This indicates that we can reject the null hypothesis that the residuals are not independent. We conclude that the residuals for the best model selected are independent and resemble a white noise process.

Forecast of Average Price

prediction = predict(best_mod,n.ahead=6)

pred_data = data.frame(
  pred = prediction,
  obs_date = c(ymd(max(fueldata_diff$obs_date)),ymd(max(fueldata_diff$obs_date)+ months(1)),ymd(max(fueldata_diff$obs_date)+ months(2)),ymd(max(fueldata_diff$obs_date)+ months(3)),ymd(max(fueldata_diff$obs_date)+ months(4)),ymd(max(fueldata_diff$obs_date)+ months(5)))
)


pred_merge = fueldata_diff %>% 
  full_join(pred_data) %>%
  mutate(
    pred_high = pred.pred+2*pred.se,
    pred_1se_high = pred.pred+1*pred.se,
    pred_low = pred.pred - 2*pred.se,
    pred_1se_low = pred.pred-1*pred.se,
    pred.pred = pred.pred,
    error = pred.pred - avg_price
    )

pred_merge %>% 
  filter(obs_date >= ymd("2000-06-01")) %>%
  ggplot()+ 
  geom_line(aes(obs_date,avg_price))+
  geom_line(aes(obs_date,pred.pred),color='blue')+
  geom_ribbon(aes(x=obs_date,ymin=pred_low,ymax = pred_high),fill='blue',alpha=0.2)+
  geom_ribbon(aes(x=obs_date,ymin=pred_1se_low,ymax = pred_1se_high),fill='blue',alpha=0.2)+
  theme_bw()+
  xlab("Date")+
  ylab("Avg Price")

The predictions from the selected best model for five months into the future starting from the end of 2021 are as per the graph above. The predicted average price for the end of 2021 seem to follow the decreasing trend observed during the last quarter of 2021 and the values predicted are not extremely different from the previous period values. Thus the conclusion is that the projected values look reasonable in the period which we have forecast for.

Time Series Modelling - Prophet Model

We will now attempt to fit a prophet model to the time series, explore seasonality and trend, and assess the model fit and performance through error metrics.

We split the time series data into train data - from January 2000 to December 2014 - and test data - from January 2015 to January 2021. We then build the prophet model on the train data and forecast. Since this series is captured at month level, the future forecast period units have been modified accordingly.

prophet_fueldata = fueldata %>% 
    rename(ds = obs_date,
    y = avg_price)  

train = prophet_fueldata %>% 
  filter(ds<ymd("2015-01-01"))

test = prophet_fueldata %>%
  filter(ds>=ymd("2015-01-01"))

model = prophet(train)
future = make_future_dataframe(model,periods = 85, freq = 'month')
forecast = predict(model,future)

plot(model,forecast)+
ylab("Average Fuel Price")+xlab("Date")+labs(title = 'Forecast: Average Price of Unleaded Regular Gasoline, FRED',subtitle = 'January 2000 - December 2021')+theme_bw()

The model has a very steep upward trend, following the trend of the original time series. Since the model follows the default additive seasonality pattern, it fails to capture some of the the higher fluctuations from 2005 onward. Judging from the light blue overlay which shows the variability in the predicted values, the maximum value of average price at the end of the forecast peiod is estimated to be above $5, which is not as observed in the actual data. We may need to tune the model in terms of seasonality to get forecast behaviour closer to reality. The actual and predicted values can be examned more closely with the ineractive plot below.

Time Series Components

We plot the components of the series separately. Since the data is at a monthly level, no weekly seasonality information is available. Additionally, we will not be able to incorporate holiday information either as it does not apply to this level of data.

prophet_plot_components(model,forecast)

The trend component shows an general upward trend in the average price. The seasonality at yearly level reveals that the average price of gasoline tends to increase around the March and October in anticipation of the increased demand in Summer and after the onset of Fall. The dip in average price is the most in the months of January to February and August to September.

Determination of Saturation Points

We look at summary statistics of the train and test data to determine whether there is a reasonable limit within which the fuel prices tend to stay. We have 20 years of data in total which should be a reliable estimate of the range of the fuel price. We will compare the forecast fuel price and the trend of the forecast to determine whether we need to apply a minimum and maximum saturation point for more realistic estimates.

# Using summary statistics to determine the range of values for train and test data
summary(train)

##        ds                   y          month_year       
##  Min.   :2000-01-01   Min.   :1.130   Length:180        
##  1st Qu.:2003-09-23   1st Qu.:1.655   Class :character  
##  Median :2007-06-16   Median :2.591   Mode  :character  
##  Mean   :2007-06-16   Mean   :2.530                     
##  3rd Qu.:2011-03-08   3rd Qu.:3.323                     
##  Max.   :2014-12-01   Max.   :4.090

summary(test)

##        ds                   y          month_year       
##  Min.   :2015-01-01   Min.   :1.767   Length:84         
##  1st Qu.:2016-09-23   1st Qu.:2.240   Class :character  
##  Median :2018-06-16   Median :2.474   Mode  :character  
##  Mean   :2018-06-16   Mean   :2.513                     
##  3rd Qu.:2020-03-08   3rd Qu.:2.775                     
##  Max.   :2021-12-01   Max.   :3.482

# Set "floor" in training data
train$floor = 1.10
train$cap = 4.10

# Set floor in forecast data
future$floor = 1.70
future$cap = 3.5
sat_model = prophet(train,growth='logistic')
sat_forecast = predict(sat_model,future)

plot(sat_model,sat_forecast)+ylim(0,4.5)+xlab("Date")+ylab("Avg Price")+labs(title = 'Forecast: Average Price of Unleaded Regular Gasoline, FRED ',subtitle = 'Saturation Points')+theme_bw()

The summary statistics reveal that the values tend to stay between $1.1 and $4.1 for the training period while the interval is tighter - between $1.7 and $3.5 - in the test data for the forecast period. Therefore, we set these saturation thresholds for the prophet model and the resultant plot reveals that the forecast values are very much different from the earlier forecast and does not show as steep an upward trend as prior to setting the limits.

Examination for Changepoints

We will examine the series for significant change points in trend of the series using the training data.

model = prophet(train, n.changepoints = 10)
future = make_future_dataframe(model,periods = 85, freq = 'month')
plot(model,forecast)+add_changepoints_to_plot(model)+xlab("Date")+ylab("Avg Price")+labs(title = 'Forecast: Average Price of Unleaded Regular Gasoline, FRED',subtitle = 'January 2000 - December 2021')+theme_bw()

We observe that only two significant change-points have been identified for this series indicating where the trend changes drastically.

Type of Seasonality - Additive or Multiplicative?

additive = prophet(train)
add_fcst = predict(additive,future)

plot(additive,add_fcst)+
ylim(0,5)+xlab("Date")+ylab("Avg Price")+labs(title = 'Forecast: Average Price of Unleaded Regular Gasoline, FRED',subtitle = 'Additive Seasonality Model')+theme_bw()

prophet_plot_components(additive,add_fcst)

The additive model does not seem to closely fit the higher seasonal fluctuations in average price values towards the later part of the training period but behaves well for the earlier time frame.

multi = prophet(train,seasonality.mode = 'multiplicative')
multi_fcst = predict(multi,future)

plot(multi,multi_fcst)+ylim(0,5) +theme_bw()+xlab("Date")+ylab("Avg Price")+labs(title = 'Forecast: Average Price of Unleaded Regular Gasoline, FRED ',subtitle = 'Multiplicative Seasonality Model')

prophet_plot_components(multi,multi_fcst)

The multiplicative model visibly seems to capture the higher fluctuations towards the end of the period but does not fit the earlier period with lower fluctuations well.

From visual assessment, the additive model seems likely to perform better overall, especially since we are not able to conclude with certainty on the existence of multiplicative seasonality.

Comparison of Additive versus Multiplicative Fit for Seasonality

mod1 = prophet(train,seasonality.mode='additive')
forecast1 = predict(mod1)
df_cv1 <- cross_validation(mod1, horizon=365,period = 365,initial = 5*365,units='days')
metrics1 = performance_metrics(df_cv1) %>% 
  mutate(model = 'mod1')

mod2 = prophet(train,seasonality.mode='multiplicative')
forecast2 = predict(mod2)
df_cv2 <- cross_validation(mod2, horizon=365,period = 365,initial = 5*365,units='days')
metrics2 = performance_metrics(df_cv2) %>% 
  mutate(model = "mod2")
metrics1 %>% 
bind_rows(metrics2) %>% 
ggplot()+
geom_line(aes(horizon,rmse,color=model)) +theme_bw()+labs(title = 'RMSE Comparison for Additive and Multiplicative Model Fits  ',subtitle = 'Data: Average Price of Unleaded Regular Gasoline, FRED')

The forecast RMSE can be observed as generally much higher for the multiplicative fit model which further supports the earlier indication that the seasonality element in the series of average fuel price tends to follow an additive behavior. Therefore, we proceed with the additive Prophet model as the appropriate choice for modelling our time series.

Prophet Model Assessment

We calculate the error metrics to assess model performance.

forecast_metric_data = forecast %>% 
  as_tibble() %>% 
  filter(ds>=ymd("2015-01-01"))
RMSE = sqrt(mean((test$y - forecast_metric_data$yhat)^2))
MAE = mean(abs(test$y - forecast_metric_data$yhat))
MAPE = mean(abs((test$y - forecast_metric_data$yhat)/test$y))
print(paste("RMSE:",round(RMSE,2)))

## [1] "RMSE: 1.81"

print(paste("MAE:",round(MAE,2)))

## [1] "MAE: 1.77"

print(paste("MAPE:",round(MAPE,2)))

## [1] "MAPE: 0.73"

The Prophet model provides an acceptable fit for the time series data, even though the frequency is at month level. The model has acceptable ranges of mean absolute percent error and rmse and also allows for a qualitative and quantitative assessment of the seasonality effect, thus enabling more accurate forecasts.

Model Comparison : ARIMA vs Prophet

Cross Validation and Error Metrics - ARIMA

xlist <-c(1,13,25,37,49,61,73,85,97)
out_list = lapply(xlist,function(x){
  train = fueldata[x:(x+59),]
  print(train)
  test = fueldata[(x+60):(x+60+11),]
  print(test)
 
  best_mod_train = arima(train$avg_price,order=c(2,1,2))
  out_data = data.frame(
    train_thru = x,
    time_ahead = 1:12,
    time_point = (x+59)+ 1:12,
    actual = test$avg_price,
    naive_fcst = train$avg_price[x],
    
    test_pred= predict(best_mod_train,12)$pred,
   
    error = test$avg_price - predict(best_mod_train,12)$pred
  )
 })
out_data = rbindlist(out_list)

out_data%>% ggplot()+
geom_point(aes(time_point,actual)) +
  geom_point(aes(time_point,test_pred,colour = factor(train_thru)))+xlab("Time Ahead")+ylab("Avg Fuel Price")+labs( title = 'Forecast: Average Price of Unleaded Regular Gasoline, FRED ',
    subtitle = 'Cross Validation- ARIMA') +theme_bw()+theme(legend.position = "None")

The plot shows how the predictions for the prescribed horizon fares against the actuals at a month level.

error_metrics = out_data %>% 
  group_by(time_ahead) %>% 
  summarize(
    RMSE = sqrt(mean(error^2)),
    MAE = mean(abs(error)),
    MAPE = mean(abs(error / actual))
    ) %>%
    ungroup()
em_melt <- melt(error_metrics, id.vars='time_ahead',variable.name = "errortype",value.name = "error")
em_melt %>% ggplot()+
geom_line(aes(time_ahead,error,color=errortype)) +theme_bw() + scale_x_continuous(breaks = c(1:12),limits = c(0, 13))+ scale_y_continuous(limits = c(0,0.75)) +xlab("Time Ahead")+ylab("Error") +labs(title = 'Error Comparison for ARIMA CV Fit',subtitle = 'Data: Average Price of Unleaded Regular Gasoline, FRED')

The root mean squared error ranges from 0.19 and 0.68 for the cross validation results on ARIMA model with a horizon of one year. The mean absolute error is 0.15 and 0.49. The mean absolute percentage error ranges from 5% to 20%.

Cross Validation and Error Metrics - Prophet

We will perform model cross-validation by making the necessary changes to accommodate monthly data. After multiple combinations, setting the initial training period as 5 years and the rolling horizon as 1 year seems to have reasonably low error metrics.

df.cv <- cross_validation(model ,horizon=365,period = 365,initial = 5*365,units='days')
head(df.cv)

## # A tibble: 6 × 6
##       y ds                   yhat yhat_lower yhat_upper cutoff             
##   <dbl> <dttm>              <dbl>      <dbl>      <dbl> <dttm>             
## 1  2.32 2006-01-01 00:00:00  2.33       2.19       2.48 2005-12-03 00:00:00
## 2  2.31 2006-02-01 00:00:00  2.42       2.26       2.57 2005-12-03 00:00:00
## 3  2.40 2006-03-01 00:00:00  2.56       2.42       2.71 2005-12-03 00:00:00
## 4  2.76 2006-04-01 00:00:00  2.65       2.52       2.79 2005-12-03 00:00:00
## 5  2.95 2006-05-01 00:00:00  2.66       2.52       2.80 2005-12-03 00:00:00
## 6  2.92 2006-06-01 00:00:00  2.66       2.52       2.82 2005-12-03 00:00:00

df.cv %>% 
  ggplot()+
  geom_point(aes(ds,y)) +
  geom_point(aes(ds,yhat,color=factor(cutoff)))+
  theme_bw()+
  xlab("Date")+
  ylab("Avg Price")+
  scale_color_discrete(name = 'Cutoff')+
  theme(legend.position = "None")+
  labs(
    title = 'Forecast: Average Price of Unleaded Regular Gasoline, FRED ',
    subtitle = 'Cross Validation - Prophet'
  )

plot_cross_validation_metric(df.cv, metric = 'rmse')

plot_cross_validation_metric(df.cv, metric = 'mape')

plot_cross_validation_metric(df.cv, metric = 'mae')

The cross validation result plot reveals how the rolling window predictions compare against actuals.
The RMSE ranges between 0.5 and 0.6 for the forecast horizon of 1 year. The mean absolute error is between 0.3 and 0.5. The mean absolute percent error ranges between 12 and 22 percent.

Based on the results, the ARIMA model seems to have lower errors for predictions in the near future while higher error later in the periods further into the future are comparable with the same for Prophet. Prophet has higher errors for near future predictions but the variability in the error is less. Based on the RMSE, and the observations above, we see that the ARIMA performs better than Prophet for this time series.

Conclusion

The forecast from ARIMA looks reasonable as it doesn’t show any concerning changes in direction, peaks or dips that indicate the model is not estimating the average price values correctly. We can conclude that the model has captured the data generating process reasonably well and can explore methods to improve the model performance to an extend using more data which also includes external regression variables that could help in explaining the variability of gasoline price better. Thus we have addressed all areas discussed such as deciphering whether the average price follow patterns based on time of the year and whether the fluctuations can be captured through the data generation process. We also fit two distinct time series models to the data to forecast prices in a set time horizon and assessed performance based on relevant metrics to identify the better model. This concludes the study focusing on time series analysis on average fuel(gasoline) prices from 2000-2021.