Intro

We will use functions from several packages. These packages have to be installed, using th install.packages(“packagename”) function. Below, we install 3 packages at once using the c() argument: readxl, for reading Excel-files; fpp2, a package that includes many forecasting techniques; and seasonal.

You only have to install a package once. We already did it on our computers, and for that reason the function is preceded by a “#”, making it a comment, or reminder.

If you want to use the functions of these packages, you do have to use the library() commands in your session to invoke these functions!

# install.packages(c("readxl","fpp2","seasonal"))

library(fpp2)
library(readxl)
library(seasonal)

Next, we read the Excel-file with the four timeseries.

The fourth series is the most interesting, since it contains both trend and season. We will select that one and skip the rest, in a data frame ts.

To use our forecasting functions, the data frame has to be formatted as a time series. We convert ts to a timeseries z with the *ts()** function. The series is quarterly, so the frequency is 4 (four quarters per year). It starts in the first quarter of 2023.

We use autoplot() to plot the series.

simplemethods <- read_excel("simplemethods.xls") # Read Excel-file
ts <- simplemethods[,c("Periods","TS")]          # We select period and one series
names(ts) <- c("Quarter","TS")                   # We rename the columns

z <- ts(ts$TS, frequency = 4, start = c(2023, 1)) 
z
##      Qtr1 Qtr2 Qtr3 Qtr4
## 2023  287  352  416  313
## 2024  359  440  520  391
## 2025  449  550  650  489
autoplot(z)

We use other visualizations to visually analyze the series.

A subseries plot shows both trend and season. The plot breaks down the series by period (here: quarter). It is easy to see the seasonal pattern, as the level of observations differs by quarter. Q3 is the quarter where the series peaks.

The trend is indicated by the systematic increase in each of the four quarters, over the 3 years of data.

ggsubseriesplot(z) + ylab("Units") + ggtitle("Demand for Item 4")

We can now decompose the time series:

(dec <- decompose(z, type="additive"))
## $x
##      Qtr1 Qtr2 Qtr3 Qtr4
## 2023  287  352  416  313
## 2024  359  440  520  391
## 2025  449  550  650  489
## 
## $seasonal
##           Qtr1      Qtr2      Qtr3      Qtr4
## 2023 -38.46875  26.90625  75.03125 -63.46875
## 2024 -38.46875  26.90625  75.03125 -63.46875
## 2025 -38.46875  26.90625  75.03125 -63.46875
## 
## $trend
##        Qtr1   Qtr2   Qtr3   Qtr4
## 2023     NA     NA 351.00 371.00
## 2024 395.00 417.75 438.75 463.75
## 2025 493.75 522.25     NA     NA
## 
## $random
##           Qtr1      Qtr2      Qtr3      Qtr4
## 2023        NA        NA -10.03125   5.46875
## 2024   2.46875  -4.65625   6.21875  -9.28125
## 2025  -6.28125   0.84375        NA        NA
## 
## $figure
## [1] -38.46875  26.90625  75.03125 -63.46875
## 
## $type
## [1] "additive"
## 
## attr(,"class")
## [1] "decomposed.ts"
autoplot(dec)

We can use decomposition for forecasting. But whether it pays off to use decomposition instead of simple(r) methods, remains to be seen! We can use the simple(r) methods, like (seasonally) naive models, as benchmarks.

We can use the mean of the timeseries (12 quarters, in 2023-2025), to forecast the four quarters of 2026.

meanf(z, 4)
##         Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 2026 Q1       434.6667 285.4854 583.8479 193.8434 675.4899
## 2026 Q2       434.6667 285.4854 583.8479 193.8434 675.4899
## 2026 Q3       434.6667 285.4854 583.8479 193.8434 675.4899
## 2026 Q4       434.6667 285.4854 583.8479 193.8434 675.4899

We can plot the forecasts from various (naive, or simple) methods.

# Plot some forecasts
autoplot(z) +
  autolayer(meanf(z, h=8),  series="Mean", PI=FALSE) +
  autolayer(naive(z, h=8),  series="Naive", PI=FALSE) +
  autolayer(snaive(z, h=8), series="Seasonal naive", PI=FALSE) +
  autolayer(rwf(z, h=8, drift=TRUE), series="Drift", PI=FALSE) +
  ggtitle("Forecasts for item 4 (TS)") +
  xlab("Year/Quarter") + ylab("#") +
  guides(colour=guide_legend(title="Forecast"))