R: Simular AR(2) estableciendo la varianza de la serie como valores iniciales

Question

1. Estoy tratando de simular ARIMA(2,0,0) que siguen esta ecuación:

2. Yo que establecer los dos primeros valores iniciales para ser la varianza de `ARIMA(2,0,0)que son:

Tengo esto r código para representar mi interés:

n <- 50
phi1 <- 0.5 
phi2 <- -0.7
e <- rnorm(n, mean=0, sd=1)
x1 <- sqrt((e[1])^2/(1-phi1^2))
x2 <- sqrt(((1-phi2^2)*e[1]^2)/((1-phi2)*((1-phi2)^2-phi1^2)))
set.seed(123)
ar1 <- arima.sim(n, model=list(ar=c(phi1, phi2)), start.innov=c(x1,x2), n.start=2, innov=c(0,0,e[-2])) 

Quiero que revises mi r código para ver si me estoy perdiendo algo.

Answer 1

Con respecto a un comentario anterior, estoy asumiendo que $\sigma^2_\varepsilon \equiv \operatorname{Var} (\varepsilon_t)$ . Si eso es cierto, entonces he visto muchos casos en los que se asume que $\sigma^2_\varepsilon = 1$ y su código debe actualizar su código a

x1 <- sqrt(1/(1-phi1^2))
x2 <- sqrt((1-phi2^2)/((1-phi2)*((1-phi2)^2-phi1^2)))

No tengo ni idea de cómo el arima.sim puede manejar la especificación de valores iniciales. He investigado un poco y he encontrado algunos notas del Dr. Henrik Madsen para la simulación de procesos AR . En las páginas 2 y 3 dice que se puede utilizar el filter para simular el proceso AR.

Así que tienes el proceso $x_t = 0.5 x_{t-1} - 0.7 x_{t-2} + \varepsilon_t$ . A partir de esta ecuación recursiva se tienen las siguientes dos ecuaciones para $x_1$ y $x_2$ : \begin{align} x_1 &= 0.5 x_0 - 0.7 x_{-1} + \varepsilon_1 \\ x_2 &= 0.5 x_1 - 0.7 x_{0} + \varepsilon_2 \end{align} Puedes ver que dependen de los valores $x_0$ y $x_{-1}$ que no existen, por lo que el filter los trata como cero y genera $$x_1 = \varepsilon_1$$ y $$x_2 = 0.5 \varepsilon_1 + \varepsilon_2$$

Así que puedes "engañarlo" generando un vector para el proceso de ruido, digamos $\left(\varepsilon_1, \varepsilon_2, \ldots, \varepsilon_n \right)$ y modificar la primera y la segunda componente del vector de la siguiente manera: $$\left(x_1,\: x_2-0.5x_1\:, \ldots, \varepsilon_n \right)$$ .

Por ejemplo, si $x_1 = 1$ , $x_2 = -2$ y $\varepsilon_t \overset{\text{iid}}{\sim} \mathcal{N} (0, 1)$ entonces puede utilizar el siguiente código:

n <- 50                       # length of series
x_1 <- 1                      # first initial value
x_2 <- -2                     # second initial value
phi <- c(0.5, -0.7)           # vector of autoregressive coefficients
e <- rnorm(n)                 # vector from standard normal noise process
e[1] <- x_1                   # updating so the first value of the series is x_1
e[2] <- x_2 - phi[1]*x_1      # updating so the second value of the series is x_2
ar2 <- filter(e, filter = phi, method = "recursive")
ar2

# output from code above
# [1]  1.00000000 -2.00000000 -1.74287046  1.89716706  1.94282186  1.15986460
# [7] -2.32879581 -1.39168937  1.05816662  1.71920744  0.49852657 -1.45650538
# ...

---

He conseguido reproducir la salida anterior utilizando el arima.sim función. Pero debo admitir que todavía no entiendo realmente cómo start.innov se maneja en la función.

ar2.sim <- arima.sim(n, model=list(ar=c(phi[1], phi[2])), 
            start.innov=c(0, 0), n.start=2, 
            innov=c(x_1,x_2 - phi[1]*x_1 ,e[-c(1,2)]), sd = 1)
ar2.sim

# output from code above
# [1]  1.00000000 -2.00000000 -1.74287046  1.89716706  1.94282186  1.15986460
# [7] -2.32879581 -1.39168937  1.05816662  1.71920744  0.49852657 -1.45650538
# ...