Objective:

To forecast West Texas Intermediate (WTI) crude oil prices over a 12-month horizon (January to December 2023).

Given the significant influence of WTI prices on Philippine fuel and commodity costs, accurate forecasting is critical for:

Budgeting and procurement planning
Strategic decision-making for cost control
Risk mitigation amid volatile market conditions

Forecast Horizon: January 2023 - December 2023

Disclaimer: This forecast is based on statistical time series models trained on historical WTI data (1980–2022). These models are designed to capture:

Long-term trends
Seasonal patterns
Internal residual structures

However, they do not account for external shocks, such as:

Global conflicts or political instability
Unexpected supply chain disruptions
OPEC+ decisions or production shifts
Natural disasters or pandemics

As such, these forecasts should be interpreted as data-driven projections, not absolute predictions.

# Dataset
library(tidyverse)
library(fpp3)
library(readr)
library(fabletools)
library(feasts)
library(purrr)
library(tseries)
library(tinytex)
library(tsibbledata)
library(dplyr)
library(forecast)
library(purrr)
library(fable)

wti <- read_csv("C:/Users/Administrator/Documents/0. Data Science/Forecasting Analytics - Data Science Assessment/wti.csv")

#Preview the dataset
glimpse(wti)

## Rows: 516
## Columns: 2
## $ DATE    <date> 1980-01-01, 1980-02-01, 1980-03-01, 1980-04-01, 1980-05-01, 1…
## $ WTISPLC <dbl> 32.50, 37.00, 38.00, 39.50, 39.50, 39.50, 39.50, 38.00, 36.00,…

summary(wti)

##       DATE               WTISPLC      
##  Min.   :1980-01-01   Min.   : 11.28  
##  1st Qu.:1990-09-23   1st Qu.: 20.80  
##  Median :2001-06-16   Median : 33.77  
##  Mean   :2001-06-16   Mean   : 44.18  
##  3rd Qu.:2012-03-08   3rd Qu.: 61.80  
##  Max.   :2022-12-01   Max.   :133.93

# Time Series Preparation and Cleaning
wti_ts <- wti |> 
  mutate(DATE = yearmonth(DATE)) |> 
  as_tsibble(index = DATE)

wti_ts

## # A tsibble: 516 x 2 [1M]
##        DATE WTISPLC
##       <mth>   <dbl>
##  1 1980 Jan    32.5
##  2 1980 Feb    37  
##  3 1980 Mar    38  
##  4 1980 Apr    39.5
##  5 1980 May    39.5
##  6 1980 Jun    39.5
##  7 1980 Jul    39.5
##  8 1980 Aug    38  
##  9 1980 Sep    36  
## 10 1980 Oct    36  
## # ℹ 506 more rows

has_gaps(wti_ts)

## # A tibble: 1 × 1
##   .gaps
##   <lgl>
## 1 FALSE

count_gaps(wti_ts)

## # A tibble: 0 × 3
## # ℹ 3 variables: .from <mth>, .to <mth>, .n <int>

interval(wti_ts)

## <interval[1]>
## [1] 1M

Time series has monthly frequency with regular intervals
No missing periods (complete and clean)
Ready for time series plotting

Exploratory Data Analysis

autoplot(wti_ts, WTISPLC) +
  labs(title = "West Texas Intermediate SPLC", y = "WTISPLC", x = "Date")

Initial Observations

Strong trends (upward: 2000s & 2020; downward: 1980s)
Oil shocks (1990 Gulf War, 2008 crash, 2020 pandemic)

Seasonality Check

wti_ts |> 
  gg_season(WTISPLC) +
  labs(title = "Seasonality Check", y = "WTISPLC", x = "Date")

Weak visible seasonality, but moderate statistical presence (to be validated using STL & ACF)

wti_ts |>
  gg_subseries(WTISPLC) +
  labs(title = "Quarterly Patterns", y = "WTISPLC", x = "Date")

Q2–Q3 tends to have higher values in certain years

STL Decomposition

Decomposition breaks a time series into three components:

Trend – the long-term direction or structural movement. (e.g., rising oil prices)
Seasonality – Repeating calendar-based patterns (e.g., quarterly or monthly effects)
Remainder (Noise) – Irregular events not explained by trend or seasonality (e.g., political events)

This helps to understand the drivers of the data, making it easier to choose and explain forecasting models.

dcmp <- wti_ts |>
  model(
    STL(WTISPLC ~ season(window = "periodic"), robust = TRUE)
  )
components(dcmp) |>
  autoplot() +
  labs(title = "STL Decomposition of WTI Spot Prices")

The STL decomposition reveals that WTI prices are primarily trend-driven, with moderate seasonality and visible irregular components.

This highlights the need for models that capture strong trends, mild seasonal cycles, and allow for occasional volatility.

wti_ts |>
  autoplot(WTISPLC, color = "gray") +
  autolayer(components(dcmp), trend, color = "#D55E00") +
  labs(title = "WTI Prices with STL Trend", y = "WTISPLC", x = "Date")

The gray line shows the raw, monthly WTI spot prices
The orange/red line is the STL trend component, which extracts and smooths the long-term direction of price movement.
The trend reveals major cycles such as:

📉 2008 Global Financial Crisis – sharp decline in oil demand.

📉 2020 COVID-19 pandemic – steep crash due to travel and economic restrictions.

📈 Post-pandemic recovery (2021–2022) – rebound in global demand.

✅ Why it matters: Useful for strategic planning, e.g., investment timing, budget forecasting, or policy shifts.

wti_ts |>
  autoplot(WTISPLC, color = "gray") +
  autolayer(components(dcmp), season_adjust, color = "#0072B2") +
  labs(title = "Seasonally Adjusted WTI Prices", y = "WTISPLC", x = "Date")

The blue line shows WTI prices with seasonality removed, making external shocks easier to see.
Compared to the trend plot, the raw line is less visually dominant, focusing attention on short-term volatility.
This plot reveals non-seasonal shocks more clearly, such as:
- The 2020 oil price crash, where demand collapsed globally.
- Sudden spikes during geopolitical tensions or OPEC+ production cuts.
✅ Why it matters: Great for impact analysis when understanding non-seasonal drivers of price movement.

