How to simulate an AR(1) process with arima.sim and an estimated model?

Question

I want to do the following two steps:

1. Based on a given time series, I want to calibrate an AR(1) process, i.e. I want to estimate the parameters.
2. Based on the estimated parameters, I want to simulate an AR(1) processes.

Here was my approach:

set.seed(123)
#Just generate random AR(1) time series; based on this, I want to estimate the parameters
ts_AR <- arima.sim(n=10000, list(ar=c(0.5)))
#1. Estimate parameters with arima()
model_AR <- arima(ts_AR, order=c(1,0,0))
#Looks actually good
model_AR
Series: ts_AR 
ARIMA(1,0,0) with non-zero mean 

Coefficients:
  ar1  intercept
0.4891    -0.0044
s.e.  0.0087     0.0195

sigma^2 estimated as 0.9974:  log likelihood=-14176.35
AIC=28358.69   AICc=28358.69   BIC=28380.32  

#2. Simulate based on model
arima.sim(model=model_AR, n = 100)
Error in arima.sim(model = model_AR, n = 100) : 
  'ar' part of model is not stationary

I'm not the biggest time-series expert, but I'm pretty sure that an AR(1) process with a persistence parameter of below one should result in a stationary model. However, the error message tells me somethings different. So do I do something stupid here? If so, why and what should I do to simulate the AR(1) process based on my estimated parameters. Or can't you just pass the output of arima as the model input into arima.sim? Then, however, I don't understand how I get such an error message...I would expect something like "model input cannot be read. It should be something like ..."

Answer 1

It's not the clearest interface in the world, but the model argument is meant to be a list giving the ARMA order, not an actual arima model.

arima.sim(model=as.list(coef(model_AR)), n=100)

This will create a simulated series with AR coefficient .489 as estimated from your starting data. Note that the intercept is ignored.

Answer 2

I don't think you are using the right approach since there's uncertainty about your coefficient estimate. The best way to achieve what you want in a proper way is to incorporate uncertainty in the generation process, there are probably parametric way to do that but I think bootstrap can be handy here.

Lets generate the AR process first

set.seed(123)
ts_AR <- arima.sim(n = 10000, list(ar = 0.5))

We'll define two helper functions that will used in the boostrap. The first one generate the statistics we need (here the coef of the AR process and the actual time series) and the second function implement our resampling scheme (it'll be based on residuals)

ar_fun <- function(ts) c(ar = coef(arima(ts, order = c(1, 0, 0),
                                       include.mean = FALSE)), ts = ts)

ar_sim <- function(res, n.sim, ran.args) {
    rg <- function(n, res) sample(res, n, replace = TRUE)
    ts <- ran.args$ts
    model <- ran.args$model
    arima.sim(model = model, n = n.sim,
              rand.gen = rg, res = c(res))
}

Now we can start our simulation

ar_fit <- arima(ts_AR, order = c(1, 0, 0), include.mean = FALSE)
ts_res <- residuals(ar_fit)
ts_res <- ts_res - mean(ts_res)
ar_model <- list(ar = coef(ar_fit))

require(boot)
set.seed(1)
ar_boot <- tsboot(ts_res, ar_fun,
                   R = 99, sim = "model",
                   n.sim = 100, orig.t = FALSE,
                   ran.gen = ar_sim,
                   ran.args = list(ts = ts_AR, model = ar_model))

If you want to get all the coefficient generated and the associated time series

coefmat <- apply(ar_boot$t, 1, "[", 1)
seriesmat <- apply(ar_boot$t, 1, "[", -1)

You can get more details in help file of tsboot and in Bootstrap Methods and Their Application, chap 8.