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:
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.
- make 21 normal random numbers with mu=100, sigma=30
- calculate the error between the mean of first 20 numbers and the last number.
- 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))