Question of the prediction interval of meanf function

[Takashi]

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.

Packages

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)

Create an autoplot

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[1]

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[3])

Not identical!

identical(ans1, ans2)

---

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

[robjhyndman]

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[1]
n <- length(mydata)
(ans1 <- pf + sd_resid * qt(0.975,n-1) * sqrt(1+1/n))
(ans2 <- fc_mean$upper[3])
identical(ans1, ans2)

[Takashi]

Thank you for your reply!!
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[1]
(ans1 <- pf + sd_resid * qnorm(0.975))
(ans2 <- fc_ses$upper[3])
identical(round(ans1,5), round(ans2,5))

[robjhyndman]

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.

[Takashi]

This question was solved and I was refreshed!
(and I learned a lot!!)
Thank you very much for your reply!!

[andresrcs]

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:

FAQ: How do I mark a solution?

If your question has been answered, don't forget to mark the solution! How do I mark a solution? Find the reply you want to mark as the solution and look for the row of small gray icons at the bottom of that reply. Click the one that looks like a box with a checkmark in it: [image] Hovering over the mark solution button shows the label, "Select if this reply solves the problem". If you don't see the mark solution button, try clicking the three dots button ( ••• ) to expand the full set of options. When a solution is chosen, the icon turns green and the hover label changes to: "Unselect if this reply no longer solves the problem". Success! [solution_reply_author] …

[Takashi]

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?

Packages

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)

Holt-Winters' additive

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[1] # 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[1]
(ans1 <- pf + sd_resid * qnorm(0.975)) # Too low!
(ans2 <- as.numeric(fc_hwm$upper[1, 2]))
identical(ans1, ans2) # Not identical!

[Takashi]

Thank you for letting me know!

[robjhyndman]

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

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

[Takashi]

I see.
Thank you very much for your reply!!