16
$\begingroup$

I create an autoregressive process "from scratch" and I set the stochastic part (noise) equal to 0. In R:

Y <- vector() 

c = 0.1
phi = 0.9

Y[1] = 5

for (i in 2:100) {
Y[i] = c + phi*Y[i-1]
}

Then I ask R to fit an AR(1) process using the arima() function :

ar <- arima(Y, order = c(1,0,0))

It estimates the ar1 coefficient to be ar1 = 0.9989 with standard error 0.0015.

Why is R not finding ar1 = 0.9 (= phi) with overwhelming small standard error?

$\endgroup$
  • 2
    $\begingroup$ I've never seen an AR process without noise. You're turning the stochastic part of the AR equation in to constant which may be the reason you're getting undesirable results. $\endgroup$ – SOULed_Outt Dec 11 '19 at 20:35
25
$\begingroup$

You have two problems--and one of them is interesting.

Without a noise term, the series is no longer stationary. Its value is increasing asymptotically, but definitely, toward $1:$

Figure 1

ARIMA applies only to stationary models--and these data are obviously not from a stationary model. That's not terribly interesting. What is interesting is that the problem persists even with noise!

What, then, happens when we add just a tiny bit of noise?

Figure 2

It's still obviously not stationary--but the reason is that the initial values are inconsistent with everything that follows.

You need to remove a "burn-in period" during which the simulated values are starting to behave like the rest of the series will. Here's what this one looks like when we strip out the first $n_0=30$ values:

Figure 3

What does arima return?

Coefficients:
         ar1  intercept
      0.9074     0.9872
s.e.  0.0309     0.0088

$0.9074 \pm 0.0309$ is a great estimate of $\phi=0.9.$

I repeated this process 99 more times, producing 100 estimates of $\phi$ along with their standard errors. Here is a plot of those estimates and crude $90\%$ confidence limits (set at $1.645$ standard errors above and below the estimates):

Figure 4

