Karsten W. - 8 months ago 49

R Question

In a recent online course on macroeconomic forecasting, there was an excercise to model

`y(t) = 3.0 + 0.55 y(t-1) + e(t)`

where

`e(t)`

`et <- c(-1.2138662, -0.2854597, 0.5902700, 0.8285463, -0.9954260, -0.3716332)`

Now I tried to do this in R (the course used EViews), however I do not come to the solution given there: 5.648 for the 5th element.

I tried (similar to this blogpost):

`y <- rep(NA,6)`

y[1] <- 0

y[2] <- 3 + 0.55*y[1]+et[1]

y[3] <- 3 + 0.55*y[2]+et[2]

y[4] <- 3 + 0.55*y[3]+et[3]

y[5] <- 3 + 0.55*y[4]+et[4]

y[6] <- 3 + 0.55*y[5]+et[5]

and then

`y <- rep(NA,6)`

y[1] <- et[1]

y[2] <- 3 + 0.55*y[1]+et[2]

y[3] <- 3 + 0.55*y[2]+et[3]

y[4] <- 3 + 0.55*y[3]+et[4]

y[5] <- 3 + 0.55*y[4]+et[5]

y[6] <- 3 + 0.55*y[5]+et[6]

and then

`arima.sim(list(order=c(1,0,0), ar=0.55), n=6, innov=head(et,6)+3)`

However all three methods give different results. I wonder why this is, I am afraid I did not understand something basic.

Answer Source

`arima.sim`

has a "burn-in" period in order to reach stationarity. Argument `n.start`

sets the length of this period, and it must has length `ar`

+ `ma`

. For an `AR(1)`

process, `n.start`

must be at least 1. Note, if you don't specify `n.start`

, `arima.sim`

will automatically compute a reasonable one for you.

Given `n.start`

(either a specified one or auto-computed one), `arima.sim`

then randomly sample `n.start`

innovations for "burn-in". It is right because of this random procedure, you will get different result from one run to another:

```
## innovations of length 6
et <- c(-1.2138662, -0.2854597, 0.5902700, 0.8285463, -0.9954260, -0.3716332)
set.seed(0)
arima.sim(list(order = c(1,0,0), ar = 0.55), n = 6, innov = et, n.start = 1)
#[1] -0.5192413 -0.5710424 0.2761967 0.9804545 -0.4561760 -0.6225300
set.seed(1)
arima.sim(list(order = c(1,0,0), ar = 0.55), n = 6, innov = et, n.start = 1)
# [1] -1.55841580 -1.14258839 -0.03815361 0.80756181 -0.55126700 -0.67483005
```

We can eliminate such randomness, by providing a set of starting innovations using `start.innov`

argument. Note, this vector must has length `n.start`

:

```
## fixing starting innovation at 0; i.e., `y[0] = 0`, so `y[1] = et[1]`.
arima.sim(list(order = c(1,0,0), ar = 0.55), n = 6, innov = et, n.start = 1,
start.innov = 0)
# [1] -1.21386620 -0.95308611 0.06607264 0.86488625 -0.51973856 -0.65748941
```

Now, compare this last `arima.sim`

with the following:

```
y <- rep(NA, 6)
y[1] <- et[1]
y[2] <- 0.55 * y[1] + et[2]
y[3] <- 0.55 * y[2] + et[3]
y[4] <- 0.55 * y[3] + et[4]
y[5] <- 0.55 * y[4] + et[5]
y[6] <- 0.55 * y[5] + et[6]
# [1] -1.21386620 -0.95308611 0.06607264 0.86488625 -0.51973856 -0.65748941
```

We have everything completely reproducible.

The above generates zero-mean time series; if you want mean at 3, shift the result by `+3`

yourself:

```
y <- y + 3
```