stl_feats <- wti_ts |>
  features(WTISPLC, feat_stl)

print(stl_feats, width = Inf)

## # A tibble: 1 × 9
##   trend_strength seasonal_strength_year seasonal_peak_year seasonal_trough_year
##            <dbl>                  <dbl>              <dbl>                <dbl>
## 1          0.965                  0.196                  7                    2
##   spikiness linearity curvature stl_e_acf1 stl_e_acf10
##       <dbl>     <dbl>     <dbl>      <dbl>       <dbl>
## 1    0.0248      422.      91.3      0.749        1.14

Trend strength = 0.965 → Strong
Seasonality strength = 0.196 → Weak to moderate
Peak = July, Trough = February
Spikiness = low → smoother series

This confirms that trend dominates the WTI series, while seasonality exists but is not overwhelming.

Stationarity and Autocorellation

Purpose of checking if the dataset is stationary,

A stationary time series has a constant mean and variance over time, and no predictable trends or seasonality.
If the series is non-stationary, it often exhibits trends, seasonality, or varying variance, which must be explicitly modeled or removed (e.g., through differencing) before forecasting.

gg_tsdisplay(wti_ts, WTISPLC, plot_type = "partial")

The time series plot shows a visible trend over time.
The ACF decays slowly and remains significant over several lags.
The PACF plot also shows persistent lags.
These patterns strongly suggest non-stationarity, likely due to trend and weak seasonality.

#ADF Test
adf.test(wti_ts$WTISPLC, alternative = "stationary")

## 
##  Augmented Dickey-Fuller Test
## 
## data:  wti_ts$WTISPLC
## Dickey-Fuller = -2.93, Lag order = 8, p-value = 0.1847
## alternative hypothesis: stationary

❌ p-value = 0.1847 > 0.05 → Fail to reject the null hypothesis of a unit root
This means the data is non-stationary.

We conducted the ADF test to statistically confirm stationarity. With a p-value greater than 0.05, we conclude that the WTI price series is non-stationary, indicating that it has underlying trends or patterns over time. This is crucial because non-stationary data can lead to unreliable forecasts unless properly transformed or modeled.
✅ Differencing is required

# First Differencing
diff1 <- wti_ts |>
  mutate(diff1 = difference(WTISPLC, 1))

gg_tsdisplay(diff1, diff1, plot_type = "partial")

The differenced time plot shows improved stationarity, especially before 2010. However, some shifts remain visible, suggesting possible structural breaks or external shocks.
The ACF drops sharply after lag 1 but some later lags still exceed the confidence bounds.
PACF same with ACF.The PACF pattern is similar.

These suggest that the series is closer to stationarity, but further confirmation is needed via a statistical test.

# First differencing (creates 1 NA at the start)
diff1_vec <- difference(wti_ts$WTISPLC)

# Remove the NA before ADF test
diff1_vec_clean <- na.omit(diff1_vec)

# Run ADF test on cleaned differenced data
adf.test(diff1_vec_clean, alternative = "stationary")

## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff1_vec_clean
## Dickey-Fuller = -8.4837, Lag order = 8, p-value = 0.01
## alternative hypothesis: stationary

✅ p-value < 0.05 → Reject the null hypothesis

📈 Interpretation: First differencing achieved stationarity

# Second difference
diff2 <- wti_ts |>
  mutate(diff2 = difference(WTISPLC, 2))

gg_tsdisplay(diff2, diff2, plot_type = "partial")  # Second difference

The time plot appears flattened and overcorrected, with little variability.
ACF and PACF plots show minimal spikes, suggesting a loss of structure (i.e., close to white noise).

This suggests that second differencing may be unnecessary and could result in over-differencing, removing important signal needed for forecasting.

acf_features <- wti_ts |> 
  features(WTISPLC, feat_acf)

print(acf_features)

## # A tibble: 1 × 7
##    acf1 acf10 diff1_acf1 diff1_acf10 diff2_acf1 diff2_acf10 season_acf1
##   <dbl> <dbl>      <dbl>       <dbl>      <dbl>       <dbl>       <dbl>
## 1 0.986  7.81      0.352       0.191     -0.314       0.143       0.767

The ACF shows a strong autocorrelation at lag 1 (0.986), which is typical for a trending series. While ACF/PACF are not formal stationarity tests, they help detect persistent structures that may require differencing or modeling.
First differencing (ACF = 0.352) autocorrelation significantly reduced, while second differencing (ACF = -0.314) Further differencing removes structure, likely over-differencing
The series exhibits trend and structure, not white noise — first differencing is appropriate.

Based on the exploratory data analysis, the WTI series displays a strong long-term trend with weak to moderate seasonality, confirmed through STL decomposition and seasonal diagnostics.

While the time series is non-stationary in its raw form (as shown by visual trends, persistent ACF, and ADF test results), applying first differencing yields a stationary series suitable for forecasting, as supported by the improved ACF/PACF behavior and the ADF test (p-value < 0.05).

Second differencing appears unnecessary and risks removing valuable structure. With these characteristics established — strong trend, moderate seasonality, and confirmed stationarity post-differencing — we are now in a statistically sound position to proceed with developing and evaluating forecasting models.

Forecasting

Why Forecast WTI Spot Prices?

Forecasting WTI (West Texas Intermediate) spot prices is critical due to their significant impact on economic policy, energy planning, and investment decisions. Accurate forecasts help stakeholders anticipate price movements, manage risk, and inform strategic actions in industries sensitive to oil price volatility.

Models

To better understand and predict the behavior of WTI spot prices, we explored a range of forecasting approaches. These models were selected to address different aspects of the time series, such as trend, seasonality, and randomness. For clarity, the methods are grouped into two broad categories:

Benchmark Models (Baselines for Comparison)