The horizontal gray line is located at $\phi=0.9$ for reference. The red confidence intervals are those that do not overlap the reference: there are $12$ of them, indicating the confidence level is around $88\%,$ agreeing (within sampling error) with the intended value of $90\%.$ The horizontal black line is the average estimate. It's a little lower than $\phi,$ perhaps because after even $200$ time steps the series still isn't quite stationary. (One also doesn't expect the distribution of the estimate to be symmetric: $1$ is an important boundary and will cause the distribution to be skewed toward the smaller values.)

Here is the same study but with $2000$ time steps in each iteration:

Figure 5

The bias in the estimate has nearly disappeared.

Another solution is to start generating the series at its asymptotic mean value (equal to cnst/(1-phi) in the R code below). But that requires knowing the asymptote, which might be harder to come by in more complex models, so it's good to know about the technique of discarding the initial segment of a simulated series.

BTW, here's a reasonably efficient and compact way to generate these datasets:

phi <- 0.9    # AR(1) coefficient
n <- 200      # Total number of time steps after the initial value
cnst <- 0.1   # Intercept
sigma <- 0.01 # Error variance

Y <- Reduce(function(y, e) y * phi + e, rnorm(n, cnst, sigma), 0, accumulate=TRUE)
n0 <- which.max(abs(Y) >= quantile(abs(Y), 0.5)) # Estimate where Y levels off
Y <- Y[-seq_len(n0)]                             # Strip the initial values
plot(Y)                                          # LOOK at Y before doing anything else...
| cite | improve this answer | |
$\endgroup$
  • 2
    $\begingroup$ Thanks whuber - really appreciate it ! Thanks to everyone else as well. $\endgroup$ – Marc_Adrien Dec 12 '19 at 0:47
  • $\begingroup$ sorry, but can we have full replication (source code disclosed)? $\endgroup$ – Maximilian Dec 13 '19 at 16:27
  • $\begingroup$ @Max I no longer have all the code, but the essentials are already here: just apply arima to Y, extract its coefficients and standard errors, and wrap a loop around everything beginning with the creation of a random Y. (Comment out the plot command, though!) Plot the results. $\endgroup$ – whuber Dec 18 '19 at 20:59
5
$\begingroup$

The code you have created is not even generating any (pseudo) random outputs, let alone an AR(1) process. If you would like to generate the output of a stationary Gaussian AR(1) process, you can use the function below. This function generates exact output from the process, using the stationary marginal distribution of the process as the starting distribution; this means that the computation does not require any removal of "burn in" iterations. As such, this code should be computationally faster ---and statistically more exact--- than methods that anchor the process to a fixed starting value and then discard burn-in iterations.

GENERATE_NAR1 <- function(n, phi = 0, mu = 0, sigma = 1) {
  if (abs(phi) >= 1) { stop('Error: This is not a stationary process --- |phi| >= 1') }
  EE    <- rnorm(n, mean = 0, sd = sigma);
  YY    <- rep(0, n);
  YY[1] <- mu + EE[1]/sqrt(1-phi^2);
  for (t in 2:n) {
     YY[t] <- mu + phi*(YY[t-1]-mu) + EE[t]; }
  YY; }

Below I will generate a series of $n=10^5$ observable values and then we can use the arima function to fit it to an ARIMA model with specified orders.

set.seed(1);

phi   <- 0.9;
mu    <- 5;
sigma <- 4;
Y     <- GENERATE_NAR1(n = 10^5, phi, mu, sigma);
MODEL <- arima(Y, order = c(1,0,0), include.mean = TRUE);

MODEL;

    Call:
arima(x = Y, order = c(1, 0, 0), include.mean = TRUE)

Coefficients:
         ar1  intercept
      0.8978     4.9099
s.e.  0.0014     0.1242

sigma^2 estimated as 16.11:  log likelihood = -280874.1,  aic = 561754.2

As you can see, with this many data points, the arima function estimates the true parameters of the AR(1) model to within a reasonable level of accuracy --- the true values are within two standard errors of the estimates. A plot of the first $n = 1,000$ values in the time-series shows that it does not require you to discard any "burn-in" iterations.

plot(Y[1:1000], type = 'l', 
     main = paste0('Gaussian AR(1) time-series \n (phi = ',
                   phi, ', mu = ' , mu, ', sigma = ' , sigma, ')'),
     xlab = 'Time', ylab = 'Value');

enter image description here

| cite | improve this answer | |
$\endgroup$
3
$\begingroup$

As also pointed out in the comments, it's not an AR(p) process any more. arima function assumes the following and fits the coefficients accordingly: $$y_t=c+\phi y_{t-1}+\epsilon_t$$ It's quite normal that you don't get near the correct $\phi$. Also, after adding noise term, try to increase the sample size to get more confident estimates.

| cite | improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ Also, be careful because how the mean is handled in the arima() function can be tricky and not what you necessarily expect. So, look at the intercept and see what it gave for that. There are atleast 3 or 4 different arima functions in R. Stoffer, whose I think is at university of pittsburgh, wrote something a long while back about the pitfalls of using arima() in R and they might be still relevant ? If you can't find his link, let me know and I can find it. It's worth reading. $\endgroup$ – mlofton Dec 11 '19 at 20:54
  • 1
    $\begingroup$ Okay, thanks guys. I added some random noise, normally distributed with mean 0 and variance $=\sigma^2$, and increased the number of points in the time series. For small values of $\sigma$ the true value of the autoregressive parameter ar1 = 0.9 is still NOT within the interval of confidence of 1 sd. For larger values of $\sigma$ this "problem" goes away. I am disappointed because i still don't have an intuition why this is so. $\endgroup$ – Marc_Adrien Dec 11 '19 at 21:19
  • 1
    $\begingroup$ stat.pitt.edu/stoffer/tsa2/Rissues.htm this is useful, thanks mlofton $\endgroup$ – Marc_Adrien Dec 11 '19 at 21:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.