Time Series forecasting for Commodities
Stephen Situ
2022-11-18
This is a project that aims to conduct time series forecasting for Brent Oil, Crude Oil WTI, Heating Oil, Natural Gas. We do some preliminary analysis and identify strong correlation between Brent Oil, Crude Oil, and Heating Oil. Afterwards, we try to identify seasonal trends and conduct time series modelling with the seasonal naive method, ets method, and arima method.
Original dataset can be found here: https://www.kaggle.com/datasets/prasertk/historical-daily-oil-and-natural-gas-prices
Install Packages
install.packages('tidyverse', repos = "http://cran.us.r-project.org")## Installing package into 'C:/Users/Steve/AppData/Local/R/win-library/4.2'
## (as 'lib' is unspecified)## package 'tidyverse' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
## C:\Users\Steve\AppData\Local\Temp\RtmpcnGyoc\downloaded_packagesinstall.packages('ggfortify', repos = "http://cran.us.r-project.org")## Installing package into 'C:/Users/Steve/AppData/Local/R/win-library/4.2'
## (as 'lib' is unspecified)## package 'ggfortify' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
## C:\Users\Steve\AppData\Local\Temp\RtmpcnGyoc\downloaded_packagesinstall.packages('forecast', repos = "http://cran.us.r-project.org")## Installing package into 'C:/Users/Steve/AppData/Local/R/win-library/4.2'
## (as 'lib' is unspecified)## package 'forecast' successfully unpacked and MD5 sums checked## Warning: cannot remove prior installation of package 'forecast'## Warning in file.copy(savedcopy, lib, recursive = TRUE):
## problem copying C:\Users\Steve\AppData\Local\R\win-
## library\4.2\00LOCK\forecast\libs\x64\forecast.dll to C:
## \Users\Steve\AppData\Local\R\win-library\4.2\forecast\libs\x64\forecast.dll:
## Permission denied## Warning: restored 'forecast'##
## The downloaded binary packages are in
## C:\Users\Steve\AppData\Local\Temp\RtmpcnGyoc\downloaded_packagesinstall.packages('lubridate', repos = "http://cran.us.r-project.org")## Installing package into 'C:/Users/Steve/AppData/Local/R/win-library/4.2'
## (as 'lib' is unspecified)## package 'lubridate' successfully unpacked and MD5 sums checked## Warning: cannot remove prior installation of package 'lubridate'## Warning in file.copy(savedcopy, lib, recursive = TRUE):
## problem copying C:\Users\Steve\AppData\Local\R\win-
## library\4.2\00LOCK\lubridate\libs\x64\lubridate.dll to C:
## \Users\Steve\AppData\Local\R\win-library\4.2\lubridate\libs\x64\lubridate.dll:
## Permission denied## Warning: restored 'lubridate'##
## The downloaded binary packages are in
## C:\Users\Steve\AppData\Local\Temp\RtmpcnGyoc\downloaded_packagesinstall.packages("corrplot", repos = "http://cran.us.r-project.org")## Installing package into 'C:/Users/Steve/AppData/Local/R/win-library/4.2'
## (as 'lib' is unspecified)## package 'corrplot' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
## C:\Users\Steve\AppData\Local\Temp\RtmpcnGyoc\downloaded_packagesinstall.packages("dplyr", repos = "http://cran.us.r-project.org")## Installing package into 'C:/Users/Steve/AppData/Local/R/win-library/4.2'
## (as 'lib' is unspecified)## package 'dplyr' successfully unpacked and MD5 sums checked## Warning: cannot remove prior installation of package 'dplyr'## Warning in file.copy(savedcopy, lib, recursive = TRUE): problem copying C:
## \Users\Steve\AppData\Local\R\win-library\4.2\00LOCK\dplyr\libs\x64\dplyr.dll
## to C:\Users\Steve\AppData\Local\R\win-library\4.2\dplyr\libs\x64\dplyr.dll:
## Permission denied## Warning: restored 'dplyr'##
## The downloaded binary packages are in
## C:\Users\Steve\AppData\Local\Temp\RtmpcnGyoc\downloaded_packageslibrary(tidyverse)## ── Attaching packages
## ───────────────────────────────────────
## tidyverse 1.3.2 ──## ✔ ggplot2 3.4.0 ✔ purrr 0.3.5
## ✔ tibble 3.1.8 ✔ dplyr 1.0.10
## ✔ tidyr 1.2.1 ✔ stringr 1.4.1
## ✔ readr 2.1.3 ✔ forcats 0.5.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()library(ggfortify)
library(forecast)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
## Registered S3 methods overwritten by 'forecast':
## method from
## autoplot.Arima ggfortify
## autoplot.acf ggfortify
## autoplot.ar ggfortify
## autoplot.bats ggfortify
## autoplot.decomposed.ts ggfortify
## autoplot.ets ggfortify
## autoplot.forecast ggfortify
## autoplot.stl ggfortify
## autoplot.ts ggfortify
## fitted.ar ggfortify
## fortify.ts ggfortify
## residuals.ar ggfortifylibrary(lubridate)##
## Attaching package: 'lubridate'
##
## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, unionlibrary(corrplot)## corrplot 0.92 loadedlibrary(dplyr)Read csv
energy_df <- read_csv('oil and gas.csv')## Rows: 23024 Columns: 8
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (2): Symbol, Currency
## dbl (5): Open, High, Low, Close, Volume
## date (1): Date
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.structure
summary(energy_df)## Symbol Date Open High
## Length:23024 Min. :2000-01-04 Min. :-14.000 Min. : 0.5085
## Class :character 1st Qu.:2005-10-03 1st Qu.: 2.741 1st Qu.: 2.7850
## Mode :character Median :2011-06-03 Median : 14.364 Median : 15.5325
## Mean :2011-05-09 Mean : 33.361 Mean : 33.8497
## 3rd Qu.:2016-12-22 3rd Qu.: 60.653 3rd Qu.: 61.5000
## Max. :2022-06-17 Max. :146.300 Max. :147.5000
## Low Close Volume Currency
## Min. :-16.740 Min. : 0.4999 Min. : 0 Length:23024
## 1st Qu.: 2.696 1st Qu.: 2.7409 1st Qu.: 32913 Class :character
## Median : 13.851 Median : 16.4140 Median : 71121 Mode :character
## Mean : 32.849 Mean : 33.2554 Mean : 128810
## 3rd Qu.: 59.782 3rd Qu.: 58.9300 3rd Qu.: 178621
## Max. :144.250 Max. :146.0800 Max. :1404916Create column for average price (High+Low)/2
energy_df1 <- mutate(energy_df, average=(High + Low)/2)Create Year, Month column
energy_df1$month <- format(as.Date(energy_df1$Date, format="%d/%m/%Y"),"%m")
energy_df1$year <- format(as.Date(energy_df1$Date, format="%d/%m/%Y"),"%Y")Group by Symbol, year, month
energy_df2 <- energy_df1 %>% group_by(Symbol,year, month) %>% summarise(average_monthly_price = mean(average))
energy_df2 <- energy_df1 %>% group_by(Symbol,year, month) %>% summarise(average_monthly_price = mean(average))## `summarise()` has grouped output by 'Symbol', 'year'. You can override using
## the `.groups` argument.Split into 4 dataframes for each commodity
brent_df <- energy_df2 %>% filter(Symbol == 'Brent Oil')
crude_df <- energy_df2 %>% filter(Symbol == 'Crude Oil WTI')
natural_gas_df <- energy_df2 %>% filter(Symbol == 'Natural Gas')
heating_oil_df <- energy_df2 %>% filter(Symbol == 'Heating Oil')Create time series object
brent_ts <- ts(brent_df[,4],start=c(2000,1),frequency = 12)
crude_df_ts <- ts(crude_df[,4],start=c(2000,1),frequency = 12)
natural_gas_ts <- ts(natural_gas_df[,4],start=c(2000,1),frequency = 12)
heating_oil_ts <- ts(heating_oil_df[,4],start=c(2000,1),frequency = 12)Create an auto plot of each
All four commodities are notoriously volatile and seem to have high correlation upon first glance
We will try identify important/significant financial/economic events
autoplot(brent_ts) + ggtitle('US Brent Oil Price') + theme(plot.title = element_text(hjust = 0.5)) +
labs(x="Year",y="Price (USD)") + annotate("rect", xmin = 2008, xmax = 2009, ymin = 35, ymax = 140, alpha = .2) +
annotate("rect", xmin = 2014, xmax = 2016, ymin = 30, ymax = 110, alpha = .2) +
annotate("rect", xmin = 2020, xmax = 2023, ymin = 20, ymax = 120, alpha = .2) +
annotate("text", x = 2003, y = 75, label = "Structural Bull Market",color="Blue") +
annotate("text", x = 2008, y = 30, label = "Financial Crisis",color="Red") +
annotate("text", x = 2003, y = 75, label = "Structural Bull Market",color="Blue") +
annotate("text", x = 2012, y = 130, label = "Partial Recovery",color="Blue") +
annotate("text", x = 2015, y = 25, label = "2014 Oil Crash",color="Red") +
annotate("text", x = 2018, y = 90, label = "Structural Bear Market",color="Blue") +
annotate("text", x = 2021, y = 15, label = "Covid 19 Crash",color="Red") +
annotate("text", x = 2021, y = 130, label = "Covid 19 Recovery",color="Blue") 
autoplot(crude_df_ts) + ggtitle('US crude Oil WTI Price') + theme(plot.title = element_text(hjust = 0.5)) +
labs(x="Year",y="Price (USD)") + annotate("rect", xmin = 2008, xmax = 2009, ymin = 35, ymax = 140, alpha = .2) +
annotate("rect", xmin = 2014, xmax = 2016, ymin = 30, ymax = 110, alpha = .2) +
annotate("rect", xmin = 2020, xmax = 2023, ymin = 20, ymax = 120, alpha = .2) +
annotate("text", x = 2003, y = 75, label = "Structural Bull Market",color="Blue") +
annotate("text", x = 2008, y = 30, label = "Financial Crisis",color="Red") +
annotate("text", x = 2003, y = 75, label = "Structural Bull Market",color="Blue") +
annotate("text", x = 2012, y = 120, label = "Partial Recovery",color="Blue") +
annotate("text", x = 2015, y = 25, label = "2014 Oil Crash",color="Red") +
annotate("text", x = 2018, y = 75, label = "Structural Bear Market",color="Blue") +
annotate("text", x = 2021, y = 15, label = "Covid 19 Crash",color="Red") +
annotate("text", x = 2021, y = 130, label = "Covid 19 Recovery",color="Blue") 
autoplot(heating_oil_ts) + ggtitle('US Heating Oil Price') + theme(plot.title = element_text(hjust = 0.5)) + theme(plot.title = element_text(hjust = 0.5)) +
labs(x="Year",y="Price (USD)") + labs(x="Year",y="Price (USD)") + annotate("rect", xmin = 2008, xmax = 2009, ymin = 1.25, ymax = 4, alpha = .2) +
annotate("rect", xmin = 2014, xmax = 2016, ymin = 1, ymax = 3.25, alpha = .2) +
annotate("rect", xmin = 2020, xmax = 2023, ymin = 0.75, ymax = 4.30, alpha = .2) +
annotate("text", x = 2003, y = 3, label = "Structural Bull Market",color="Blue") +
annotate("text", x = 2008, y = 1, label = "Financial Crisis",color="Red") +
annotate("text", x = 2012, y = 3.5, label = "Partial Recovery",color="Blue") +
annotate("text", x = 2015, y = 0.75, label = "2014 Oil Crash",color="Red") +
annotate("text", x = 2018, y = 2.5, label = "Structural Bear Market",color="Blue") +
annotate("text", x = 2020, y = 0.50, label = "Covid 19 Crash",color="Red") +
annotate("text", x = 2021, y = 4.6, label = "Covid 19 Recovery",color="Blue") 
autoplot(natural_gas_ts) + ggtitle('US Natural Gas Price') + theme(plot.title = element_text(hjust = 0.5)) +
labs(x="Year",y="Price (USD)") + labs(x="Year",y="Price (USD)") + annotate("rect", xmin = 2008, xmax = 2009, ymin = 3, ymax = 13, alpha = .2) +
annotate("rect", xmin = 2020, xmax = 2023, ymin = 1, ymax = 9, alpha = .2) +
annotate("text", x = 2003, y = 10, label = "Structural Bull Market",color="Blue") +
annotate("text", x = 2008, y = 2, label = "Financial Crisis",color="Red") +
annotate("text", x = 2014, y = 7, label = "Structural Bear Market",color="Blue") +
annotate("text", x = 2020, y = 0.5, label = "Covid 19 Crash",color="Red") +
annotate("text", x = 2021, y = 10, label = "Covid 19 Recovery",color="Blue") 
Building a Correlation Matrix and Visualization
We see a strong correlation between brent,crude, and heating oil
corr_df <- data.frame(matrix(NA,nrow = 270,ncol = 0))
corr_df$crude_price = crude_df$average_monthly_price
corr_df$brent_price = brent_df$average_monthly_price
corr_df$heatingoil_price = heating_oil_df$average_monthly_price
corr_df$naturalgas_price = natural_gas_df$average_monthly_price
corr <- round(cor(corr_df), digits = 2)
corr## crude_price brent_price heatingoil_price naturalgas_price
## crude_price 1.00 0.98 0.97 0.25
## brent_price 0.98 1.00 0.99 0.14
## heatingoil_price 0.97 0.99 1.00 0.17
## naturalgas_price 0.25 0.14 0.17 1.00 corrplot(corr, method="circle")
# De trend data with, take first difference
brent_diff <- diff(brent_ts)
crude_diff <- diff(crude_df_ts)
ng_diff <- diff(natural_gas_ts)
heating_oil_diff <- diff(heating_oil_ts)Plot the differences and Look for seasonality, Very hard to see seasonal patterns
ggseasonplot(brent_diff)
ggseasonplot(crude_diff)
ggseasonplot(ng_diff)
ggseasonplot(heating_oil_diff)
Seasonal sub series plot allows us to visualize more
Brent, Crude, and Heating Oil drop off during Aug-Dec, and rise during Jan-Jul
Natural Gas rise during Sep-Dec + Apr-Jun, they fall during Jan-Mar + Jul-Aug
ggsubseriesplot(brent_diff)
ggsubseriesplot(crude_diff)
ggsubseriesplot(ng_diff)
ggsubseriesplot(heating_oil_diff)
Use a benchmark seasonal naive method: y_t = y_{t-s} + e_t
Residual SD = 7.8673
brent_s <- snaive(brent_diff)
print(summary(brent_s))##
## Forecast method: Seasonal naive method
##
## Model Information:
## Call: snaive(y = brent_diff)
##
## Residual sd: 7.8673
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.1778596 7.867284 5.583062 253.8256 421.7488 1 0.4107349
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jul 2022 0.8493182 -9.233011 10.931648 -14.570274 16.268911
## Aug 2022 -3.7218182 -13.804148 6.360511 -19.141411 11.697774
## Sep 2022 4.2943182 -5.788011 14.376648 -11.125274 19.713911
## Oct 2022 8.8435498 -1.238780 18.925879 -6.576043 24.263142
## Nov 2022 -2.3999134 -12.482243 7.682416 -17.819506 13.019679
## Dec 2022 -6.4531818 -16.535511 3.629148 -21.872774 8.966411
## Jan 2023 10.5538095 0.471480 20.636139 -4.865783 25.973402
## Feb 2023 8.8974405 -1.184889 18.979770 -6.522152 24.317033
## Mar 2023 18.6754891 8.593160 28.757819 3.255897 34.095082
## Apr 2023 -7.1955487 -17.277878 2.886781 -22.615141 8.224044
## May 2023 5.0229004 -5.059429 15.105230 -10.396692 20.442493
## Jun 2023 8.9197552 -1.162574 19.002085 -6.499837 24.339348
## Jul 2023 0.8493182 -13.409249 15.107885 -20.957279 22.655915
## Aug 2023 -3.7218182 -17.980385 10.536749 -25.528415 18.084779
## Sep 2023 4.2943182 -9.964249 18.552885 -17.512279 26.100915
## Oct 2023 8.8435498 -5.415017 23.102117 -12.963047 30.650146
## Nov 2023 -2.3999134 -16.658481 11.858654 -24.206510 19.406683
## Dec 2023 -6.4531818 -20.711749 7.805385 -28.259779 15.353415
## Jan 2024 10.5538095 -3.704758 24.812377 -11.252787 32.360406
## Feb 2024 8.8974405 -5.361127 23.156008 -12.909156 30.704037
## Mar 2024 18.6754891 4.416922 32.934056 -3.131108 40.482086
## Apr 2024 -7.1955487 -21.454116 7.063019 -29.002145 14.611048
## May 2024 5.0229004 -9.235667 19.281468 -16.783696 26.829497
## Jun 2024 8.9197552 -5.338812 23.178322 -12.886841 30.726352checkresiduals(brent_s)
##
## Ljung-Box test
##
## data: Residuals from Seasonal naive method
## Q* = 147.86, df = 24, p-value < 2.2e-16
##
## Model df: 0. Total lags used: 24Use ETS Model
Use ETS(M,A,N) with Residual SD = 0.0867
brent_ets <- ets(brent_ts)
print(summary(brent_ets))## ETS(M,A,N)
##
## Call:
## ets(y = brent_ts)
##
## Smoothing parameters:
## alpha = 0.9999
## beta = 1e-04
##
## Initial states:
## l = 24.3194
## b = 0.61
##
## sigma: 0.0867
##
## AIC AICc BIC
## 2394.773 2395.000 2412.765
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -0.2532313 5.655817 4.000612 -1.02114 6.735565 0.2237797 0.3866065checkresiduals(brent_ets)
##
## Ljung-Box test
##
## data: Residuals from ETS(M,A,N)
## Q* = 44.192, df = 20, p-value = 0.001418
##
## Model df: 4. Total lags used: 24Use Arima model on stationary data
Use ARIMA(0,1,2)(2,1,0)[12] with residual SD = 6.1212
brent_arima <- auto.arima(brent_ts,d=1,D=1,stepwise=FALSE,approximation=FALSE,trace=TRUE)##
## ARIMA(0,1,0)(0,1,0)[12] : 1791.584
## ARIMA(0,1,0)(0,1,1)[12] : Inf
## ARIMA(0,1,0)(0,1,2)[12] : Inf
## ARIMA(0,1,0)(1,1,0)[12] : 1745.896
## ARIMA(0,1,0)(1,1,1)[12] : Inf
## ARIMA(0,1,0)(1,1,2)[12] : Inf
## ARIMA(0,1,0)(2,1,0)[12] : 1711.648
## ARIMA(0,1,0)(2,1,1)[12] : Inf
## ARIMA(0,1,0)(2,1,2)[12] : Inf
## ARIMA(0,1,1)(0,1,0)[12] : 1755.492
## ARIMA(0,1,1)(0,1,1)[12] : Inf
## ARIMA(0,1,1)(0,1,2)[12] : Inf
## ARIMA(0,1,1)(1,1,0)[12] : 1713.704
## ARIMA(0,1,1)(1,1,1)[12] : Inf
## ARIMA(0,1,1)(1,1,2)[12] : Inf
## ARIMA(0,1,1)(2,1,0)[12] : 1680.386
## ARIMA(0,1,1)(2,1,1)[12] : Inf
## ARIMA(0,1,1)(2,1,2)[12] : Inf
## ARIMA(0,1,2)(0,1,0)[12] : 1743.489
## ARIMA(0,1,2)(0,1,1)[12] : Inf
## ARIMA(0,1,2)(0,1,2)[12] : Inf
## ARIMA(0,1,2)(1,1,0)[12] : 1705.709
## ARIMA(0,1,2)(1,1,1)[12] : Inf
## ARIMA(0,1,2)(1,1,2)[12] : Inf
## ARIMA(0,1,2)(2,1,0)[12] : 1672.69
## ARIMA(0,1,2)(2,1,1)[12] : Inf
## ARIMA(0,1,3)(0,1,0)[12] : 1745.228
## ARIMA(0,1,3)(0,1,1)[12] : Inf
## ARIMA(0,1,3)(0,1,2)[12] : Inf
## ARIMA(0,1,3)(1,1,0)[12] : 1707.242
## ARIMA(0,1,3)(1,1,1)[12] : Inf
## ARIMA(0,1,3)(2,1,0)[12] : 1674.323
## ARIMA(0,1,4)(0,1,0)[12] : 1747.307
## ARIMA(0,1,4)(0,1,1)[12] : Inf
## ARIMA(0,1,4)(1,1,0)[12] : 1709.308
## ARIMA(0,1,5)(0,1,0)[12] : 1749.385
## ARIMA(1,1,0)(0,1,0)[12] : 1746.202
## ARIMA(1,1,0)(0,1,1)[12] : Inf
## ARIMA(1,1,0)(0,1,2)[12] : Inf
## ARIMA(1,1,0)(1,1,0)[12] : 1707.05
## ARIMA(1,1,0)(1,1,1)[12] : Inf
## ARIMA(1,1,0)(1,1,2)[12] : Inf
## ARIMA(1,1,0)(2,1,0)[12] : 1673.884
## ARIMA(1,1,0)(2,1,1)[12] : Inf
## ARIMA(1,1,0)(2,1,2)[12] : Inf
## ARIMA(1,1,1)(0,1,0)[12] : 1748.25
## ARIMA(1,1,1)(0,1,1)[12] : Inf
## ARIMA(1,1,1)(0,1,2)[12] : Inf
## ARIMA(1,1,1)(1,1,0)[12] : 1709.104
## ARIMA(1,1,1)(1,1,1)[12] : Inf
## ARIMA(1,1,1)(1,1,2)[12] : Inf
## ARIMA(1,1,1)(2,1,0)[12] : 1675.961
## ARIMA(1,1,1)(2,1,1)[12] : Inf
## ARIMA(1,1,2)(0,1,0)[12] : 1745.243
## ARIMA(1,1,2)(0,1,1)[12] : Inf
## ARIMA(1,1,2)(0,1,2)[12] : Inf
## ARIMA(1,1,2)(1,1,0)[12] : 1707.347
## ARIMA(1,1,2)(1,1,1)[12] : Inf
## ARIMA(1,1,2)(2,1,0)[12] : 1674.394
## ARIMA(1,1,3)(0,1,0)[12] : Inf
## ARIMA(1,1,3)(0,1,1)[12] : Inf
## ARIMA(1,1,3)(1,1,0)[12] : 1709.318
## ARIMA(1,1,4)(0,1,0)[12] : Inf
## ARIMA(2,1,0)(0,1,0)[12] : 1748.25
## ARIMA(2,1,0)(0,1,1)[12] : Inf
## ARIMA(2,1,0)(0,1,2)[12] : Inf
## ARIMA(2,1,0)(1,1,0)[12] : 1709.099
## ARIMA(2,1,0)(1,1,1)[12] : Inf
## ARIMA(2,1,0)(1,1,2)[12] : Inf
## ARIMA(2,1,0)(2,1,0)[12] : Inf
## ARIMA(2,1,0)(2,1,1)[12] : Inf
## ARIMA(2,1,1)(0,1,0)[12] : 1748.947
## ARIMA(2,1,1)(0,1,1)[12] : Inf
## ARIMA(2,1,1)(0,1,2)[12] : Inf
## ARIMA(2,1,1)(1,1,0)[12] : 1710.516
## ARIMA(2,1,1)(1,1,1)[12] : Inf
## ARIMA(2,1,1)(2,1,0)[12] : 1677.431
## ARIMA(2,1,2)(0,1,0)[12] : 1747.309
## ARIMA(2,1,2)(0,1,1)[12] : Inf
## ARIMA(2,1,2)(1,1,0)[12] : 1709.262
## ARIMA(2,1,3)(0,1,0)[12] : Inf
## ARIMA(3,1,0)(0,1,0)[12] : 1745.909
## ARIMA(3,1,0)(0,1,1)[12] : Inf
## ARIMA(3,1,0)(0,1,2)[12] : Inf
## ARIMA(3,1,0)(1,1,0)[12] : 1708.073
## ARIMA(3,1,0)(1,1,1)[12] : Inf
## ARIMA(3,1,0)(2,1,0)[12] : 1674.943
## ARIMA(3,1,1)(0,1,0)[12] : 1747.47
## ARIMA(3,1,1)(0,1,1)[12] : Inf
## ARIMA(3,1,1)(1,1,0)[12] : 1709.859
## ARIMA(3,1,2)(0,1,0)[12] : Inf
## ARIMA(4,1,0)(0,1,0)[12] : 1747.399
## ARIMA(4,1,0)(0,1,1)[12] : Inf
## ARIMA(4,1,0)(1,1,0)[12] : 1709.745
## ARIMA(4,1,1)(0,1,0)[12] : 1749.49
## ARIMA(5,1,0)(0,1,0)[12] : 1749.485
##
##
##
## Best model: ARIMA(0,1,2)(2,1,0)[12]print(summary(brent_arima))## Series: brent_ts
## ARIMA(0,1,2)(2,1,0)[12]
##
## Coefficients:
## ma1 ma2 sar1 sar2
## 0.3968 0.2070 -0.5493 -0.3727
## s.e. 0.0627 0.0638 0.0607 0.0597
##
## sigma^2 = 37.47: log likelihood = -831.23
## AIC=1672.45 AICc=1672.69 BIC=1690.2
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.1428664 5.925587 4.309494 0.2114035 7.203224 0.2410575
## ACF1
## Training set -0.006963124checkresiduals(brent_arima)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,2)(2,1,0)[12]
## Q* = 24.103, df = 20, p-value = 0.2379
##
## Model df: 4. Total lags used: 24ETS model has much lower SD, but more unnatural shape. We will plot ETS and Arima model
forecast_brent_ets <-forecast(brent_ets,h=24)
autoplot(forecast_brent_ets) + ggtitle('US Brent Oil Price forecast with ETS method') + theme(plot.title = element_text(hjust = 0.5)) +
labs(x="Year",y="Price (USD)")
forecast_arima_ets <-forecast(brent_arima,h=24)
autoplot(forecast_arima_ets) + ggtitle('US Brent Oil Price forecast with Arima method') + theme(plot.title = element_text(hjust = 0.5)) +
labs(x="Year",y="Price (USD)")
We do the same for Crude Oil
crude_s <- snaive(crude_diff)
print(summary(crude_s))##
## Forecast method: Seasonal naive method
##
## Model Information:
## Call: snaive(y = crude_diff)
##
## Residual sd: 7.8191
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.1967223 7.819116 5.591141 -382.3438 896.0021 1 0.4181407
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jul 2022 0.8595356 -9.1610647 10.880136 -14.465650 16.18472
## Aug 2022 -2.5674901 -12.5880904 7.453110 -17.892676 12.75770
## Sep 2022 3.8924901 -6.1281102 13.913090 -11.432696 19.21768
## Oct 2022 7.2778778 -2.7427225 17.298478 -8.047308 22.60306
## Nov 2022 -1.5978680 -11.6184683 8.422732 -16.923054 13.72732
## Dec 2022 -3.7102273 -13.7308276 6.310373 -19.035413 11.61496
## Jan 2023 9.8343182 -0.1862821 19.854918 -5.490868 25.15950
## Feb 2023 6.9040152 -3.1165852 16.924615 -8.421171 22.22920
## Mar 2023 14.0388406 4.0182403 24.059441 -1.286345 29.36403
## Apr 2023 -0.4421739 -10.4627742 9.578426 -15.767360 14.88301
## May 2023 8.4834783 -1.5371221 18.504079 -6.841707 23.80866
## Jun 2023 9.8207525 -0.1998478 19.841353 -5.504433 25.14594
## Jul 2023 0.8595356 -13.3117333 15.030804 -20.813550 22.53262
## Aug 2023 -2.5674901 -16.7387590 11.603779 -24.240576 19.10560
## Sep 2023 3.8924901 -10.2787787 18.063759 -17.780595 25.56558
## Oct 2023 7.2778778 -6.8933910 21.449147 -14.395208 28.95096
## Nov 2023 -1.5978680 -15.7691368 12.573401 -23.270953 20.07522
## Dec 2023 -3.7102273 -17.8814961 10.461042 -25.383313 17.96286
## Jan 2024 9.8343182 -4.3369507 24.005587 -11.838767 31.50740
## Feb 2024 6.9040152 -7.2672537 21.075284 -14.769070 28.57710
## Mar 2024 14.0388406 -0.1324283 28.210109 -7.634245 35.71193
## Apr 2024 -0.4421739 -14.6134428 13.729095 -22.115259 21.23091
## May 2024 8.4834783 -5.6877906 22.654747 -13.189607 30.15656
## Jun 2024 9.8207525 -4.3505164 23.992021 -11.852333 31.49384checkresiduals(crude_s)
##
## Ljung-Box test
##
## data: Residuals from Seasonal naive method
## Q* = 156.24, df = 24, p-value < 2.2e-16
##
## Model df: 0. Total lags used: 24crude_ets <- ets(crude_df_ts)
print(summary(crude_ets))## ETS(M,A,N)
##
## Call:
## ets(y = crude_df_ts)
##
## Smoothing parameters:
## alpha = 0.9999
## beta = 1e-04
##
## Initial states:
## l = 27.1766
## b = 0.7252
##
## sigma: 0.0926
##
## AIC AICc BIC
## 2412.398 2412.625 2430.390
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.3819662 5.540071 3.894827 -1.364705 6.950221 0.2290028
## ACF1
## Training set 0.3895083checkresiduals(crude_ets)
##
## Ljung-Box test
##
## data: Residuals from ETS(M,A,N)
## Q* = 39.408, df = 20, p-value = 0.00593
##
## Model df: 4. Total lags used: 24crude_arima <- auto.arima(crude_df_ts,d=1,D=1,stepwise=FALSE,approximation=FALSE,trace=TRUE)##
## ARIMA(0,1,0)(0,1,0)[12] : 1788.428
## ARIMA(0,1,0)(0,1,1)[12] : Inf
## ARIMA(0,1,0)(0,1,2)[12] : Inf
## ARIMA(0,1,0)(1,1,0)[12] : 1736.605
## ARIMA(0,1,0)(1,1,1)[12] : Inf
## ARIMA(0,1,0)(1,1,2)[12] : Inf
## ARIMA(0,1,0)(2,1,0)[12] : 1710.497
## ARIMA(0,1,0)(2,1,1)[12] : Inf
## ARIMA(0,1,0)(2,1,2)[12] : Inf
## ARIMA(0,1,1)(0,1,0)[12] : 1751.481
## ARIMA(0,1,1)(0,1,1)[12] : Inf
## ARIMA(0,1,1)(0,1,2)[12] : Inf
## ARIMA(0,1,1)(1,1,0)[12] : 1706.136
## ARIMA(0,1,1)(1,1,1)[12] : Inf
## ARIMA(0,1,1)(1,1,2)[12] : Inf
## ARIMA(0,1,1)(2,1,0)[12] : 1676.729
## ARIMA(0,1,1)(2,1,1)[12] : Inf
## ARIMA(0,1,1)(2,1,2)[12] : Inf
## ARIMA(0,1,2)(0,1,0)[12] : 1740.153
## ARIMA(0,1,2)(0,1,1)[12] : Inf
## ARIMA(0,1,2)(0,1,2)[12] : Inf
## ARIMA(0,1,2)(1,1,0)[12] : 1697.242
## ARIMA(0,1,2)(1,1,1)[12] : Inf
## ARIMA(0,1,2)(1,1,2)[12] : Inf
## ARIMA(0,1,2)(2,1,0)[12] : 1668.219
## ARIMA(0,1,2)(2,1,1)[12] : Inf
## ARIMA(0,1,3)(0,1,0)[12] : 1741.967
## ARIMA(0,1,3)(0,1,1)[12] : Inf
## ARIMA(0,1,3)(0,1,2)[12] : Inf
## ARIMA(0,1,3)(1,1,0)[12] : 1699.318
## ARIMA(0,1,3)(1,1,1)[12] : Inf
## ARIMA(0,1,3)(2,1,0)[12] : 1670.27
## ARIMA(0,1,4)(0,1,0)[12] : 1743.966
## ARIMA(0,1,4)(0,1,1)[12] : Inf
## ARIMA(0,1,4)(1,1,0)[12] : 1701.344
## ARIMA(0,1,5)(0,1,0)[12] : 1745.871
## ARIMA(1,1,0)(0,1,0)[12] : 1740.967
## ARIMA(1,1,0)(0,1,1)[12] : Inf
## ARIMA(1,1,0)(0,1,2)[12] : Inf
## ARIMA(1,1,0)(1,1,0)[12] : 1698.235
## ARIMA(1,1,0)(1,1,1)[12] : Inf
## ARIMA(1,1,0)(1,1,2)[12] : Inf
## ARIMA(1,1,0)(2,1,0)[12] : 1668.314
## ARIMA(1,1,0)(2,1,1)[12] : Inf
## ARIMA(1,1,0)(2,1,2)[12] : Inf
## ARIMA(1,1,1)(0,1,0)[12] : 1742.91
## ARIMA(1,1,1)(0,1,1)[12] : Inf
## ARIMA(1,1,1)(0,1,2)[12] : Inf
## ARIMA(1,1,1)(1,1,0)[12] : 1700.193
## ARIMA(1,1,1)(1,1,1)[12] : Inf
## ARIMA(1,1,1)(1,1,2)[12] : Inf
## ARIMA(1,1,1)(2,1,0)[12] : 1670.352
## ARIMA(1,1,1)(2,1,1)[12] : Inf
## ARIMA(1,1,2)(0,1,0)[12] : 1742.01
## ARIMA(1,1,2)(0,1,1)[12] : Inf
## ARIMA(1,1,2)(0,1,2)[12] : Inf
## ARIMA(1,1,2)(1,1,0)[12] : 1699.318
## ARIMA(1,1,2)(1,1,1)[12] : Inf
## ARIMA(1,1,2)(2,1,0)[12] : Inf
## ARIMA(1,1,3)(0,1,0)[12] : 1744.009
## ARIMA(1,1,3)(0,1,1)[12] : Inf
## ARIMA(1,1,3)(1,1,0)[12] : Inf
## ARIMA(1,1,4)(0,1,0)[12] : Inf
## ARIMA(2,1,0)(0,1,0)[12] : 1742.855
## ARIMA(2,1,0)(0,1,1)[12] : Inf
## ARIMA(2,1,0)(0,1,2)[12] : Inf
## ARIMA(2,1,0)(1,1,0)[12] : 1700.134
## ARIMA(2,1,0)(1,1,1)[12] : Inf
## ARIMA(2,1,0)(1,1,2)[12] : Inf
## ARIMA(2,1,0)(2,1,0)[12] : Inf
## ARIMA(2,1,0)(2,1,1)[12] : Inf
## ARIMA(2,1,1)(0,1,0)[12] : 1743.826
## ARIMA(2,1,1)(0,1,1)[12] : Inf
## ARIMA(2,1,1)(0,1,2)[12] : Inf
## ARIMA(2,1,1)(1,1,0)[12] : 1701.125
## ARIMA(2,1,1)(1,1,1)[12] : Inf
## ARIMA(2,1,1)(2,1,0)[12] : 1671.842
## ARIMA(2,1,2)(0,1,0)[12] : 1743.853
## ARIMA(2,1,2)(0,1,1)[12] : Inf
## ARIMA(2,1,2)(1,1,0)[12] : 1701.364
## ARIMA(2,1,3)(0,1,0)[12] : Inf
## ARIMA(3,1,0)(0,1,0)[12] : 1741.768
## ARIMA(3,1,0)(0,1,1)[12] : Inf
## ARIMA(3,1,0)(0,1,2)[12] : Inf
## ARIMA(3,1,0)(1,1,0)[12] : 1699.259
## ARIMA(3,1,0)(1,1,1)[12] : Inf
## ARIMA(3,1,0)(2,1,0)[12] : 1670.029
## ARIMA(3,1,1)(0,1,0)[12] : 1743.839
## ARIMA(3,1,1)(0,1,1)[12] : Inf
## ARIMA(3,1,1)(1,1,0)[12] : 1701.354
## ARIMA(3,1,2)(0,1,0)[12] : 1745.924
## ARIMA(4,1,0)(0,1,0)[12] : 1743.839
## ARIMA(4,1,0)(0,1,1)[12] : Inf
## ARIMA(4,1,0)(1,1,0)[12] : 1701.354
## ARIMA(4,1,1)(0,1,0)[12] : Inf
## ARIMA(5,1,0)(0,1,0)[12] : 1745.93
##
##
##
## Best model: ARIMA(0,1,2)(2,1,0)[12]print(summary(crude_arima))## Series: crude_df_ts
## ARIMA(0,1,2)(2,1,0)[12]
##
## Coefficients:
## ma1 ma2 sar1 sar2
## 0.3949 0.2102 -0.5527 -0.3453
## s.e. 0.0612 0.0621 0.0596 0.0591
##
## sigma^2 = 36.89: log likelihood = -828.99
## AIC=1667.98 AICc=1668.22 BIC=1685.73
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.1509411 5.879033 4.391245 0.1460376 7.703617 0.2581906
## ACF1
## Training set 0.003493412checkresiduals(crude_arima)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,2)(2,1,0)[12]
## Q* = 22.846, df = 20, p-value = 0.2964
##
## Model df: 4. Total lags used: 24forecast_crude_ets <-forecast(crude_ets,h=24)
autoplot(forecast_crude_ets) + ggtitle('US Crude Oil WTI Price forecast with ETS method') + theme(plot.title = element_text(hjust = 0.5)) +
labs(x="Year",y="Price (USD)")
forecast_crude_arrima <-forecast(crude_arima,h=24)
autoplot(forecast_crude_arrima) + ggtitle('US Crude Oil WTI Price forecast with Arima method') + theme(plot.title = element_text(hjust = 0.5)) +
labs(x="Year",y="Price (USD)")
Heating Oil
heating_oil_s <- snaive(heating_oil_diff)
print(summary(heating_oil_s))##
## Forecast method: Seasonal naive method
##
## Model Information:
## Call: snaive(y = heating_oil_diff)
##
## Residual sd: 0.2136
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.008388734 0.2135632 0.1566725 479.5766 784.813 1 0.3495745
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jul 2022 0.01152579 -0.26216644 0.2852180 -0.40705035 0.4301019
## Aug 2022 -0.06048715 -0.33417938 0.2132051 -0.47906330 0.3580890
## Sep 2022 0.12281976 -0.15087247 0.3965120 -0.29575638 0.5413959
## Oct 2022 0.31735859 0.04366636 0.5910508 -0.10121755 0.7359347
## Nov 2022 -0.12505790 -0.39875013 0.1486343 -0.54363405 0.2935182
## Dec 2022 -0.13837500 -0.41206723 0.1353172 -0.55695115 0.2802011
## Jan 2023 0.36082500 0.08713277 0.6345172 -0.05775115 0.7794011
## Feb 2023 0.23050552 -0.04318671 0.5041977 -0.18807063 0.6490817
## Mar 2023 0.81306770 0.53937547 1.0867599 0.39449156 1.2316438
## Apr 2023 0.16728944 -0.10640279 0.4409817 -0.25128670 0.5858656
## May 2023 0.11145839 -0.16223385 0.3851506 -0.30711776 0.5300345
## Jun 2023 0.40616304 0.13247081 0.6798553 -0.01241310 0.8247392
## Jul 2023 0.01152579 -0.37553347 0.3985851 -0.58043027 0.6034819
## Aug 2023 -0.06048715 -0.44754642 0.3265721 -0.65244322 0.5314689
## Sep 2023 0.12281976 -0.26423950 0.5098790 -0.46913630 0.7147758
## Oct 2023 0.31735859 -0.06970067 0.7044179 -0.27459747 0.9093147
## Nov 2023 -0.12505790 -0.51211716 0.2620014 -0.71701396 0.4668982
## Dec 2023 -0.13837500 -0.52543426 0.2486843 -0.73033106 0.4535811
## Jan 2024 0.36082500 -0.02623426 0.7478843 -0.23113106 0.9527811
## Feb 2024 0.23050552 -0.15655374 0.6175648 -0.36145054 0.8224616
## Mar 2024 0.81306770 0.42600844 1.2001270 0.22111164 1.4050238
## Apr 2024 0.16728944 -0.21976982 0.5543487 -0.42466662 0.7592455
## May 2024 0.11145839 -0.27560088 0.4985176 -0.48049768 0.7034144
## Jun 2024 0.40616304 0.01910378 0.7932223 -0.18579302 0.9981191checkresiduals(heating_oil_s)
##
## Ljung-Box test
##
## data: Residuals from Seasonal naive method
## Q* = 123.81, df = 24, p-value = 1.998e-15
##
## Model df: 0. Total lags used: 24heating_oil_ets <- ets(heating_oil_ts)
print(summary(heating_oil_ets))## ETS(M,N,N)
##
## Call:
## ets(y = heating_oil_ts)
##
## Smoothing parameters:
## alpha = 0.9999
##
## Initial states:
## l = 0.7748
##
## sigma: 0.0829
##
## AIC AICc BIC
## 446.5427 446.6329 457.3380
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.01318218 0.1558693 0.1095806 0.2840788 6.277259 0.2138429
## ACF1
## Training set 0.3471525checkresiduals(heating_oil_ets)
##
## Ljung-Box test
##
## data: Residuals from ETS(M,N,N)
## Q* = 54.907, df = 22, p-value = 0.0001224
##
## Model df: 2. Total lags used: 24heating_oil_arima <- auto.arima(heating_oil_ts,d=1,D=1,stepwise=FALSE,approximation=FALSE,trace=TRUE)##
## ARIMA(0,1,0)(0,1,0)[12] : -62.1747
## ARIMA(0,1,0)(0,1,1)[12] : Inf
## ARIMA(0,1,0)(0,1,2)[12] : Inf
## ARIMA(0,1,0)(1,1,0)[12] : -101.0199
## ARIMA(0,1,0)(1,1,1)[12] : Inf
## ARIMA(0,1,0)(1,1,2)[12] : Inf
## ARIMA(0,1,0)(2,1,0)[12] : -135.9307
## ARIMA(0,1,0)(2,1,1)[12] : Inf
## ARIMA(0,1,0)(2,1,2)[12] : Inf
## ARIMA(0,1,1)(0,1,0)[12] : -84.33342
## ARIMA(0,1,1)(0,1,1)[12] : Inf
## ARIMA(0,1,1)(0,1,2)[12] : Inf
## ARIMA(0,1,1)(1,1,0)[12] : -122.9344
## ARIMA(0,1,1)(1,1,1)[12] : Inf
## ARIMA(0,1,1)(1,1,2)[12] : Inf
## ARIMA(0,1,1)(2,1,0)[12] : -157.7046
## ARIMA(0,1,1)(2,1,1)[12] : Inf
## ARIMA(0,1,1)(2,1,2)[12] : Inf
## ARIMA(0,1,2)(0,1,0)[12] : -98.31443
## ARIMA(0,1,2)(0,1,1)[12] : Inf
## ARIMA(0,1,2)(0,1,2)[12] : Inf
## ARIMA(0,1,2)(1,1,0)[12] : -133.6174
## ARIMA(0,1,2)(1,1,1)[12] : Inf
## ARIMA(0,1,2)(1,1,2)[12] : Inf
## ARIMA(0,1,2)(2,1,0)[12] : -168.4654
## ARIMA(0,1,2)(2,1,1)[12] : Inf
## ARIMA(0,1,3)(0,1,0)[12] : -96.25611
## ARIMA(0,1,3)(0,1,1)[12] : Inf
## ARIMA(0,1,3)(0,1,2)[12] : Inf
## ARIMA(0,1,3)(1,1,0)[12] : -131.5377
## ARIMA(0,1,3)(1,1,1)[12] : Inf
## ARIMA(0,1,3)(2,1,0)[12] : -166.5195
## ARIMA(0,1,4)(0,1,0)[12] : -95.22559
## ARIMA(0,1,4)(0,1,1)[12] : Inf
## ARIMA(0,1,4)(1,1,0)[12] : -130.2402
## ARIMA(0,1,5)(0,1,0)[12] : -93.48771
## ARIMA(1,1,0)(0,1,0)[12] : -93.97742
## ARIMA(1,1,0)(0,1,1)[12] : Inf
## ARIMA(1,1,0)(0,1,2)[12] : Inf
## ARIMA(1,1,0)(1,1,0)[12] : -131.2266
## ARIMA(1,1,0)(1,1,1)[12] : Inf
## ARIMA(1,1,0)(1,1,2)[12] : Inf
## ARIMA(1,1,0)(2,1,0)[12] : -166.1977
## ARIMA(1,1,0)(2,1,1)[12] : Inf
## ARIMA(1,1,0)(2,1,2)[12] : Inf
## ARIMA(1,1,1)(0,1,0)[12] : -94.52147
## ARIMA(1,1,1)(0,1,1)[12] : Inf
## ARIMA(1,1,1)(0,1,2)[12] : Inf
## ARIMA(1,1,1)(1,1,0)[12] : -131.0916
## ARIMA(1,1,1)(1,1,1)[12] : Inf
## ARIMA(1,1,1)(1,1,2)[12] : Inf
## ARIMA(1,1,1)(2,1,0)[12] : -166.3786
## ARIMA(1,1,1)(2,1,1)[12] : Inf
## ARIMA(1,1,2)(0,1,0)[12] : -96.2625
## ARIMA(1,1,2)(0,1,1)[12] : Inf
## ARIMA(1,1,2)(0,1,2)[12] : Inf
## ARIMA(1,1,2)(1,1,0)[12] : -131.5384
## ARIMA(1,1,2)(1,1,1)[12] : Inf
## ARIMA(1,1,2)(2,1,0)[12] : -166.584
## ARIMA(1,1,3)(0,1,0)[12] : -94.36505
## ARIMA(1,1,3)(0,1,1)[12] : Inf
## ARIMA(1,1,3)(1,1,0)[12] : -129.61
## ARIMA(1,1,4)(0,1,0)[12] : -94.04817
## ARIMA(2,1,0)(0,1,0)[12] : -96.01226
## ARIMA(2,1,0)(0,1,1)[12] : Inf
## ARIMA(2,1,0)(0,1,2)[12] : Inf
## ARIMA(2,1,0)(1,1,0)[12] : -131.9769
## ARIMA(2,1,0)(1,1,1)[12] : Inf
## ARIMA(2,1,0)(1,1,2)[12] : Inf
## ARIMA(2,1,0)(2,1,0)[12] : -167.4028
## ARIMA(2,1,0)(2,1,1)[12] : Inf
## ARIMA(2,1,1)(0,1,0)[12] : -97.28405
## ARIMA(2,1,1)(0,1,1)[12] : Inf
## ARIMA(2,1,1)(0,1,2)[12] : Inf
## ARIMA(2,1,1)(1,1,0)[12] : -132.3607
## ARIMA(2,1,1)(1,1,1)[12] : Inf
## ARIMA(2,1,1)(2,1,0)[12] : -166.8888
## ARIMA(2,1,2)(0,1,0)[12] : -95.67125
## ARIMA(2,1,2)(0,1,1)[12] : Inf
## ARIMA(2,1,2)(1,1,0)[12] : -130.4752
## ARIMA(2,1,3)(0,1,0)[12] : -93.65252
## ARIMA(3,1,0)(0,1,0)[12] : -96.74054
## ARIMA(3,1,0)(0,1,1)[12] : Inf
## ARIMA(3,1,0)(0,1,2)[12] : Inf
## ARIMA(3,1,0)(1,1,0)[12] : -131.679
## ARIMA(3,1,0)(1,1,1)[12] : Inf
## ARIMA(3,1,0)(2,1,0)[12] : -166.7906
## ARIMA(3,1,1)(0,1,0)[12] : -95.60593
## ARIMA(3,1,1)(0,1,1)[12] : Inf
## ARIMA(3,1,1)(1,1,0)[12] : -130.4437
## ARIMA(3,1,2)(0,1,0)[12] : Inf
## ARIMA(4,1,0)(0,1,0)[12] : -95.52051
## ARIMA(4,1,0)(0,1,1)[12] : Inf
## ARIMA(4,1,0)(1,1,0)[12] : -130.4666
## ARIMA(4,1,1)(0,1,0)[12] : -93.57767
## ARIMA(5,1,0)(0,1,0)[12] : -93.52304
##
##
##
## Best model: ARIMA(0,1,2)(2,1,0)[12]print(summary(heating_oil_arima))## Series: heating_oil_ts
## ARIMA(0,1,2)(2,1,0)[12]
##
## Coefficients:
## ma1 ma2 sar1 sar2
## 0.3131 0.2186 -0.5480 -0.3858
## s.e. 0.0622 0.0591 0.0616 0.0600
##
## sigma^2 = 0.02898: log likelihood = 89.35
## AIC=-168.7 AICc=-168.47 BIC=-150.96
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.006373093 0.1647865 0.1196768 0.2368422 6.917821 0.2335453
## ACF1
## Training set 0.005184373checkresiduals(heating_oil_arima)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,2)(2,1,0)[12]
## Q* = 22.729, df = 20, p-value = 0.3023
##
## Model df: 4. Total lags used: 24forecast_heating_oil_ets <-forecast(heating_oil_ets,h=24)
autoplot(forecast_heating_oil_ets) + ggtitle('Heating Oil Price forecast with ETS method') + theme(plot.title = element_text(hjust = 0.5)) +
labs(x="Year",y="Price (USD)")
forecast_heating_oil_arrima <-forecast(heating_oil_arima,h=24)
autoplot(forecast_heating_oil_arrima) + ggtitle('Heating Oil Price forecast with Arima method') + theme(plot.title = element_text(hjust = 0.5)) +
labs(x="Year",y="Price (USD)")
Natural Gas
natural_gas_s <- snaive(ng_diff)
print(summary(natural_gas_s))##
## Forecast method: Seasonal naive method
##
## Model Information:
## Call: snaive(y = ng_diff)
##
## Residual sd: 0.9913
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -0.001600486 0.991289 0.6723625 245.7073 580.7684 1 0.1368467
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jul 2022 0.5326907 -0.73769721 1.80307863 -1.4101999 2.4755814
## Aug 2022 0.2294229 -1.04096499 1.49981084 -1.7134677 2.1723136
## Sep 2022 1.0438597 -0.22652823 2.31424760 -0.8990310 2.9867503
## Oct 2022 0.5299607 -0.74042725 1.80034858 -1.4129300 2.4728513
## Nov 2022 -0.4625931 -1.73298099 0.80779484 -2.4054837 1.4802976
## Dec 2022 -1.2407727 -2.51116064 0.02961519 -3.1836634 0.7021179
## Jan 2023 0.2573182 -1.01306974 1.52770610 -1.6855725 2.2002088
## Feb 2023 0.3579762 -0.91241173 1.62836411 -1.5849145 2.3008668
## Mar 2023 0.4369369 -0.83345106 1.70732477 -1.5059538 2.3798275
## Apr 2023 1.7686120 0.49822404 3.03899987 -0.1742787 3.7115026
## May 2023 1.4665620 0.19617404 2.73694987 -0.4763287 3.4094526
## Jun 2023 0.2398746 -1.03051334 1.51026250 -1.7030161 2.1827652
## Jul 2023 0.5326907 -1.26390911 2.32929053 -2.2149716 3.2803530
## Aug 2023 0.2294229 -1.56717690 2.02602275 -2.5182394 2.9770852
## Sep 2023 1.0438597 -0.75274014 2.84045951 -1.7038026 3.7915220
## Oct 2023 0.5299607 -1.26663916 2.32656048 -2.2177017 3.2776230
## Nov 2023 -0.4625931 -2.25919290 1.33400675 -3.2102554 2.2850692
## Dec 2023 -1.2407727 -3.03737255 0.55582709 -3.9884350 1.5068896
## Jan 2024 0.2573182 -1.53928164 2.05391800 -2.4903441 3.0049805
## Feb 2024 0.3579762 -1.43862363 2.15457601 -2.3896861 3.1056385
## Mar 2024 0.4369369 -1.35966297 2.23353667 -2.3107255 3.1845992
## Apr 2024 1.7686120 -0.02798787 3.56521178 -0.9790504 4.5162743
## May 2024 1.4665620 -0.33003787 3.26316178 -1.2811004 4.2142243
## Jun 2024 0.2398746 -1.55672524 2.03647440 -2.5077877 2.9875369checkresiduals(natural_gas_s)
##
## Ljung-Box test
##
## data: Residuals from Seasonal naive method
## Q* = 98.269, df = 24, p-value = 5.918e-11
##
## Model df: 0. Total lags used: 24natural_gas_ets <- ets(natural_gas_ts)
print(summary(natural_gas_ets))## ETS(M,N,N)
##
## Call:
## ets(y = natural_gas_ts)
##
## Smoothing parameters:
## alpha = 0.9999
##
## Initial states:
## l = 2.3436
##
## sigma: 0.1257
##
## AIC AICc BIC
## 1152.136 1152.227 1162.932
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.02252958 0.7083146 0.4623718 -0.2557278 9.450893 0.307449
## ACF1
## Training set 0.1568284checkresiduals(natural_gas_ets)
##
## Ljung-Box test
##
## data: Residuals from ETS(M,N,N)
## Q* = 33.677, df = 22, p-value = 0.05296
##
## Model df: 2. Total lags used: 24natural_gas_arima <- auto.arima(natural_gas_ts,d=1,D=1,stepwise=FALSE,approximation=FALSE,trace=TRUE)##
## ARIMA(0,1,0)(0,1,0)[12] : 726.853
## ARIMA(0,1,0)(0,1,1)[12] : Inf
## ARIMA(0,1,0)(0,1,2)[12] : Inf
## ARIMA(0,1,0)(1,1,0)[12] : 659.6415
## ARIMA(0,1,0)(1,1,1)[12] : Inf
## ARIMA(0,1,0)(1,1,2)[12] : 586.8115
## ARIMA(0,1,0)(2,1,0)[12] : 629.4426
## ARIMA(0,1,0)(2,1,1)[12] : 587.2666
## ARIMA(0,1,0)(2,1,2)[12] : Inf
## ARIMA(0,1,1)(0,1,0)[12] : 723.4769
## ARIMA(0,1,1)(0,1,1)[12] : Inf
## ARIMA(0,1,1)(0,1,2)[12] : 581.4574
## ARIMA(0,1,1)(1,1,0)[12] : 654.5283
## ARIMA(0,1,1)(1,1,1)[12] : 581.7518
## ARIMA(0,1,1)(1,1,2)[12] : 582.0649
## ARIMA(0,1,1)(2,1,0)[12] : 622.5578
## ARIMA(0,1,1)(2,1,1)[12] : 582.5242
## ARIMA(0,1,1)(2,1,2)[12] : Inf
## ARIMA(0,1,2)(0,1,0)[12] : 724.938
## ARIMA(0,1,2)(0,1,1)[12] : Inf
## ARIMA(0,1,2)(0,1,2)[12] : 583.5375
## ARIMA(0,1,2)(1,1,0)[12] : 656.5919
## ARIMA(0,1,2)(1,1,1)[12] : 583.8321
## ARIMA(0,1,2)(1,1,2)[12] : 584.1592
## ARIMA(0,1,2)(2,1,0)[12] : 624.637
## ARIMA(0,1,2)(2,1,1)[12] : 584.6181
## ARIMA(0,1,3)(0,1,0)[12] : 722.7901
## ARIMA(0,1,3)(0,1,1)[12] : Inf
## ARIMA(0,1,3)(0,1,2)[12] : 585.633
## ARIMA(0,1,3)(1,1,0)[12] : 656.6917
## ARIMA(0,1,3)(1,1,1)[12] : 585.9171
## ARIMA(0,1,3)(2,1,0)[12] : 626.5874
## ARIMA(0,1,4)(0,1,0)[12] : 722.3971
## ARIMA(0,1,4)(0,1,1)[12] : Inf
## ARIMA(0,1,4)(1,1,0)[12] : 655.326
## ARIMA(0,1,5)(0,1,0)[12] : 721.8086
## ARIMA(1,1,0)(0,1,0)[12] : 723.9236
## ARIMA(1,1,0)(0,1,1)[12] : Inf
## ARIMA(1,1,0)(0,1,2)[12] : 581.6218
## ARIMA(1,1,0)(1,1,0)[12] : 654.5443
## ARIMA(1,1,0)(1,1,1)[12] : 581.9192
## ARIMA(1,1,0)(1,1,2)[12] : 582.3083
## ARIMA(1,1,0)(2,1,0)[12] : 622.8369
## ARIMA(1,1,0)(2,1,1)[12] : 582.7678
## ARIMA(1,1,0)(2,1,2)[12] : 583.988
## ARIMA(1,1,1)(0,1,0)[12] : 724.3137
## ARIMA(1,1,1)(0,1,1)[12] : Inf
## ARIMA(1,1,1)(0,1,2)[12] : 583.5375
## ARIMA(1,1,1)(1,1,0)[12] : 656.562
## ARIMA(1,1,1)(1,1,1)[12] : 583.8322
## ARIMA(1,1,1)(1,1,2)[12] : Inf
## ARIMA(1,1,1)(2,1,0)[12] : Inf
## ARIMA(1,1,1)(2,1,1)[12] : Inf
## ARIMA(1,1,2)(0,1,0)[12] : Inf
## ARIMA(1,1,2)(0,1,1)[12] : Inf
## ARIMA(1,1,2)(0,1,2)[12] : 585.6309
## ARIMA(1,1,2)(1,1,0)[12] : 658.3649
## ARIMA(1,1,2)(1,1,1)[12] : 585.9255
## ARIMA(1,1,2)(2,1,0)[12] : Inf
## ARIMA(1,1,3)(0,1,0)[12] : 723.9258
## ARIMA(1,1,3)(0,1,1)[12] : Inf
## ARIMA(1,1,3)(1,1,0)[12] : 657.9806
## ARIMA(1,1,4)(0,1,0)[12] : 715.1417
## ARIMA(2,1,0)(0,1,0)[12] : 725.7714
## ARIMA(2,1,0)(0,1,1)[12] : Inf
## ARIMA(2,1,0)(0,1,2)[12] : 583.5722
## ARIMA(2,1,0)(1,1,0)[12] : 656.608
## ARIMA(2,1,0)(1,1,1)[12] : 583.8753
## ARIMA(2,1,0)(1,1,2)[12] : 584.1879
## ARIMA(2,1,0)(2,1,0)[12] : 624.7534
## ARIMA(2,1,0)(2,1,1)[12] : Inf
## ARIMA(2,1,1)(0,1,0)[12] : 726.3763
## ARIMA(2,1,1)(0,1,1)[12] : Inf
## ARIMA(2,1,1)(0,1,2)[12] : 585.6342
## ARIMA(2,1,1)(1,1,0)[12] : 658.507
## ARIMA(2,1,1)(1,1,1)[12] : Inf
## ARIMA(2,1,1)(2,1,0)[12] : Inf
## ARIMA(2,1,2)(0,1,0)[12] : Inf
## ARIMA(2,1,2)(0,1,1)[12] : Inf
## ARIMA(2,1,2)(1,1,0)[12] : 659.9851
## ARIMA(2,1,3)(0,1,0)[12] : Inf
## ARIMA(3,1,0)(0,1,0)[12] : 721.0856
## ARIMA(3,1,0)(0,1,1)[12] : Inf
## ARIMA(3,1,0)(0,1,2)[12] : 585.4517
## ARIMA(3,1,0)(1,1,0)[12] : 656.0038
## ARIMA(3,1,0)(1,1,1)[12] : 585.6901
## ARIMA(3,1,0)(2,1,0)[12] : 626.1129
## ARIMA(3,1,1)(0,1,0)[12] : 723.1651
## ARIMA(3,1,1)(0,1,1)[12] : Inf
## ARIMA(3,1,1)(1,1,0)[12] : 658.1007
## ARIMA(3,1,2)(0,1,0)[12] : Inf
## ARIMA(4,1,0)(0,1,0)[12] : 723.1642
## ARIMA(4,1,0)(0,1,1)[12] : Inf
## ARIMA(4,1,0)(1,1,0)[12] : 658.1005
## ARIMA(4,1,1)(0,1,0)[12] : Inf
## ARIMA(5,1,0)(0,1,0)[12] : 722.9748
##
##
##
## Best model: ARIMA(0,1,1)(0,1,2)[12]print(summary(natural_gas_arima))## Series: natural_gas_ts
## ARIMA(0,1,1)(0,1,2)[12]
##
## Coefficients:
## ma1 sma1 sma2
## 0.1648 -0.9795 0.0971
## s.e. 0.0619 0.0726 0.0713
##
## sigma^2 = 0.5113: log likelihood = -286.65
## AIC=581.3 AICc=581.46 BIC=595.49
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.02017634 0.6935477 0.4605215 -0.8374983 10.05886 0.3062187
## ACF1
## Training set -0.002685253checkresiduals(natural_gas_arima)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,1)(0,1,2)[12]
## Q* = 29.428, df = 21, p-value = 0.1041
##
## Model df: 3. Total lags used: 24forecast_natural_gas_ets <-forecast(natural_gas_ets,h=24)
autoplot(forecast_natural_gas_ets) + ggtitle('Natural Gas Price forecast with ETS method') + theme(plot.title = element_text(hjust = 0.5)) +
labs(x="Year",y="Price (USD)")
forecast_natural_gas_arrima <-forecast(natural_gas_arima,h=24)
autoplot(forecast_natural_gas_arrima) + ggtitle('Natural Gas Price forecast with Arima method') + theme(plot.title = element_text(hjust = 0.5)) +
labs(x="Year",y="Price (USD)")