These simple models serve as reference points to evaluate the added value of more complex forecasting techniques.
- Naïve Model (NAIVE): Assumes the next value equals the most recent observation — suitable for random walk behavior.
- Drift Model (Random Walk with drift): Repeats values from the same season in the previous year — effective for highly seasonal data.
- Seasonal Naïve Model (SNAIVE): Assumes that each future value repeats the value from the same period in the previous year. Useful for data with strong seasonal patterns.
- Seasonal Naïve with Drift: This hybrid approach is not statistically coherent due to conflicting assumptions between drift and seasonality. It is not supported in the {fable} framework and was excluded from the analysis.
Advanced Models: These models attempt to capture deeper structure in the data — including autoregressive behavior, trends, and seasonal patterns — to improve both forecast accuracy and interpretability.
- Exponential Smoothing (ETS): Automatically selects the optimal combination of error, trend, and seasonal components. Flexible for both additive and multiplicative dynamics.
- ARIMA (AutoRegressive Integrated Moving Average): Combines autoregressive terms, differencing, and moving averages. Can be extended to SARIMA to account for seasonal effects.
- ARDL (AutoRegressive Distributed Lag) (For Future Work): A model that incorporates both past values of the target variable and lagged exogenous variables (e.g., inflation, interest rates). While not yet implemented, ARDL offers potential for richer forecasting once relevant macroeconomic variables are available.

Modeling Strategy

We adopted a stepwise approach, starting with benchmark models (Naïve, Drift, Seasonal Naïve) to set a performance baseline. These were followed by advanced models (ETS, ARIMA), designed to capture more complex dynamics like trend, seasonality, and autocorrelation. This progression allows a clear comparison of model improvements and ensures model complexity is justified by measurable gains in accuracy.

Model Evaluation Criteria

To assess model performance, we used a combination of forecast accuracy metrics and residual diagnostics:

Accuracy Metrics
- RMSE (Root Mean Squared Error): Penalizes large errors more heavily
- MAE (Mean Absolute Error): Measures average magnitude of errors
- MAPE (Mean Absolute Percentage Error): Useful for comparing relative accuracy
Residual Diagnostics
- Ljung-Box test: Checks for autocorrelation in residuals
- ACF plots: Visualize lag dependence
- Normality & Homoscedasticity: Ensures residuals behave like white noise

Benchmark Models

Naive

Assumptions:

Assumes the next value is equal to the most recent observation.
Suitable as a benchmark model for time series data with random walk behavior (no trend or seasonality).

model_naive <- wti_ts |>
  model(NAIVE(WTISPLC))

model_naive |> gg_tsresiduals()

Residuals: Show visible structure and autocorrelation.
ACF Plot: Several bars extend beyond the blue dashed lines (confidence bounds). This indicates autocorrelation in residuals, meaning there’s still predictable structure the model didn’t capture — violating the white noise assumption.
Histogram of Residuals: The histogram is skewed, bumpy, or multimodal. Residuals are not normally distributed, reinforcing that the model doesn’t account for all patterns — possibly due to missing trend/seasonality.

model_naive |>
  augment() |>  # adds `.resid`, `.innov`, etc.
  features(.innov, ljung_box, lag = 10)

## # A tibble: 1 × 3
##   .model         lb_stat lb_pvalue
##   <chr>            <dbl>     <dbl>
## 1 NAIVE(WTISPLC)    99.3  1.11e-16

Residual Diagnostics

Ljung-Box test p-value: 1.11e-16
✅ Result: ❌ Residuals are not white noise
Naive model fails the Ljung-Box test, confirming that there’s still structure in the series.

Insights

❌ Fails to capture trends or seasonal structures
✅ Useful as a baseline to compare performance of more complex models
⚠️ Significant autocorrelation in residuals → suggests unmodeled structure

Drift

Assumptions:

Forecasts the next value as the last observed value + a constant drift term
Captures linear trend over time
Does not account for seasonality

model_drift <- wti_ts |>
  model(RW(WTISPLC ~ drift()))

model_drift |> gg_tsresiduals()

Residuals: Autocorrelated and non-random.
ACF Plot: Several lags exceed bounds → autocorrelation is present.
Histogram of Residuals: Histogram is bumpy or skewed → residuals are not normally distributed.

model_drift |>
  augment() |>  # adds `.resid`, `.innov`, etc.
  features(.innov, ljung_box, lag = 10)

## # A tibble: 1 × 3
##   .model                lb_stat lb_pvalue
##   <chr>                   <dbl>     <dbl>
## 1 RW(WTISPLC ~ drift())    99.3  1.11e-16

Residual Diagnostics

Ljung-Box test p-value: 1.11e-1
✅ Result: ❌ Residuals are not white noise

Insights

✅ Handles linear trends better than Naive
❌ Does not capture seasonal effects or nonlinear trend changes
❗ Still exhibits autocorrelation in residuals, indicating unmodeled structure remains

Seasonal Naive

Assumptions:

Forecasts each period (e.g., month) as equal to the value from the same period in the previous year
Suitable for strong, stable seasonal patterns
Does not account for trend

model_snaive <- wti_ts |>
  model(SNAIVE(WTISPLC))

model_snaive |> gg_tsresiduals()

Residuals: Non-random, strong autocorrelation.
ACF Plot: Autocorrelation evident in ACF
Histogram of Residuals: Histogram not normally distributed

model_snaive |>
  augment() |>  # adds `.resid`, `.innov`, etc.
  features(.innov, ljung_box, lag = 10)

## # A tibble: 1 × 3
##   .model          lb_stat lb_pvalue
##   <chr>             <dbl>     <dbl>
## 1 SNAIVE(WTISPLC)   1421.         0

Residual Diagnostics

Ljung-Box test p-value: 0
✅ Result: ❌ Residuals are not white noise

Insights

✅ Effective baseline when seasonality dominates
❌ Fails to account for trends or evolving seasonal effects
❗ Significant autocorrelation remains → structure is left unmodeled

Seasonal Naive with Drift

model(wti_ts, SNAIVE(WTISPLC ~ drift()))

## # A mable: 1 x 1
##   `SNAIVE(WTISPLC ~ drift())`
##                       <model>
## 1           <SNAIVE w/ drift>

Why This Model Was Not Used

Although both Drift and Seasonal Naive capture important dynamics — trend and seasonality respectively — their assumptions are incompatible:
- Drift assumes a linear progression over time.
- Seasonal Naive assumes a strict repeating seasonal pattern.
The combination implies both a consistent seasonal pattern and a trend that evolves linearly, which is statistically incoherent.

Framework Limitation

The {fable} package does not support this model.
Attempting to fit SNAIVE(WTISPLC ~ drift()) results in an error due to this conceptual conflict.

