ar(1) simulation with non-zero mean

Question

I can't seem to find the correct way to simulate an AR(1) time series with a mean that is not zero. I need 53 data points, rho = .8, mean = 300.

However, arima.sim(list(order=c(1,0,0), ar=.8), n=53, mean=300, sd=21) gives me values in the 1500s. For example:

1480.099 1480.518 1501.794 1509.464 1499.965 1489.545 1482.367 1505.103 (and so on)

I have also tried arima.sim(n=52, model=list(ar=c(.8)), start.innov=300, n.start=1) but then it just counts down like this:

238.81775870 190.19203239 151.91292491 122.09682547 96.27074057 [6] 77.17105923 63.15148491 50.04211711 39.68465916 32.46837830 24.78357345 21.27437183 15.93486092 13.40199333 10.99762449 8.70208879 5.62264196 3.15086491 2.13809323 1.30009732

and I have tried arima.sim(list(order=c(1,0,0), ar=.8), n=53,sd=21) + 300 which seems to give a correct answer. For example:

280.6420 247.3219 292.4309 289.8923 261.5347 279.6198 290.6622 295.0501 264.4233 273.8532 261.9590 278.0217 300.6825 291.4469 291.5964 293.5710 285.0330 274.5732 285.2396 298.0211 319.9195 324.0424 342.2192 353.8149 and so on..

However, I am in doubt that this is doing the correct thing? Is it still auto-correlating on the correct number then?

Answer 1

Your last option is okay to get the desired mean, "mu". It generates data from the model:

(y[t] - mu) = phi * (y[t-1] - mu) + \epsilon[t], epsilon[t] ~ N(0, sigma=21), t=1,2,...,n.

Your first approach sets an intercept, "alpha", rather than a mean:

y[t] = alpha + phi * y[t-1] + epsilon[t].

Your second option sets the starting value y[0] equal to 300. As long as |phi|<1 the influence of this initial value will vanish after a few periods and will have no effect on the level of the series.

Edit

The value of the standard deviation that you observe in the simulated data is correct. Be aware that the variance of the AR(1) process, y[t], is not equal the variance of the innovations, epsilon[t]. The variance of the AR(1) process, sigma^2_y, can be obtained obtained as follows:

Var(y[t]) = Var(alpha) + phi^2 * Var(y[t-1]) + Var(epsilon[t])

As the process is stationary Var(y[t]) = Var(t[t-1]) which we call sigma^2_y. Thus, we get:

sigma^2_y = 0 + phi^2 * sigma^2_y + sigma^2_epsilon sigma^2_y = sigma^2_epsilon / (1 - phi^2)

For the values of the parameters that you are using you have:

sigma_y = sqrt(21^2 / (1 - 0.8^2)) = 35.

Answer 2

Use the rGARMA function in the ts.extend package

You can generate random vectors from any stationary Gaussian ARMA model using the ts.extend package. This package generates random vectors directly form the multivariate normal distribution using the computed autocorrelation matrix for the random vector, so it gives random vectors from the exact distribution and does not require "burn-in" iterations. Here is an example of generating multiple independent time-series vectors all from an AR(1) model.

#Load the package
library(ts.extend)

#Set parameters
MEAN     <- 300
ERRORVAR <- 21^2
AR       <- 0.8
m        <- 53

#Generate n = 16 random vectors from this model
set.seed(1)
SERIES <- rGARMA(n = 16, m = m, mean = MEAN, ar = AR, errorvar = ERRORVAR)

#Plot the series using ggplot2 graphics
library(ggplot2)
plot(SERIES)

As you can see, the generated time-series vectors in this plot use the appropriate mean and error variance that were specified in the inputs.