Consider the the number of pigs slaughtered in Victoria, available in the
aus_livestockdataset.
- Use the
ETS()function in R to estimate the equivalent model for simple exponential smoothing. Find the optimal values of \(\alpha\) and \(\ell_0\), and generate forecasts for the next four months.
fit <- aus_livestock |>
filter(Animal == "Pigs", State == "Victoria") |>
model(ses = ETS(Count ~ error("A") + trend("N") + season("N")))
report(fit)Series: Count
Model: ETS(A,N,N)
Smoothing parameters:
alpha = 0.3221247
Initial states:
l[0]
100646.6
sigma^2: 87480760
AIC AICc BIC
13737.10 13737.14 13750.07 Optimal values are \(\alpha = 0.3221247\) and \(\ell_0 = 100646.6\)
fc <- fit |> forecast(h = "4 months")
fc# A fable: 4 x 6 [1M]
# Key: Animal, State, .model [1]
Animal State .model Month
<fct> <fct> <chr> <mth>
1 Pigs Victoria ses 2019 Jan
2 Pigs Victoria ses 2019 Feb
3 Pigs Victoria ses 2019 Mar
4 Pigs Victoria ses 2019 Apr
# ℹ 2 more variables: Count <dist>, .mean <dbl>fc |>
autoplot(filter(aus_livestock, Month >= yearmonth("2010 Jan")))
- Compute a 95% prediction interval for the first forecast using \(\hat{y} \pm 1.96s\) where \(s\) is the standard deviation of the residuals. Compare your interval with the interval produced by R.
s <- augment(fit) |>
pull(.resid) |>
sd()
yhat <- fc |>
pull(.mean) |>
head(1)
yhat + c(-1, 1) * 1.96 * s[1] 76871.01 113502.10fc |>
head(1) |>
mutate(interval = hilo(Count, 95)) |>
pull(interval)<hilo[1]>
[1] [76854.79, 113518.3]95The intervals are close but not identical. This is because R estimates the variance of the residuals differently, taking account of the degrees of freedom properly (and also using a more accurate critical value rather than just 1.96).
Try the following.
res <- augment(fit) |> pull(.resid)
s <- sqrt(sum(res^2) / (length(res) - NROW(tidy(fit))))
yhat + c(-1, 1) * qnorm(0.975) * s[1] 76854.79 113518.33Write your own function to implement simple exponential smoothing. The function should take arguments
y(the response data),alpha(the smoothing parameter \(\alpha\)) andlevel(the initial level \(\ell_0\)). It should return the forecast of the next observation in the series. Does it give the same forecast asETS()?
my_ses <- function(y, alpha, level) {
yhat <- numeric(length(y) + 1)
yhat[1] <- level
for (i in 2:(length(yhat))) {
yhat[i] <- alpha * y[i - 1] + (1 - alpha) * yhat[i - 1]
}
return(last(yhat))
}
vic_pigs_vec <- aus_livestock |>
filter(Animal == "Pigs", State == "Victoria") |>
pull(Count)
ses_fc <- vic_pigs_vec |>
my_ses(alpha = 0.3221, level = 100646.6)
c(my_ses = ses_fc, fable = fc$.mean[1]) my_ses fable
95186.59 95186.56 Yes, the same forecasts are obtained. The slight differences are due to rounding of \(\alpha\) and \(\ell_0\).
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 \(\alpha\) and \(\ell_0\). Do you get the same values as theETS()function?
my_ses_sse <- function(par, y) {
alpha <- par[1]
level <- par[2]
n <- length(y)
yhat <- numeric(n)
yhat[1] <- level
for (i in 2:n) {
yhat[i] <- alpha * y[i - 1] + (1 - alpha) * yhat[i - 1]
}
return(sum((y - yhat)^2))
}
optim(c(0.1, vic_pigs_vec[1]), my_ses_sse, y = vic_pigs_vec)$par[1] 3.220344e-01 1.005254e+05tidy(fit)# A tibble: 2 × 5
Animal State .model term estimate
<fct> <fct> <chr> <chr> <dbl>
1 Pigs Victoria ses alpha 0.322
2 Pigs Victoria ses l[0] 100647. Similar, but not identical estimates. This is due to different starting values being used.
Combine your previous two functions to produce a function that both finds the optimal values of \(\alpha\) and \(\ell_0\), and produces a forecast of the next observation in the series.
my_ses <- function(y) {
par <- optim(c(0.1, y[1]), my_ses_sse, y = y)$par
alpha <- par[1]
level <- par[2]
yhat <- numeric(length(y) + 1)
yhat[1] <- level
for (i in 2:(length(yhat))) {
yhat[i] <- alpha * y[i - 1] + (1 - alpha) * yhat[i - 1]
}
return(last(yhat))
}
my_ses(vic_pigs_vec)[1] 95186.69