Conclusion

❌ Seasonal Naive with Drift is excluded from the analysis.
✅ Instead, Drift and SNAIVE are evaluated separately as valid benchmarks for trend and seasonality, respectively.

Benchmark Models (Forecast)

# Split data: use data until end of 2021 as training
train <- wti_ts |> filter(DATE < yearmonth("2022 Jan"))

# Refit benchmark models on training set
models <- train |> 
  model(
    Naive = NAIVE(WTISPLC),
    Drift = RW(WTISPLC ~ drift()),
    SeasonalNaive = SNAIVE(WTISPLC)
  )

# Forecast next 12 months (Jan–Dec 2022)
fc <- models |> forecast(h = "12 months")

# Visualize forecast with actual data for context
fc |> 
  autoplot(wti_ts) +  # <--- overlay with actual full time series
  labs(
    title = "Forecast: Benchmark Models (2022)",
    subtitle = "Naive, Drift, and Seasonal Naive vs Actual WTI Prices",
    y = "WTI Spot Price (WTISPLC)",
    x = "Date"
  ) +
  facet_wrap(~ .model, ncol = 1) +
  theme_minimal()

Advanced Models - Exponential Smoothing

Exponential Smoothing – Simple

Assumptions:

Additive error: Forecast error is independent of level.
No trend: Assumes the series fluctuates around a constant mean.
No seasonality: Ignores any recurring seasonal behavior.

model_ets_simple <- wti_ts |>
  model(ETS(WTISPLC ~ error("A") + trend("N") + season("N")))

model_ets_simple |> gg_tsresiduals()

Residual Time Plot: Fluctuates around zero but shows some clustering → not entirely random.
ACF Plot: Multiple spikes exceed bounds → autocorrelation.
Histogram: Skewed and slightly bumpy → residuals not normally distributed.

model_ets_simple |>
  augment() |>
  features(.innov, ljung_box, lag = 10)

## # A tibble: 1 × 3
##   .model                                                       lb_stat lb_pvalue
##   <chr>                                                          <dbl>     <dbl>
## 1 "ETS(WTISPLC ~ error(\"A\") + trend(\"N\") + season(\"N\"))"    99.5  1.11e-16

Residual Diagnostics

Ljung-Box test p-value: 1.11e-16
✅ Result: ❌ Residuals are not white noise

Insights:

✅ Useful as a simple benchmark when data lacks visible patterns
❌ Ignores trend and seasonality — both are significant in WTI
❗ Residual autocorrelation implies poor model fit
⚠️ Too simplistic for medium- to long-term economic forecasting

Exponential Smoothing – Additive

Assumptions:

Additive error: Noise is added directly to the forecast.
Additive trend: Captures a consistent upward or downward linear movement.
No seasonality: Assumes no repeating seasonal pattern (valid if deseasonalized beforehand).

model_ets_add <- wti_ts |>
  model(ETS(WTISPLC ~ error("A") + trend("A") + season("N")))

model_ets_add |> gg_tsresiduals()

Residual Time Plot: Shows pattern and clustering → unmodeled structure remains.
ACF Plot: Several spikes exceed bounds → autocorrelation.
Histogram: Not perfectly normal → potential bias.

model_ets_add |>
  augment() |>
  features(.innov, ljung_box, lag = 10)

## # A tibble: 1 × 3
##   .model                                                       lb_stat lb_pvalue
##   <chr>                                                          <dbl>     <dbl>
## 1 "ETS(WTISPLC ~ error(\"A\") + trend(\"A\") + season(\"N\"))"    98.1  1.11e-16

Residual Diagnostics

Ljung-Box test p-value: 1.11e-16
✅ Result: ❌ Residuals are not white noise

Insights

✅ Captures linear trend better than baseline models
❌ Misses potential seasonal patterns
❗ Residual autocorrelation suggests unmodeled structure

Exponential Smoothing – Holt’s Linear

Assumptions:

Additive error: Errors are added to forecasts, assuming constant variance.
Additive trend: Captures a consistent linear increase or decrease over time.
No seasonality: Ignores recurring seasonal effects.

model_ets_holt <- wti_ts |>
  model(ETS(WTISPLC ~ error("A") + trend("A") + season("N")))

model_ets_holt |> gg_tsresiduals()

Residual Time Plot: Residuals show gradual shifts → trend is not fully captured.
ACF Plot: Significant autocorrelation at low lags.
Histogram: Deviates from normality — some skewness.

model_ets_holt |>
  augment() |>
  features(.innov, ljung_box, lag = 10)

## # A tibble: 1 × 3
##   .model                                                       lb_stat lb_pvalue
##   <chr>                                                          <dbl>     <dbl>
## 1 "ETS(WTISPLC ~ error(\"A\") + trend(\"A\") + season(\"N\"))"    98.1  1.11e-16

Residual Diagnostics

Ljung-Box test p-value: 1.11e-16
✅ Result: ❌ Residuals are not white noise

Insights:

✅ Effectively models linear trends
❌ Fails to address nonlinear trends or seasonality
❗ Residual autocorrelation suggests missing dynamics in the series
⚠️ May underperform on data with evolving or seasonal structure (e.g., WTI)

Exponential Smoothing – Damped Trend

Assumptions:

Additive error: Errors are added to the forecast.
Damped additive trend: Trend is present but slows (dampens) over time.
No seasonality: Assumes no recurring seasonal component.

model_ets_damped <- wti_ts |>
  model(ETS(WTISPLC ~ error("A") + trend("Ad") + season("N")))

model_ets_damped |> gg_tsresiduals()

Residual Time Plot: Less aggressive drift, but still shows patterns → unmodeled structure remains.
ACF Plot: Lags exceed bounds → residuals are autocorrelated.
Histogram: Irregular shape → residuals not perfectly normal.

model_ets_damped |>
  augment() |>
  features(.innov, ljung_box, lag = 10)

## # A tibble: 1 × 3
##   .model                                                       lb_stat lb_pvalue
##   <chr>                                                          <dbl>     <dbl>
## 1 "ETS(WTISPLC ~ error(\"A\") + trend(\"Ad\") + season(\"N\")…    40.0 0.0000171

Residual Diagnostics

Ljung-Box test p-value: < 0.05
✅ Result: ❌ Residuals are not white noise

