how to use arima to do mean model

Question

I am learning arima by this site:

http://people.duke.edu/~rnau/411home.htm

and I want to get the same result as following notes:

http://people.duke.edu/~rnau/Review_of_basic_statistics_and_the_mean_model_for_forecasting--Robert_Nau.pdf

I was thinking that an arima with order [0, 0, 0] is the mean model, but the results is different from the notes, Here is the code:

require(forecast);
x <- c(114, 126, 123, 112, 68, 116, 50, 108, 163, 79,
      67, 98, 131, 83, 56, 109, 81, 61, 90, 92);
m <- arima(x, order=c(0, 0, 0));
print(m);
print(forecast(m, 1));
print(predict(m)$se);

the output:

> print(m);
Series: x 
ARIMA(0,0,0) with non-zero mean 

Coefficients:
      intercept
        96.3500
s.e.     6.3124

sigma^2 estimated as 796.9:  log likelihood=-95.19
AIC=194.37   AICc=195.08   BIC=196.36

> print(forecast(m, 1));
   Point Forecast    Lo 80    Hi 80   Lo 95    Hi 95
21          96.35 60.17192 132.5281 41.0204 151.6796

> print(predict(m)$se);
Time Series:
Start = 21 
End = 21 
Frequency = 1 
[1] 28.2299

but the results in the notes are:

SE_fcst = 29.68  (R result: 28.2299)
95% confidence intervals = 34, 158  (R result: 41, 152)

Where am I wrong?

edit

I do the simulation with random numbers, and the result is the same as the notes.

1. make 21 normal random numbers with mu=100, sigma=30
2. calculate the error between the mean of first 20 numbers and the last number.
3. repeat 1 & 2 for 100000 times

Here is the python code that to do the simulation:

import numpy as np
N = 1000000
n = 20
x = np.random.normal(100, 30, (N, n))
p = np.mean(x, axis=1)
nx = np.random.normal(100, 30, N)

err = p - nx
print (err**2).mean()**0.5

s = np.std(x, axis=1, ddof=1)
SE_mean = s / n**0.5
print (s**2 + SE_mean**2).mean()**0.5

the output is:

30.7480552149 (the real standard error of forecast)
30.7375157915 (the estimated standard error of forecast by sqrt(s**2 + SE_mean**2))

Answer 1

The first difference that I notice comes from the way in which the standard deviation of the residuals is obtained. In arima, it is obtained as follows (given that the residuals have zero mean):

n <- length(x)
sqrt(sum(residuals(m)^2) / n)
#[1] 28.2299

Sometimes, the fact that the mean is unknown and must be estimated (in this case as the sample mean of the residuals which is zero) is reflected by removing one degree of freedom in the denominator. See for example this link for further discussion.

sqrt(sum(residuals(m)^2) / (n-1))
#[1] 28.96327

This latter expression is used in the reference that you give. You can do some preliminary calculations to choose the factor that you wish when computing the standard errors of the forecasts:

s2 <- sqrt(sum(residuals(m)^2) / (n-1))
KFvar <- KalmanForecast(n.ahead = 1, m$model)[[2]]
KFvar * s2
#[1] 28.96327

Edit

The notes that you give makes a different calculus. Two sources of uncertainty are added up, the variability intrinsic to the data and the error of the estimated mean. The standard errors of the forecasts based on arima are obtained as shown in this post; the errors of the parameter estimates are not accounted in those calculations.