SES and LES
Jiantong Wang
Simple Exponential Smoothing
WFJ Example
In this example, we need load in the dataset of WFJ sales.
WFJ_Full <- read.csv("WFJ_Sales_Full.csv")
WFJ_Full.ts <- ts(WFJ_Full$WFJ.Sales,start = 1,frequency = 1)The only parameter we need for SES is \(\aplha\). The parameter can be user specified, or generated by R. Here for example, we want to make SES with alpha = 0.2, 0.5, and 0.7.
alpha_set <- c(0.2, 0.5, 0.7)
alpha_set## [1] 0.2 0.5 0.7Then let’s create a new matrix to store SES forecasts.
WFJ.frc.SES <- matrix(NA, nrow = 62, ncol = 3)
colnames(WFJ.frc.SES) <- paste("alpha=", alpha_set)
head(WFJ.frc.SES)## alpha= 0.2 alpha= 0.5 alpha= 0.7
## [1,] NA NA NA
## [2,] NA NA NA
## [3,] NA NA NA
## [4,] NA NA NA
## [5,] NA NA NA
## [6,] NA NA NAUse the iteration to calculate the predicted values
Let’s recall the formula from our lecture:
\[ F_{t+1}=L_{t+1}=(1-a) L_t+a y_t=L_t+a e_t \]
So, for each time, to calculate the prediction, all we need are last time period’s prediction error and smoothing parameter.
WFJ.frc.SES_iter <- WFJ.frc.SES
for(i in 2:62){
if(i == 2) {
predicted <- WFJ_Full.ts[1]
WFJ.frc.SES_iter[2,] <- rep(predicted,3)
}else{
for(k in 1:3){
alpha <- alpha_set[k]
error <- WFJ_Full.ts[(i-1)] - WFJ.frc.SES_iter[(i-1),k]
predicted <- WFJ.frc.SES_iter[(i-1),k] + alpha * error
WFJ.frc.SES_iter[i,k] <- predicted
}
}
}
head(WFJ.frc.SES_iter)## alpha= 0.2 alpha= 0.5 alpha= 0.7
## [1,] NA NA NA
## [2,] 23056.13 23056.13 23056.13
## [3,] 23408.37 23936.72 24288.96
## [4,] 23586.60 24118.14 24296.37
## [5,] 23517.78 23680.30 23558.64
## [6,] 23386.68 23271.30 23071.20use ets() function
ets() function is used to fit exponential smoothing state space model. We will discuss this kind of models later this semester.
library(forecast)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoofor (i in 2:62){
k <- i
ts.ref <- WFJ_Full.ts[1:(k-1)]
for (j in 1:3){ # we have 3 alphas
fit.model <- ets(ts.ref,model="ANN",alpha=alpha_set[j])
# in model part, the first letter A stands for our error term is additive, which means the random error is additive the prediction.
# the second letter denotes trend. For here, we do not have trend in our model, so it is N: None
# the third letter denotes season. For here, we do not have seasonal pattern in our model, so it's N.
WFJ.frc.SES[i,j] <- forecast::forecast(fit.model,h=1)$mean
}
}
head(WFJ.frc.SES)## alpha= 0.2 alpha= 0.5 alpha= 0.7
## [1,] NA NA NA
## [2,] 23056.13 23056.13 23056.13
## [3,] 23408.37 23936.72 24288.96
## [4,] 23586.60 24118.14 24296.37
## [5,] 23517.78 23680.30 23558.64
## [6,] 23386.68 23271.30 23071.20The results are slightly different from what we get through the manual calculation. It is because of the algorithm the function used, it did some approximation to boost the calculation.
Make time series plot
WFJ_SES_0.2.ts <- ts(WFJ.frc.SES[,1],start = 2,end = 62,frequency = 1)
WFJ_SES_0.5.ts <- ts(WFJ.frc.SES[,2],start = 2,end = 62,frequency = 1)
WFJ_SES_0.7.ts <- ts(WFJ.frc.SES[,3],start = 2,end = 62,frequency = 1)
autoplot(WFJ_Full.ts,ylab="Weekly Sales", main = "WFJ Sales & SES Prediction",xlab = "Time/Week", series = "Actual WFJ Sales") +
autolayer(WFJ_SES_0.2.ts, col = 2, series = "SES(0.2)") +
autolayer(WFJ_SES_0.5.ts, col = 3, series = "SES(0.5)") +
autolayer(WFJ_SES_0.7.ts, col = 3, series = "SES(0.7)") ## Warning: Removed 1 row containing missing values (`geom_line()`).
## Removed 1 row containing missing values (`geom_line()`).
## Removed 1 row containing missing values (`geom_line()`).
Compare the performance of each alpha
We can calculate the RMSE of the 3 models.
error_SES_0.2 <- WFJ_Full.ts - WFJ_SES_0.2.ts
error_SES_0.5 <- WFJ_Full.ts - WFJ_SES_0.5.ts
error_SES_0.7 <- WFJ_Full.ts - WFJ_SES_0.7.ts
RMSE_SES_0.2 <- sqrt(mean(error_SES_0.2^2, na.rm = T))
RMSE_SES_0.5 <- sqrt(mean(error_SES_0.5^2, na.rm = T))
RMSE_SES_0.7 <- sqrt(mean(error_SES_0.7^2, na.rm = T))
c(RMSE_SES_0.2, RMSE_SES_0.5, RMSE_SES_0.7)## [1] 4871.759 4490.625 4424.289Among 0.2, 0.5, 0.7, \(\alpha = 0.7\) is the best smoothing parameter.
Select Optimal Paramter by ets()
ets(WFJ_Full.ts,model="ANN")$par## alpha l
## 7.310805e-01 2.347550e+04The optimal \(\alpha\) selected by ets() is 0.731.
Linear Exponential Smoothing
Read in data and basic visualization
When a time series has a long-term trend, we should consider using linear exponential smoothing model. Here we have the total revenues of Netflix. First of all, let’s draw a time series plot of it.
Netflix <- read.csv("Netflix.csv")
Netflix.ts <- ts(Netflix$Total.Revenues,start = c(2000,1), frequency = 4)
autoplot(Netflix.ts)
From the time series plot, obviously, there is an increasing trend. An LES model should be adopted.
Find optimal parameters
We can use ets() function to help us find optimal parameters.
# Find optimal alpha and beta
ets(Netflix.ts ,model="AAN")$par## alpha beta l b
## 0.9999000 0.4420412 -0.3568988 3.7055700alpha_opt <- ets(Netflix.ts ,model="AAN")$par[1]We can find the optimal parameters for the LES model are \(\alpha = 1, \beta = 0.44\).
We can also specify \(\alpha\) by ourselves. For example, we also set \(\alpha\) = 0.2, 0.5, or 0.8. To find the corresponding optimal \(\beta\), we can also use ets() function.
alpha <- c(0.2,0.5,0.8,alpha_opt)
k <- length(alpha)
# Find optimal betas.
beta <- vector("numeric",k)
for (i in 1:k){
pars <- ets(Netflix.ts,model="AAN",alpha=alpha[i])$par
beta[i] <- pars[which(names(pars)=="beta")]
}Next is to create a matrix to store the predictions we are going to make. We set the column names as their parameters.
frc <- matrix(NA,nrow = 64,ncol = 4)
colnames(frc) <- paste("alpha=", alpha_set,"beta = ", round(beta,2))
head(frc)## alpha= 0.2 beta = 0.2 alpha= 0.5 beta = 0.5 alpha= 0.7 beta = 0.55
## [1,] NA NA NA
## [2,] NA NA NA
## [3,] NA NA NA
## [4,] NA NA NA
## [5,] NA NA NA
## [6,] NA NA NA
## alpha= 0.2 beta = 0.44
## [1,] NA
## [2,] NA
## [3,] NA
## [4,] NA
## [5,] NA
## [6,] NAUse ets() to make predictions
We have the alpha and the corresponding optimal beta, so we can build a linear exponential smoothing model by using function ets(). Please note here, we need to set the argument of model in ets() as “AAN”, which will make the function to build a LES model. For an LES model, we have the basic equation as F = L + T + E. T is standing for trend, and E is standing for random error.
n <- length(Netflix.ts)
# Produce forecasts
for (i in 3:64){
fitsample <- Netflix.ts[1:(i -1)]
fitsample <- ts(fitsample,frequency=frequency(Netflix.ts),start=start(Netflix.ts))
for (j in 1:k){
# First we fit ets with the fixed alpha value and model="AAN" - we will let it optimise the initial trend
fit <- ets(fitsample,model="AAN",alpha=alpha[j],beta=beta[j])
# Produce the 1-step ahead forecast
frc[i,j] <- forecast(fit,h=1)$mean
}
}
head(frc)## alpha= 0.2 beta = 0.2 alpha= 0.5 beta = 0.5 alpha= 0.7 beta = 0.55
## [1,] NA NA NA
## [2,] NA NA NA
## [3,] 10.18068 8.750063 8.362294
## [4,] 12.01510 11.180453 11.935743
## [5,] 14.17947 14.552785 15.862556
## [6,] 16.75927 18.698005 20.110058
## alpha= 0.2 beta = 0.44
## [1,] NA
## [2,] NA
## [3,] 8.633334
## [4,] 12.351960
## [5,] 16.019719
## [6,] 20.143650Compare performance
LES_1 <- ts(frc[,1],start = c(2000,1),frequency = 4)
LES_2 <- ts(frc[,2],start = c(2000,1),frequency = 4)
LES_3 <- ts(frc[,3],start = c(2000,1),frequency = 4)
LES_4 <- ts(frc[,4],start = c(2000,1),frequency = 4)
autoplot(Netflix.ts,series = "Actual", main = "Netflix Revenue", ylab = "Total Revenue", xlab = "Time")+
autolayer(LES_1) +
autolayer(LES_2) +
autolayer(LES_3) +
autolayer(LES_4,series = "Optimal one")## Warning: Removed 2 rows containing missing values (`geom_line()`).
## Removed 2 rows containing missing values (`geom_line()`).
## Removed 2 rows containing missing values (`geom_line()`).
## Removed 2 rows containing missing values (`geom_line()`).
sqrt(mean((Netflix.ts - LES_1)^2, na.rm = T))## [1] 42.317sqrt(mean((Netflix.ts - LES_2)^2, na.rm = T))## [1] 24.92471sqrt(mean((Netflix.ts - LES_3)^2, na.rm = T))## [1] 20.68453sqrt(mean((Netflix.ts - LES_4)^2, na.rm = T))## [1] 19.33488