 # Question of the prediction interval of meanf function

I work for an auditing firm and perform time series analysis with reference to the online training("Forecasting Using R" at DataCamp) and the books.

In order to deepen my understanding, I have recalculated various methods myself, but I cannot reproduce the prediction interval of the meanf function.
Perhaps my formula is wrong, but it would be helpful if someone could teach me the correct formula.

library(fpp2)
library(zoo)

# Create a data

mydata <- c(10, 20, 30, 40)

# Create a ts object called myts

myts <- ts(mydata, start = c(2019, 1), frequency = 4)

autoplot(myts)

# meanf

fc_mean <- meanf(myts, h = 2)
autoplot(fc_mean) +
scale_x_yearqtr(format = "%Y-%q") +
labs(x = "", y = "")

# SD of Residuals

sd_resid <- fc_mean\$model\$sd

# Value of Point forecast

pf <- fc_mean\$mean

# Culculate the upper of 95% prediction interval

(ans1 <- pf + sd_resid * qnorm(0.975))

# the upper of 95% prediction interval of the model

(ans2 <- fc_mean\$upper)

# Not identical!

identical(ans1, ans2)

Referred here by Forecasting: Principles and Practice, by Rob J Hyndman and George Athanasopoulos

The `meanf()` function uses a t-distribution for producing prediction intervals. It also accounts for the parameter estimation error. Try this:

``````library(fpp2)
mydata <- c(10, 20, 30, 40)
myts <- ts(mydata, start = c(2019, 1), frequency = 4)
fc_mean <- meanf(myts, h = 2)
sd_resid <- fc_mean\$model\$sd
pf <- fc_mean\$mean
n <- length(mydata)
(ans1 <- pf + sd_resid * qt(0.975,n-1) * sqrt(1+1/n))
(ans2 <- fc_mean\$upper)
identical(ans1, ans2)
``````

I understand very well.

Another question came to mind, but in the case of the ses method, it seems that the normal distribution is used and the parameter estimation error is not taken into account.
Is the same consideration as the meanf method unnecessary?

``````library(fpp2)
mydata <- c(10, 30, 20, 40)
myts <- ts(mydata, start = c(2019, 1), frequency = 4)
fc_ses <- ses(myts, h = 2)
sd_resid <- sqrt(fc_ses\$model\$sigma2)
pf <- fc_ses\$mean
(ans1 <- pf + sd_resid * qnorm(0.975))
(ans2 <- fc_ses\$upper)
identical(round(ans1,5), round(ans2,5))
``````

It would probably be better to use qt rather than qnorm in all these functions.
Taking parameter estimation into account is very difficult except for some simple methods such as linear models. So we usually ignore it.

1 Like

This question was solved and I was refreshed!
(and I learned a lot!!)

If your question's been answered (even by you!), would you mind choosing a solution? It helps other people see which questions still need help, or find solutions if they have similar problems. Here’s how to do it:

I have another question.
I cannot reproduce the prediction interval of the holt-winters' multiplicative function.
(Although, I can do that of the holt-winters' additive function.)
What is the cause?

library(fpp2)

# Create a data

mydata <- c(10,30,20,40,20,50,40,60,50,70,60,90)

# Create a ts object called myts

myts <- ts(mydata, start = c(2019, 1), frequency = 4)

fc_hwa <- hw(myts, seasonal = "additive", h = 4)
autoplot(fc_hwa) # Create an autoplot
sd_resid <- sqrt(fc_hwa\$model\$sigma2) # SD of Residuals
pf <- fc_hwa\$mean # Value of Point forecast
(ans1 <- pf + sd_resid * qnorm(0.975)) # Culculate the upper of 95% prediction interval
(ans2 <- as.numeric(fc_hwa\$upper[1, 2])) # The upper of 95% prediction interval of the model
identical(ans1, ans2) # Identical

# Holt-Winters' multiplicative

fc_hwm <- hw(myts, seasonal = "multiplicative", h = 4)
autoplot(fc_hwm)
sd_resid <- sqrt(fc_hwm\$model\$sigma2)
pf <- fc_hwm\$mean
(ans1 <- pf + sd_resid * qnorm(0.975)) # Too low!
(ans2 <- as.numeric(fc_hwm\$upper[1, 2]))
identical(ans1, ans2) # Not identical!

Thank you for letting me know!

With multiplicative errors, the one-step forecast error variance is yhat^2*sigma^2

``````> (ans1 <- pf + pf * sd_resid * qnorm(0.975))
 128.0357
> (ans2 <- as.numeric(fc_hwm\$upper[1, 2]))
 128.0357
``````

I see.