how to use arima to do mean model

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))


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.

• Thanks for the answer, I added some simulation code to the question, It's seems that SE_fcst = sqrt(s^2 + SE_mean^2) in the notes is correct. I am not sure whether ARIMA with order (0, 0, 0) is the same as the mean model in the notes or not? – HYRY Oct 9 '14 at 0:38