Insights:

✅ Better long-term realism by avoiding indefinitely increasing/decreasing trends
✅ Suitable for series with trend that slows down over time
❌ Still misses seasonal patterns
❗ Residual autocorrelation indicates remaining unexplained structure

Exponential Smoothing – Multiplicative

Assumptions:

Multiplicative error: Forecast error scales with the level of the series (i.e., larger values lead to proportionally larger deviations).
Additive trend: Models consistent upward or downward movement.
No seasonality: Ignores repeating seasonal patterns.

model_ets_mult <- wti_ts |>
  model(ETS(WTISPLC ~ error("M") + trend("A") + season("N")))

model_ets_mult |> gg_tsresiduals()

model_ets_mult |>
  augment() |>
  features(.innov, ljung_box, lag = 10)

## # A tibble: 1 × 3
##   .model                                                       lb_stat lb_pvalue
##   <chr>                                                          <dbl>     <dbl>
## 1 "ETS(WTISPLC ~ error(\"M\") + trend(\"A\") + season(\"N\"))"    44.8   2.39e-6

Residual Diagnostics

Ljung-Box test p-value: 2.39e-06
✅ Result: ❌ Residuals are not white noise

Insights:

✅ Handles heteroskedasticity better than additive error models
✅ Captures linear trend component
❌ Still leaves autocorrelation in residuals → indicates missed structure
⚠️ May be improved with seasonal or damped variants

Exponential Smoothing – Seasonal Additive Seasonality

Assumptions:

Additive error: Forecast error is constant regardless of level.
Additive trend: Linear increase or decrease over time.
Additive seasonality: Seasonal effect adds a constant amount (fixed size) for each season.

model_ets_add_season <- wti_ts |>
  model(ETS(WTISPLC ~ error("A") + trend("A") + season("A")))

model_ets_add_season |> gg_tsresiduals()

Residual Time Plot: Less clustering, though some mild drift remains.
ACF Plot: Fewer but still notable autocorrelation spikes.
Histogram: Closer to normal shape but still slightly irregular.

model_ets_add_season |>
  augment() |>
  features(.innov, ljung_box, lag = 10)

## # A tibble: 1 × 3
##   .model                                                       lb_stat lb_pvalue
##   <chr>                                                          <dbl>     <dbl>
## 1 "ETS(WTISPLC ~ error(\"A\") + trend(\"A\") + season(\"A\"))"    89.2  7.66e-15

Residual Diagnostics

Ljung-Box test p-value: < 0.05
✅ Result: ❌ Residuals are not white noise

Insights:

✅ Captures both trend and seasonal structures effectively
✅ Suitable for series with stable seasonal patterns and linear trend
❌ Some autocorrelation remains → unmodeled complexity
❗ May struggle with non-linear trend or changing seasonality

Exponential Smoothing – Multiplicative Seasonality

Assumptions:

Multiplicative error: Forecast uncertainty increases as the level increases.
Additive trend: Captures steady, linear upward/downward progression.
Multiplicative seasonality: Seasonal effects change proportionally with the level of the series (e.g., seasonality amplifies when prices are high).

model_ets_mult_season <- wti_ts |>
  model(ETS(WTISPLC ~ error("M") + trend("A") + season("M")))

model_ets_mult_season |> gg_tsresiduals()

Residual Time Plot: Best structure among ETS models — fewer patterns, relatively stable.
ACF Plot: Minimal spikes; some autocorrelation remains.
Histogram: Closest to bell-shaped → better normality than others.

model_ets_mult_season |>
  augment() |>
  features(.innov, ljung_box, lag = 10)

## # A tibble: 1 × 3
##   .model                                                       lb_stat lb_pvalue
##   <chr>                                                          <dbl>     <dbl>
## 1 "ETS(WTISPLC ~ error(\"M\") + trend(\"A\") + season(\"M\"))"    50.0   2.65e-7

Residual Diagnostics

Ljung-Box test p-value: < 0.05
✅ Result: ❌ Residuals are not white noise

Insights:

✅ Best for data where seasonal variation grows/shrinks with level (e.g., energy prices)
✅ Handles multiplicative dynamics better than additive seasonal models
❌ Residual autocorrelation remains → model still leaves structure unexplained
❗ Strong candidate when volatility varies with trend and seasonal effects

Advanced Models - ETS (Forecast and Accuracy)

# 1. Split data: training before 2022, testing for 2022
train <- wti_ts |> filter(DATE < yearmonth("2022 Jan"))
test_2022 <- wti_ts |> filter_index("2022 Jan" ~ "2022 Dec")

# 2. Fit multiple ETS models on training data
ets_models <- train |> 
  model(
    ETS_Simple      = ETS(WTISPLC ~ error("A") + trend("N") + season("N")),
    ETS_Add         = ETS(WTISPLC ~ error("A") + trend("A") + season("N")),
    ETS_Holt        = ETS(WTISPLC ~ error("A") + trend("A") + season("N")),
    ETS_Damped      = ETS(WTISPLC ~ error("A") + trend("Ad") + season("N")),
    ETS_Mult        = ETS(WTISPLC ~ error("M") + trend("A") + season("N")),
    ETS_Add_Season  = ETS(WTISPLC ~ error("A") + trend("A") + season("A")),
    ETS_Mult_Season = ETS(WTISPLC ~ error("M") + trend("A") + season("M"))
  )

fc_ets <- forecast(ets_models, h = "12 months")

fc_ets |>
  autoplot(wti_ts) +
  labs(
    title = "Forecast: ETS Models",
    subtitle = "Various Models",
    y = "WTI Spot Price (WTISPLC)",
    x = "Date"
  ) +
  facet_wrap(~ .model, ncol = 1, scales = "free_y") +  # Display one model per row
  theme_minimal() +
  theme(
    strip.text = element_text(face = "bold"),  # Model titles
    plot.title = element_text(face = "bold", size = 14)
  )

Auto Regressive Integrated Moving Average (ARIMA)

Assumptions:

Autoregressive terms (AR): Past values influence the current observation
Differencing (I): Removes trends to ensure stationarity
Moving average terms (MA): Past errors influence the current observation
Seasonality: Detected automatically and included if statistically significant

