library(fpp2)
## Loading required package: ggplot2
## Loading required package: forecast
## Loading required package: fma
## Loading required package: expsmooth
  1. Consider the pigs series — the number of pigs slaughtered in Victoria each month.
  1. Use the ses() function in R to find the optimal values of α and ℓ, and generate forecasts for the next four months.
ses_m1<- ses(pigs, h = 4)

# finding optimal parameters
ses_m1$model$par
##        alpha            l 
## 2.971488e-01 7.726006e+04
summary(ses_m1)
## 
## Forecast method: Simple exponential smoothing
## 
## Model Information:
## Simple exponential smoothing 
## 
## Call:
##  ses(y = pigs, h = 4) 
## 
##   Smoothing parameters:
##     alpha = 0.2971 
## 
##   Initial states:
##     l = 77260.0561 
## 
##   sigma:  10308.58
## 
##      AIC     AICc      BIC 
## 4462.955 4463.086 4472.665 
## 
## Error measures:
##                    ME    RMSE      MAE       MPE     MAPE      MASE
## Training set 385.8721 10253.6 7961.383 -0.922652 9.274016 0.7966249
##                    ACF1
## Training set 0.01282239
## 
## Forecasts:
##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Sep 1995       98816.41 85605.43 112027.4 78611.97 119020.8
## Oct 1995       98816.41 85034.52 112598.3 77738.83 119894.0
## Nov 1995       98816.41 84486.34 113146.5 76900.46 120732.4
## Dec 1995       98816.41 83958.37 113674.4 76092.99 121539.8

We can see that ses() has calculated the following optimal parametrs: α = 0.2971 and l0 = 77260.1 As α is small, more weight is given to observations from the more distant past.

autoplot(ses_m1) + autolayer(fitted(ses_m1), series = "Fitted") +
  ylab("Number Of Pigs Slaughtered")

  1. Compute a 95% prediction interval for the first forecast using ^y±1.96s where s is the standard deviation of the residuals. Compare your interval with the interval produced by R.
# calculating 95% prediction interval for the first forecast produced by R
ses_m1$upper[1, "95%"]
##      95% 
## 119020.8
ses_m1$lower[1, "95%"]
##      95% 
## 78611.97
# calculating the standard deviation of the residuals
s <- sd(ses_m1$residuals)

# calculating 95% prediction interval for the first forecast using formula ^y±1.96s
ses_m1$mean[1] + 1.96*s
## [1] 118952.8
ses_m1$mean[1] - 1.96*s
## [1] 78679.97

A 95% prediction interval calculated using formula ^y±1.96s and produced by R are similar.

  1. Modify your function from the previous exercise to return the sum of squared errors rather than the forecast of the next observation. Then use the optim() function to find the optimal values of α and ℓ0. Do you get the same values as the ses() function?

Solution:

The unknown parameters and the initial values for any exponential smoothing method can be estimated by minimising the SSE. Hence, we find the values of the unknown parameters and the initial values that minimise the SSE.

#  function takes argument as a vector, because optim() requires vector as the first argument of the function.
sse_fn <- function( pars=c(alpha,l0), y){
  
 # assigning the initial sse value
  sse <- 0
   #  assigning first parameter to alpha
  alpha <- pars[1]
   #  assigning second parametr to l0 
  l0 <- pars[2]
  y_hat <- l0
  
  #  calculating sse 
  for(index in 1:length(y)){
    sse <- sse + (y[index] - y_hat)^2
    y_hat <- alpha*y[index] + (1 - alpha)*y_hat
  }
  
  return(sse)
}
# finding the optimal values of α and ℓ0 using the optim(). Setting initial value of alpha = 0, l0 = first observation value of pigs data set
opt_par<- optim(par = c(0, pigs[1]), y = pigs, fn = sse_fn)
opt_par
## $par
## [1]     0.297086 77272.075070
## 
## $value
## [1] 19765613447
## 
## $counts
## function gradient 
##      179       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL
# finding the optimal values of α and ℓ0 using the optim()
ses_m1$model$par
##        alpha            l 
## 2.971488e-01 7.726006e+04

Optimal parameters calculated using optim() based on sum of squared errors and produced by R are similar.