# Simulating the error term with an AR(2) model

I'm working on some simulations in R where I have a time series in the format $$y_t=\mu_t+e_t$$ where $$\mu_t$$ is a trend function that I have defined as $$\mu_t=5*t+4*sin(2t)$$ and where $$t$$ goes from $$\{ 0,4\pi \}$$ with $$n=100$$ datapoints. Generally $$e_t \sim N(0, \sigma^2)$$ however in this case of these simulations I wish to replace $$e_t$$ with an AR(2) process. The reasoning behind this is that normally if $$e_t \sim N(0, \sigma^2)$$ then we assume that each $$e_t$$ is i.i.d however by limiting the error term to a known process the i.i.d is dropped from our assumptions, which is what I am trying to simulate here and I would like to know how I would go about changing the error term to be AR(2).

I've posted a short R-code snippet that shows the function defined above with a normal $$e_t$$ term which would be replaced with an AR(2) process.

n_sim <- 100 # number of data points
t_sim <- seq(0,4*pi,,100)
c_norm_sim <- rnorm(n_sim)
y2_sim <-5*t_sim+4*sin(2*t_sim)+2*c_norm_sim # Gaussian/normal error
plot(t_sim, y2_sim, t="l", ylim=range(y1_sim,y2_sim)*c(1,1.2))



Question: How do I simulate a vector of error values from the AR(2) process?

• What is your question? Commented Jun 30, 2022 at 20:00
• Hello Hardy, my question is how would I estimate the error terms as an AR(2) process. Commented Jun 30, 2022 at 20:02
• Estimate or simulate? Commented Jul 1, 2022 at 6:56
• Simulate the AR(2) process first (using the definition of an auto-regressive AR(2) process) and then create your sequence $(y_t)$. Commented Jul 1, 2022 at 9:11
• Sorry, I meant to say simulate. Commented Jul 1, 2022 at 14:02

#### Use the rGARMA function in the ts.extend package

You can simulate an $$\text{AR}(2)$$ model using the rGARMA function in the ts.extend package (see O'Neill 2021 for details). This package allows you to simulate from the marginal or conditional distribution from any stationary Gaussian ARMA process.$$^\dagger$$ To use the function you need to stipulate the parameters of your model (and any conditioning values you want to use).

Below I give an example of some code where you first stipulate the parameters of the model and then simulate $$n=12$$ vectors of $$m=100$$ values from the process. (By default the function uses a process with zero mean and unit error variance, but you can change these values if you want.) The generated object ERRORS is a $$12 \times 100$$ matrix showing the simulated time-series vectors. If you just want to simulate a single time-series vector you can set $$n=1$$ in the simulation.

#Load the package
library(ts.extend)

#Set the model parameters
AR    <- c(0.8, -0.1)
SIGMA <- 5

#Simulate and plot a vector of 100 values from an AR(2) process
set.seed(1)
ERRORS <- rGARMA(n = 12, m = 100, ar = AR, errorvar = SIGMA)
plot(ERRORS)


The plot generated by the package shows each of the simulated time-series vectors; the points in the background are the points from all $$n=12$$ vectors. (If you want to get rid of these then just set background = FALSE when you call the plot.)

$$^\dagger$$ The "GARMA" in rGARMA refers to the Gaussian Auto-Regressive Moving Average model.

• But how does rGARMA differ from the more standard arima.sim function in the stats package? Commented Jul 4, 2022 at 8:19
• There are a couple of key differences --- arima.sim uses a burn-in method, whereas rGARMA uses a direct generation method using the joint distribution. Because it uses this method, the latter function also allows you to specify conditioning values whereas the former does not, and it does not require any burn-in. Nevertheless, another answer pointing out how to use arima.sim for this question would round things out nicely.
– Ben
Commented Jul 4, 2022 at 8:28