train <- wti_ts |> filter(DATE < yearmonth("2022 Jan"))
test  <- wti_ts |> filter(DATE >= yearmonth("2022 Jan") & DATE <= yearmonth("2022 Dec"))

# Fit multiple manual ARIMA models
models_manual <- train |> model(
  ARIMA_210 = ARIMA(WTISPLC ~ pdq(2,1,0)),
  ARIMA_013 = ARIMA(WTISPLC ~ pdq(0,1,3)),
  ARIMA_111 = ARIMA(WTISPLC ~ pdq(1,1,1)),
  stepwise = ARIMA(WTISPLC),
  search = ARIMA(WTISPLC, stepwise = FALSE)
)

glance(models_manual) |> arrange(AICc) |> select(.model, AICc, BIC)

## # A tibble: 5 × 3
##   .model     AICc   BIC
##   <chr>     <dbl> <dbl>
## 1 search    2786. 2811.
## 2 stepwise  2794. 2819.
## 3 ARIMA_013 2800. 2817.
## 4 ARIMA_210 2801. 2813.
## 5 ARIMA_111 2801. 2813.

models_manual |>
  select(search) |>  # or your best model
  gg_tsresiduals()

glance(models_manual)

## # A tibble: 5 × 8
##   .model    sigma2 log_lik   AIC  AICc   BIC ar_roots  ma_roots  
##   <chr>      <dbl>   <dbl> <dbl> <dbl> <dbl> <list>    <list>    
## 1 ARIMA_210   15.2  -1397. 2801. 2801. 2813. <cpl [2]> <cpl [0]> 
## 2 ARIMA_013   15.1  -1396. 2800. 2800. 2817. <cpl [0]> <cpl [3]> 
## 3 ARIMA_111   15.2  -1397. 2801. 2801. 2813. <cpl [1]> <cpl [1]> 
## 4 stepwise    14.9  -1391. 2793. 2794. 2819. <cpl [2]> <cpl [3]> 
## 5 search      14.7  -1387. 2786. 2786. 2811. <cpl [2]> <cpl [25]>

augment(models_manual) |>
  filter(.model == 'search') |>
  features(.innov, ljung_box, lag = 10, dof = 3)

## # A tibble: 1 × 3
##   .model lb_stat lb_pvalue
##   <chr>    <dbl>     <dbl>
## 1 search    3.74     0.809

forecast_search <- models_manual |> 
  forecast(new_data = test)  # test covers Jan–Dec 2022

forecast_search |>
  autoplot(wti_ts) +  # Plot on full data (train + test)
  autolayer(test, WTISPLC, color = "red") +
  labs(
    title = "Forecast: ARIMA Models",
    subtitle = "Comparison Across Selected Models",
    y = "WTI Spot Price (WTISPLC)",
    x = "Date"
  ) +
  facet_wrap(~ .model, ncol = 1, scales = "free_y") +  # One model per row
  theme_minimal() +
  theme(
    strip.text = element_text(face = "bold"),
    plot.title = element_text(face = "bold", size = 14)
  )

accuracy(forecast_search, test)

## # A tibble: 5 × 10
##   .model    .type    ME  RMSE   MAE   MPE  MAPE  MASE RMSSE  ACF1
##   <chr>     <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ARIMA_013 Test   27.3  29.8  27.3  27.7  27.7   NaN   NaN 0.583
## 2 ARIMA_111 Test   27.1  29.6  27.1  27.5  27.5   NaN   NaN 0.580
## 3 ARIMA_210 Test   27.0  29.5  27.0  27.4  27.4   NaN   NaN 0.581
## 4 search    Test   28.0  30.1  28.0  28.5  28.5   NaN   NaN 0.556
## 5 stepwise  Test   31.6  33.8  31.6  32.3  32.3   NaN   NaN 0.557

models_sarima <- train |> model(
  SARIMA_011_011 = ARIMA(WTISPLC ~ pdq(0,1,1) + PDQ(0,1,1)),
  SARIMA_111_011 = ARIMA(WTISPLC ~ pdq(1,1,1) + PDQ(0,1,1)),
  SARIMA_210_011 = ARIMA(WTISPLC ~ pdq(2,1,0) + PDQ(0,1,1)),
  SARIMA_stepwise = ARIMA(WTISPLC),
  SARIMA_search = ARIMA(WTISPLC, stepwise = FALSE)
)

glance(models_sarima) |> arrange(AICc) |> select(.model, AICc, BIC)

## # A tibble: 5 × 3
##   .model           AICc   BIC
##   <chr>           <dbl> <dbl>
## 1 SARIMA_210_011  2782. 2798.
## 2 SARIMA_111_011  2782. 2799.
## 3 SARIMA_search   2786. 2811.
## 4 SARIMA_011_011  2788. 2801.
## 5 SARIMA_stepwise 2794. 2819.

augment(models_sarima) |>
  filter(.model == 'SARIMA_search') |>
  features(.innov, ljung_box, lag = 10, dof = 3)

## # A tibble: 1 × 3
##   .model        lb_stat lb_pvalue
##   <chr>           <dbl>     <dbl>
## 1 SARIMA_search    3.74     0.809

forecast_sarima <- models_sarima |> 
  forecast(new_data = test)

forecast_sarima |> 
  autoplot(train, level = NULL) + 
  autolayer(test, WTISPLC, color = "red") +
  labs(title = "SARIMA Forecast vs Actuals (2022)",
       y = "Price", x = "Date")

accuracy(forecast_sarima, test)

## # A tibble: 5 × 10
##   .model          .type    ME  RMSE   MAE   MPE  MAPE  MASE RMSSE  ACF1
##   <chr>           <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 SARIMA_011_011  Test   22.3  25.0  22.3  22.4  22.4   NaN   NaN 0.602
## 2 SARIMA_111_011  Test   24.0  26.6  24.0  24.2  24.2   NaN   NaN 0.595
## 3 SARIMA_210_011  Test   23.9  26.5  23.9  24.1  24.1   NaN   NaN 0.596
## 4 SARIMA_search   Test   28.0  30.1  28.0  28.5  28.5   NaN   NaN 0.556
## 5 SARIMA_stepwise Test   31.6  33.8  31.6  32.3  32.3   NaN   NaN 0.557

