Forecasting Workshop
Karla Gonzabay and Javier Calderon
2026-06-29
Introduction
Forecasting is an essential tool for analyzing time series data and making informed decisions. In this practice, three smoothing methods were applied, including Simple Exponential Smoothing and the Holt–Winters methods (Additive and Multiplicative). These techniques allow us to identify and model different components of a time series, such as level, trend, and seasonality, providing a basis for generating accurate future predictions. The goal of this exercise is to compare these methods and understand how the characteristics of the data influence the choice of the most suitable forecasting technique.
1. Simple Exponential Smoothing
We start reading the data set:
library(readxl)
ds1 <- read_excel("~/Simple Exponential smoothing data.xlsx",
range = "A5:C22", col_names = c("Year", "Time", "Observation"))
ds1## # A tibble: 18 × 3
## Year Time Observation
## <dbl> <dbl> <dbl>
## 1 1996 1 445.
## 2 1997 2 453.
## 3 1998 3 454.
## 4 1999 4 422.
## 5 2000 5 456.
## 6 2001 6 440.
## 7 2002 7 425.
## 8 2003 8 486.
## 9 2004 9 500.
## 10 2005 10 521.
## 11 2006 11 509.
## 12 2007 12 489.
## 13 2008 13 510.
## 14 2009 14 457.
## 15 2010 15 474.
## 16 2011 16 526.
## 17 2012 17 550.
## 18 2013 18 542.Then, we obtain the time series with ts() and plot it:
ts1 <- ts(ds1$Observation, start = 1996)
plot(ts1,
main = "Yearly demand from 1996 to 2013",
ylab = "Demand (in units)",
xlab = "Year",
type = "o")
Here we can see the time series from 1996 to 2013 and analyzing it there is no trend or seasonality. This time series only has level so we should use the Simple Exponential Smoothing.
Model Formulation
The Simple Exponential Smoothing (SES) method is used to forecast time series data that show no trend or seasonality, only a stable level with random variations. The smoothing and forecast equation are the following:
\[l_{t} = \alpha x_{t} + (1 - \alpha) l_{t-1}\]
\[\hat{x}_{t+h} = l_{t}\] lt is the level of the series at time t xt is the observation alpha is the smoothing parameter for the level
Forecasting
Now we have the right model, so we fit the model and next five years forecast, using different smoothing parameters:
- Optimal smoothing parameters:
library(fpp2)
fit1o <- ses(ts1, initial = "optimal", h = 5)
fit1o$model## Simple exponential smoothing
##
## Call:
## ses(y = ts1, h = 5, initial = "optimal")
##
## Smoothing parameters:
## alpha = 0.8339
##
## Initial states:
## l = 446.5649
##
## sigma: 29.8295
##
## AIC AICc BIC
## 178.1445 179.8588 180.8156- α = 0.1:
fit11 <- ses(ts1, initial = "optimal", h = 5, alpha = 0.1)
fit11$model$par## l alpha
## 470.416 0.100- α = 0.2:
fit12 <- ses(ts1, initial = "optimal", h = 5, alpha = 0.2)
fit12$model$par## l alpha
## 456.6441 0.2000- α = 0.3:
fit13 <- ses(ts1, initial = "optimal", h = 5, alpha = 0.3)
fit13$model$par## l alpha
## 450.2134 0.3000We can observe the level update per year until 2013 for each parameter:
fit1o$model$states## Time Series:
## Start = 1995
## End = 2013
## Frequency = 1
## l
## [1,] 446.5649
## [2,] 445.5602
## [3,] 451.9308
## [4,] 453.9981
## [5,] 427.6329
## [6,] 451.3206
## [7,] 442.2059
## [8,] 428.0169
## [9,] 476.5421
## [10,] 496.4614
## [11,] 517.1568
## [12,] 510.3134
## [13,] 492.4492
## [14,] 506.9758
## [15,] 465.0692
## [16,] 472.3662
## [17,] 517.0479
## [18,] 544.3838
## [19,] 542.6795fit11$model$states## Time Series:
## Start = 1995
## End = 2013
## Frequency = 1
## l
## [1,] 470.4160
## [2,] 467.9104
## [3,] 466.4394
## [4,] 465.2365
## [5,] 460.9508
## [6,] 460.4597
## [7,] 458.4528
## [8,] 455.1265
## [9,] 458.2348
## [10,] 462.4544
## [11,] 468.3369
## [12,] 472.3982
## [13,] 474.0474
## [14,] 477.6297
## [15,] 475.5387
## [16,] 475.3668
## [17,] 480.4251
## [18,] 487.3656
## [19,] 492.8631fit12$model$states## Time Series:
## Start = 1995
## End = 2013
## Frequency = 1
## l
## [1,] 456.6441
## [2,] 454.3873
## [3,] 454.1498
## [4,] 454.2019
## [5,] 447.8375
## [6,] 449.4780
## [7,] 447.6604
## [8,] 443.1663
## [9,] 451.7751
## [10,] 461.5060
## [11,] 473.4608
## [12,] 480.5587
## [13,] 482.2249
## [14,] 487.7539
## [15,] 481.5472
## [16,] 480.0017
## [17,] 489.1914
## [18,] 501.3191
## [19,] 509.5233fit13$model$states## Time Series:
## Start = 1995
## End = 2013
## Frequency = 1
## l
## [1,] 450.2134
## [2,] 448.7574
## [3,] 450.0902
## [4,] 451.3861
## [5,] 442.6843
## [6,] 446.6910
## [7,] 444.8007
## [8,] 438.9175
## [9,] 453.1052
## [10,] 467.3027
## [11,] 483.4959
## [12,] 491.1321
## [13,] 490.4595
## [14,] 496.2826
## [15,] 484.4138
## [16,] 481.2357
## [17,] 494.6500
## [18,] 511.2040
## [19,] 520.5448The specific fitted model with an \(\alpha\) of 0.8339 is:
\[l_{t} = 0.8339 x_{t} + (1 - 0.8339) l_{t-1}\]
\[\hat{x}_{t+h} = l_{t-1}\]
We can also see the fitted model and forecast:
fit1o$fitted #fitted values## Time Series:
## Start = 1996
## End = 2013
## Frequency = 1
## [1] 446.5649 445.5602 451.9308 453.9981 427.6329 451.3206 442.2059 428.0169
## [9] 476.5421 496.4614 517.1568 510.3134 492.4492 506.9758 465.0692 472.3662
## [17] 517.0479 544.3838fit1o$mean #Point forecast to h=5## Time Series:
## Start = 2014
## End = 2018
## Frequency = 1
## [1] 542.6795 542.6795 542.6795 542.6795 542.6795autoplot(ts1) +
autolayer(fit1o, series = "Optimized", PI=FALSE) +
autolayer(fit11, series = "α = 0.1", PI=FALSE) +
autolayer(fit12, series = "α = 0.2", PI=FALSE) +
autolayer(fit13, series = "α = 0.3", PI=FALSE) +
ggtitle("Yearly demand from 1996 to 2013") + ylab("Demand (in units)") +
xlab("Year") + guides (colour=guide_legend(title = "Forecast"))
Each forecast has a constant value since the next year observations we hope they are the same last value of the estimated level.
Residuals
We plot the residuals to analyze the correlations:
checkresiduals(fit1o)
##
## Ljung-Box test
##
## data: Residuals from Simple exponential smoothing
## Q* = 3.8065, df = 4, p-value = 0.4328
##
## Model df: 0. Total lags used: 4By analyzing the residuals, we can conclude that there is no correlation between the residuals and their mean is close to zero.
The last thing is to calculate the forecast metrics:
errors1 <- as.data.frame(rbind(
"α=0.8339" = round(accuracy(fit1o)[, c("MAE", "RMSE", "MAPE")], 2),
"α=0.1" = round(accuracy(fit11)[, c("MAE", "RMSE", "MAPE")], 2),
"α=0.2" = round(accuracy(fit12)[, c("MAE", "RMSE", "MAPE")], 2),
"α=0.3" = round(accuracy(fit13)[, c("MAE", "RMSE", "MAPE")], 2)
))
errors1## MAE RMSE MAPE
## α=0.8339 22.26 28.12 4.61
## α=0.1 31.97 36.97 6.50
## α=0.2 27.42 33.48 5.52
## α=0.3 25.49 31.28 5.16With this results, the model fitted with the optimal parameter give us the lowest error values, making this model the best compared to the others.
2. Additive Model
We start reading the data set:
ds3 <- read_excel("~/Holt Winter Additive data.xlsx",
range = "A7:C50", col_names = c("Year", "Time", "Observation"))Then, we obtain the time series with ts() and plot it:
ts3 <- ts(data=ds3$Observation, start = 2005, frequency = 4)
plot(ts3,
main = "Quarterly demand from 2005 to 2015",
ylab = "Demand (in units)",
xlab = "Year (in quarters)",
xaxt = "n",
ylim = c(15,80),
type = "o", pch = 0)
axis(1, at = seq(2005, 2015, by = 1), labels = seq(2005, 2015, by = 1))
Here we can see the time series from 2005 to 2015 by quarters and analyzing the time series shows an upward trend and a seasonality repeating each year. The seasonality shows a drop in the demand during the first quarter of each year, followed by a peak towards the end of each year. As we have here a level, an upward trend and a seasonality repeating each year with a constant magnitude we applied the Additive Triple Exponential Smoothing or Holt Winters Additive method.
Model Formulation
The Additive Holt–Winters Method, also known as Triple Exponential Smoothing, is used to forecast time series data that show both trend and seasonal patterns. It decomposes the series into three components: level, trend, and seasonality. The additive form is appropriate when the seasonal variations remain approximately constant over time. The forecasting, level, trend and seasonality equations are:
\[\hat{X}_{t+h} = l_{t} + hb_t + s_{t+h-m(k+1)}\]
\[l_{t} = \alpha (x_{t}-s_{t+h-m}) + (1 - \alpha) (l_{t-1} + b_{t-1})\] \[b_{t} = \beta (l_t - l_{t-1}) + (1 - \beta) b_{t-1}\] \[s_{t} = \gamma (x_t-l_{t-1}-b_{t-1}) + (1 - \gamma) s_{t-m}\]
Forecasting
Now we have the right model, so we fit the model and next eight quarters forecast, using different smoothing and trend parameters:
- Optimal smoothing parameters:
fit3o <- hw(ts3, h = 8, seasonal = "additive")
fit3o$model## Holt-Winters' additive method
##
## Call:
## hw(y = ts3, h = 8, seasonal = "additive")
##
## Smoothing parameters:
## alpha = 0.3063
## beta = 1e-04
## gamma = 0.4263
##
## Initial states:
## l = 32.2597
## b = 0.7014
## s = 1.3106 -1.6935 -9.3132 9.6962
##
## sigma: 1.9494
##
## AIC AICc BIC
## 234.4171 239.7112 250.4748- α = 0.2, beta = 0.002, gamma = 0.3:
fit31 <- hw(ts3, h = 8, alpha = 0.2, beta = 0.002, gamma = 0.3, seasonal = "additive")
fit31$model$par## l b s0 s1 s2 alpha beta
## 32.4610287 0.6831573 1.2309645 -1.7813503 -9.3442306 0.2000000 0.0020000
## gamma
## 0.3000000- α = 0.5, beta = 0.004, gamma = 0.1:
fit32 <- hw(ts3, h = 8, alpha = 0.5, beta = 0.004, gamma = 0.1, seasonal = "additive")
fit32$model$par## l b s0 s1 s2 alpha
## 31.7540544 0.7249231 1.1585945 -1.5112339 -10.4325933 0.5000000
## beta gamma
## 0.0040000 0.1000000We can also see the fitted model and forecast for the next two years:
fit3o$fitted #fitted values## Qtr1 Qtr2 Qtr3 Qtr4
## 2005 42.65723 24.21081 32.66618 36.37206
## 2006 45.53781 27.51872 36.21481 40.33713
## 2007 49.17134 32.18209 39.31365 43.50887
## 2008 49.89746 32.85434 39.71474 43.48280
## 2009 53.65576 35.82730 43.37632 45.34947
## 2010 56.83932 37.31366 45.42533 47.91400
## 2011 60.42249 37.71212 47.59100 49.91681
## 2012 63.19432 40.58637 49.32571 53.47546
## 2013 65.75739 44.06492 54.19500 55.53369
## 2014 68.24659 43.48545 54.81621 58.70628
## 2015 69.05311 47.59377 59.24376 64.22408fit3o$mean #Point forecast for h=8## Qtr1 Qtr2 Qtr3 Qtr4
## 2016 76.09837 51.60333 63.96867 68.37170
## 2017 78.90404 54.40899 66.77434 71.17737autoplot(ts3) +
autolayer(fit3o, series = "Optimized", PI=FALSE) +
autolayer(fit31, series = "α = 0.2, beta = 0.002, gamma = 0.3", PI=FALSE) +
autolayer(fit32, series = "α = 0.5, beta = 0.004, gamma = 0.1", PI=FALSE) +
ggtitle("Forecast with Additive Holt Winters Method") + ylab("Demand (in units)") +
xlab("Year") + guides (colour=guide_legend(title = "Forecast"))
Each forecast continue with the time series and we expect to see the same behavior for the next eight quarters (two years).
Residuals
We plot the residuals to analyze the correlations:
checkresiduals(fit3o)
##
## Ljung-Box test
##
## data: Residuals from Holt-Winters' additive method
## Q* = 4.5557, df = 8, p-value = 0.8038
##
## Model df: 0. Total lags used: 8By analyzing the residuals, we can conclude that there is a random behavior and no correlation between the residuals also their mean is close to zero.
The last thing is to calculate the forecast metrics:
errors3 <- as.data.frame(rbind(
"α=0.3063, beta=0.0001, gamma = 0.4263" = round(accuracy(fit3o)[, c("MAE", "RMSE", "MAPE")], 2),
"α=0.2, beta=0.002, gamma = 0.3" = round(accuracy(fit31)[, c("MAE", "RMSE", "MAPE")], 2),
"α=0.5, beta=0.004, gamma = 0.1" = round(accuracy(fit32)[, c("MAE", "RMSE", "MAPE")], 2)
))
errors3## MAE RMSE MAPE
## α=0.3063, beta=0.0001, gamma = 0.4263 1.37 1.76 2.97
## α=0.2, beta=0.002, gamma = 0.3 1.43 1.80 3.09
## α=0.5, beta=0.004, gamma = 0.1 1.67 2.00 3.75Based on the results obtained, the model fitted using the optimal parameters consistently produces the lowest error values among all the models tested. This indicates that, compared to the other approaches, it provides the most accurate and reliable forecasts, making it the best-performing model for this data set.
3. Multiplicative Model
We start reading the data set:
ds4 <- read_excel("~/Holt Winter multiplicative data.xlsx",
range = "A7:C50", col_names = c("Quarter", "Time", "Observation"))Then, we obtain the time series with ts() and plot it:
ts4 <- ts(data=ds4$Observation, start = 2005, frequency = 4)
plot(ts4,
main = "Quarterly demand from 2005 to 2015",
ylab = "Demand (in units)",
xlab = "Year (in quarters)",
xaxt = "n",
ylim = c(15,80),
type = "o", pch = 0)
axis(1, at = seq(2005, 2015, by = 1), labels = seq(2005, 2015, by = 1))
Here we can see the time series from 2005 to 2015, presented quarterly. The series exhibits an upward trend along with a seasonality pattern that repeats each year. Specifically, the seasonality shows a drop in demand during the first quarter, followed by a peak towards the end of the year. Since the series includes a level, an increasing trend, and seasonal fluctuations whose amplitude changes proportionally with the level, we applied the Multiplicative Triple Exponential Smoothing, also known as the Holt–Winters Multiplicative method.
Model Formulation
The Holt–Winters Multiplicative Method is a forecasting technique for time series with both trend and seasonal patterns, where seasonal fluctuations change proportionally with the level. It decomposes the series into level, trend, and seasonality, updating them with smoothing parameters to produce accurate forecasts for data with varying seasonal intensity. Here we have the model:
\[\hat{X}_{t+h} = (l_{t} + hb_t) s_{t+h-m(k+1)}\]
\[l_{t} = \alpha (x_{t}/s_{t+h-m}) + (1 - \alpha) (l_{t-1} + b_{t-1})\] \[b_{t} = \beta (l_t - l_{t-1}) + (1 - \beta) b_{t-1}\] \[s_{t} = \gamma (x_t/(l_{t-1}-b_{t-1})) + (1 - \gamma) s_{t-m}\]
Forecasting
Now we have the right model, so we fit the model and next eight quarters forecast, using different smoothing and trend parameters:
- Optimal smoothing parameters:
fit4o <- hw(ts4, h = 8, seasonal = "multiplicative")
fit4o$model## Holt-Winters' multiplicative method
##
## Call:
## hw(y = ts4, h = 8, seasonal = "multiplicative")
##
## Smoothing parameters:
## alpha = 0.4406
## beta = 0.0134
## gamma = 0.0023
##
## Initial states:
## l = 32.4875
## b = 0.6974
## s = 1.0237 0.9618 0.7704 1.2442
##
## sigma: 0.0367
##
## AIC AICc BIC
## 221.1313 226.4254 237.1890- α = 0.5, beta = 0.03, gamma = 0.001:
fit41 <- hw(ts4, h = 8, alpha = 0.5, beta = 0.03, gamma = 0.001,
seasonal = "multiplicative")
fit41$model$par## l b s0 s1 s2 alpha beta
## 32.8627034 0.7538879 1.0235106 0.9648253 0.7679187 0.5000000 0.0300000
## gamma
## 0.0010000- α = 0.6, beta = 0.02, gamma = 0.004:
fit42 <- hw(ts4, h = 8, alpha = 0.6, beta = 0.02, gamma = 0.004,
seasonal = "multiplicative")
fit42$model$par## l b s0 s1 s2 alpha beta
## 32.8823477 0.7639318 1.0232190 0.9649544 0.7685409 0.6000000 0.0200000
## gamma
## 0.0040000We can also see the fitted model and forecast for the next two years:
fit4o$fitted #fitted values## Qtr1 Qtr2 Qtr3 Qtr4
## 2005 41.28689 26.36042 32.62011 35.43591
## 2006 44.91855 28.43989 36.70598 39.63640
## 2007 50.24528 31.40864 39.84476 42.20533
## 2008 50.30134 31.91456 40.44309 44.16408
## 2009 53.76967 34.34469 43.00433 46.12632
## 2010 56.35302 36.28628 45.13587 48.57174
## 2011 58.97855 37.27806 47.78102 51.14259
## 2012 62.76659 39.02701 49.48525 54.86171
## 2013 67.26005 42.02854 52.46503 56.18995
## 2014 69.83611 42.43956 53.92965 58.43127
## 2015 72.59030 45.62125 58.77277 64.38368fit4o$mean #Point forecast for h=8## Qtr1 Qtr2 Qtr3 Qtr4
## 2016 80.08894 50.15482 63.34322 68.17810
## 2017 83.80112 52.45291 66.21274 71.23205autoplot(ts4) +
autolayer(fit4o, series = "Optimized", PI=FALSE) +
autolayer(fit41, series = "α = 0.5, beta = 0.03, gamma = 0.001", PI=FALSE) +
autolayer(fit42, series = "α = 0.6, beta = 0.02, gamma = 0.004", PI=FALSE) +
ggtitle("Forecast with Multiplicative Holt Winters Method") + ylab("Demand (in units)") +
xlab("Year (in quarters)") + guides (colour=guide_legend(title = "Forecast"))
The forecasts project the time series forward, anticipating a continuation of the observed behavior for the next eight quarters (two years).
Residuals
We plot the residuals to analyze the correlations:
checkresiduals(fit4o)
##
## Ljung-Box test
##
## data: Residuals from Holt-Winters' multiplicative method
## Q* = 4.8579, df = 8, p-value = 0.7727
##
## Model df: 0. Total lags used: 8From the residual analysis, we can conclude that the residuals are random, uncorrelated, and their average value is approximately zero.
The last thing is to calculate the forecast metrics:
errors4 <- as.data.frame(rbind(
"α=0.4406, beta=0.00134, gamma = 0.0023" = round(accuracy(fit4o)[, c("MAE", "RMSE", "MAPE")], 2),
"α=0.5, beta=0.03, gamma = 0.001" = round(accuracy(fit41)[, c("MAE", "RMSE", "MAPE")], 2),
"α=0.6, beta=0.02, gamma = 0.004" = round(accuracy(fit42)[, c("MAE", "RMSE", "MAPE")], 2)
))
errors4## MAE RMSE MAPE
## α=0.4406, beta=0.00134, gamma = 0.0023 1.25 1.58 2.71
## α=0.5, beta=0.03, gamma = 0.001 1.28 1.59 2.78
## α=0.6, beta=0.02, gamma = 0.004 1.30 1.61 2.82Based on the results, the model with the best parameters gives the lowest errors compared to the other models. This shows that it predicts the data most accurately and can be considered the best model for this data set.
Conclusion
In this practice, three different forecasting methods were applied: Simple Exponential Smoothing, Holt–Winters Additive, and Holt–Winters Multiplicative. Each method allowed us to explore different components of a time series, such as level, trend, and seasonality.
The analysis showed that selecting the appropriate method and tuning its parameters is essential to obtain accurate forecasts. Overall, the exercise highlights the importance of understanding the characteristics of a time series in order to choose the most suitable forecasting technique and generate reliable predictions.