AutoRegressive Distributed Lag (ARDL)

Assumptions:

Autoregressive terms (AR): Current values are influenced by their own past values.
Distributed lag terms (DL): Past values of other explanatory variables can affect the current value.
Can accommodate a mix of stationary (I(0)) and non-stationary (I(1)) variables.
Stationarity not strictly required, but stability in variance is still beneficial.

Comparison of all models

# Create unified model list (individual ETS and ARIMA models)
model_list <- list(
  Naive             = model_naive,
  Drift             = model_drift,
  Seasonal_Naive    = model_snaive,
  
  # ETS models
  ETS_Simple        = model_ets_simple,
  ETS_Add           = model_ets_add,
  ETS_Holt          = model_ets_holt,
  ETS_Damped        = model_ets_damped,
  ETS_Add_Season    = model_ets_add_season,
  ETS_Mult          = model_ets_mult,
  ETS_Mult_Season   = model_ets_mult_season,
  
  # ARIMA models (extracted from models_manual collection)
  ARIMA_210         = models_manual |> select(ARIMA_210),
  ARIMA_013         = models_manual |> select(ARIMA_013),
  ARIMA_111         = models_manual |> select(ARIMA_111),
  ARIMA_stepwise    = models_manual |> select(stepwise),
  ARIMA_search      = models_manual |> select(search)
)

# Safe Ljung-Box function wrapper
safe_ljung <- safely(function(m) {
  m |> augment() |> features(.innov, ljung_box, lag = 10, dof = 3)
})

# Compute p-values for each model
ljung_pvalues <- imap_dfr(model_list, function(mdl, name) {
  res <- safe_ljung(mdl)
  if (!is.null(res$result)) {
    tibble(
      Model = name,
      `Ljung-Box p-value` = round(res$result$lb_pvalue, 4)
    )
  } else {
    tibble(
      Model = name,
      `Ljung-Box p-value` = NA
    )
  }
})

# Print the summary table
print(ljung_pvalues)

## # A tibble: 15 × 2
##    Model           `Ljung-Box p-value`
##    <chr>                         <dbl>
##  1 Naive                         0    
##  2 Drift                         0    
##  3 Seasonal_Naive                0    
##  4 ETS_Simple                    0    
##  5 ETS_Add                       0    
##  6 ETS_Holt                      0    
##  7 ETS_Damped                    0    
##  8 ETS_Add_Season                0    
##  9 ETS_Mult                      0    
## 10 ETS_Mult_Season               0    
## 11 ARIMA_210                     0.125
## 12 ARIMA_013                     0.266
## 13 ARIMA_111                     0.127
## 14 ARIMA_stepwise                0.901
## 15 ARIMA_search                  0.809

❗ Most models fail the white noise test → residuals are still autocorrelated.

✅ Only ARIMA passed → best at capturing underlying structure statistically.

# Split into train and test
train <- wti_ts |> filter(DATE < yearmonth("2022 Jan"))
test  <- wti_ts |> filter(DATE >= yearmonth("2022 Jan") & DATE <= yearmonth("2022 Dec"))

# Fit multiple ETS and ARIMA models
models_all <- train |> model(
  Naive           = NAIVE(WTISPLC),
  Drift           = RW(WTISPLC ~ drift()),
  SeasonalNaive   = SNAIVE(WTISPLC),
  
  # ETS Models
  ETS_Simple      = ETS(WTISPLC ~ error("A") + trend("N") + season("N")),
  ETS_Add         = ETS(WTISPLC ~ error("A") + trend("A") + season("N")),
  ETS_Holt        = ETS(WTISPLC ~ error("A") + trend("A") + season("N")),
  ETS_Damped      = ETS(WTISPLC ~ error("A") + trend("Ad") + season("N")),
  ETS_Add_Season  = ETS(WTISPLC ~ error("A") + trend("A") + season("A")),
  ETS_Mult        = ETS(WTISPLC ~ error("M") + trend("A") + season("N")),
  ETS_Mult_Season = ETS(WTISPLC ~ error("M") + trend("A") + season("M")),
  
  # ARIMA Models
  ARIMA_stepwise  = ARIMA(WTISPLC),
  ARIMA_search    = ARIMA(WTISPLC, stepwise = FALSE),
  ARIMA_manual_210 = ARIMA(WTISPLC ~ pdq(2,1,0)),
  ARIMA_manual_013 = ARIMA(WTISPLC ~ pdq(0,1,3)),
  ARIMA_manual_111 = ARIMA(WTISPLC ~ pdq(1,1,1))
)

# Forecast next 12 months (Jan–Dec 2022)
fc_all <- forecast(models_all, new_data = test)

# Accuracy Table: Sorted by MAPE on Test Set
# Compute accuracy against actuals (test set)
accuracy_all <- accuracy(fc_all, test)

# Tabulate and sort results: focus on MAPE
accuracy_all |>
  select(.model, .type, ME, RMSE, MAE, MPE, MAPE, MASE, RMSSE, ACF1) |>
  filter(.type == "Test") |>
  arrange(MAPE)

## # A tibble: 15 × 10
##    .model           .type    ME  RMSE   MAE   MPE  MAPE  MASE RMSSE  ACF1
##    <chr>            <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
##  1 ETS_Mult_Season  Test   18.0  21.3  18.0  17.9  17.9   NaN   NaN 0.653
##  2 ETS_Add_Season   Test   21.2  23.9  21.2  21.2  21.2   NaN   NaN 0.587
##  3 ETS_Mult         Test   21.5  24.7  21.5  21.4  21.4   NaN   NaN 0.607
##  4 ETS_Add          Test   22.5  25.5  22.5  22.6  22.6   NaN   NaN 0.594
##  5 ETS_Holt         Test   22.5  25.5  22.5  22.6  22.6   NaN   NaN 0.594
##  6 Drift            Test   22.6  25.5  22.6  22.6  22.6   NaN   NaN 0.593
##  7 ETS_Simple       Test   23.1  25.9  23.1  23.2  23.2   NaN   NaN 0.587
##  8 Naive            Test   23.1  25.9  23.1  23.2  23.2   NaN   NaN 0.587
##  9 SeasonalNaive    Test   26.8  30.8  26.8  27.0  27.0   NaN   NaN 0.787
## 10 ARIMA_manual_210 Test   27.0  29.5  27.0  27.4  27.4   NaN   NaN 0.581
## 11 ARIMA_manual_111 Test   27.1  29.6  27.1  27.5  27.5   NaN   NaN 0.580
## 12 ARIMA_manual_013 Test   27.3  29.8  27.3  27.7  27.7   NaN   NaN 0.583
## 13 ARIMA_search     Test   28.0  30.1  28.0  28.5  28.5   NaN   NaN 0.556
## 14 ETS_Damped       Test   30.8  32.9  30.8  31.5  31.5   NaN   NaN 0.535
## 15 ARIMA_stepwise   Test   31.6  33.8  31.6  32.3  32.3   NaN   NaN 0.557

# Note: MASE and RMSSE are NaN since baseline scaling may not apply to seasonal models

print(accuracy_all)

## # A tibble: 15 × 10
##    .model           .type    ME  RMSE   MAE   MPE  MAPE  MASE RMSSE  ACF1
##    <chr>            <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
##  1 ARIMA_manual_013 Test   27.3  29.8  27.3  27.7  27.7   NaN   NaN 0.583
##  2 ARIMA_manual_111 Test   27.1  29.6  27.1  27.5  27.5   NaN   NaN 0.580
##  3 ARIMA_manual_210 Test   27.0  29.5  27.0  27.4  27.4   NaN   NaN 0.581
##  4 ARIMA_search     Test   28.0  30.1  28.0  28.5  28.5   NaN   NaN 0.556
##  5 ARIMA_stepwise   Test   31.6  33.8  31.6  32.3  32.3   NaN   NaN 0.557
##  6 Drift            Test   22.6  25.5  22.6  22.6  22.6   NaN   NaN 0.593
##  7 ETS_Add          Test   22.5  25.5  22.5  22.6  22.6   NaN   NaN 0.594
##  8 ETS_Add_Season   Test   21.2  23.9  21.2  21.2  21.2   NaN   NaN 0.587
##  9 ETS_Damped       Test   30.8  32.9  30.8  31.5  31.5   NaN   NaN 0.535
## 10 ETS_Holt         Test   22.5  25.5  22.5  22.6  22.6   NaN   NaN 0.594
## 11 ETS_Mult         Test   21.5  24.7  21.5  21.4  21.4   NaN   NaN 0.607
## 12 ETS_Mult_Season  Test   18.0  21.3  18.0  17.9  17.9   NaN   NaN 0.653
## 13 ETS_Simple       Test   23.1  25.9  23.1  23.2  23.2   NaN   NaN 0.587
## 14 Naive            Test   23.1  25.9  23.1  23.2  23.2   NaN   NaN 0.587
## 15 SeasonalNaive    Test   26.8  30.8  26.8  27.0  27.0   NaN   NaN 0.787

ETS_Mult_Season wins in forecast accuracy across all error metrics.
Despite its strong residual diagnostics, ARIMA fails to capture seasonal variation, resulting in higher forecast error compared to ETS_Mult_Season.

Summary Insights

✅ Best Forecast Accuracy: ETS(WTISPLC ~ error(“M”) + trend(“A”) + season(“M”))
- Captures both multiplicative seasonal effects and trend.
- Performs well despite some autocorrelation.
✅ Best Statistical Fit: ARIMA(WTISPLC)
- Residuals are white noise.
- But lags behind in forecast accuracy due to less emphasis on seasonal structure.
❗ All Other Models: Either too simple (Naive, Drift) or miss key dynamics (e.g., trend dampening fails without seasonality).

Final Recommendation

➡️ Use ARIMA for forecasting WTI prices over the next 12 months.

# Final model fitted on entire dataset
model_final <- wti_ts |>
  model(ARIMA_search    = ARIMA(WTISPLC, stepwise = FALSE))

# Forecast next 12 months
forecast_final <- forecast(model_final, h = "12 months")

# Convert to tibble and extract forecasted values
forecast_table <- forecast_final |> 
  as_tibble() |> 
  select(DATE, .mean) |> 
  rename(
    `Forecast Month` = DATE,
    `Forecasted WTI Price` = .mean
  )

# Print the table
forecast_table

## # A tibble: 12 × 2
##    `Forecast Month` `Forecasted WTI Price`
##               <mth>                  <dbl>
##  1         2023 Jan                   73.9
##  2         2023 Feb                   72.4
##  3         2023 Mar                   72.7
##  4         2023 Apr                   72.4
##  5         2023 May                   73.2
##  6         2023 Jun                   73.6
##  7         2023 Jul                   73.2
##  8         2023 Aug                   73.7
##  9         2023 Sep                   73.0
## 10         2023 Oct                   72.4
## 11         2023 Nov                   72.6
## 12         2023 Dec                   73.2

forecast_final |>
  autoplot(wti_ts) +
  labs(
    title = "Final 12-Month Forecast for WTI Prices",
    y = "WTISPLC",
    x = "Date"
  ) +
  theme_minimal()

forecast_table |>
  mutate(Quarter = quarter(`Forecast Month`, with_year = TRUE)) |>
  group_by(Quarter) |>
  summarise(`Avg Forecast Price` = mean(`Forecasted WTI Price`)) |>
  arrange(Quarter)

## # A tibble: 4 × 2
##   Quarter `Avg Forecast Price`
##     <dbl>                <dbl>
## 1   2023.                 73.0
## 2   2023.                 73.1
## 3   2023.                 73.3
## 4   2023.                 72.7

Executive Summary

Best Forecasting Model: ARIMA (Auto-Selected)

Forecasted Trend: Jan–Dec 2023

Q1 Dip: $73.0 average
Q3 Peak: $73.3 average
Q4 Dip: $72.7 average
Gradual upward trend throughout the year

Business actions:

Lock fuel prices/contracts early to avoid Q3–Q4 surge
Plan for fuel and import inflation in budgeting
Adjust cost escalation clauses in supplier agreements

Forecasting Analytics - Data Science Assessment

Toby G. Mejia

June 27, 2025

Objective:

Exploratory Data Analysis

Forecasting

Why Forecast WTI Spot Prices?

Comparison of all models

Final Recommendation

Executive Summary