Confuso sobre el proceso autorregresivo AR(1)

Question

Creo un proceso autorregresivo "desde cero" y pongo la parte estocástica (ruido) igual a 0. En R :

Y <- vector() 

c = 0.1
phi = 0.9

Y[1] = 5

for (i in 2:100) {
Y[i] = c + phi*Y[i-1]
}

Entonces pregunto R para ajustar un proceso AR(1) utilizando el arima() función :

ar <- arima(Y, order = c(1,0,0))

Estima que el coeficiente ar1 es ar1 = 0,9989 con un error estándar de 0,0015.

¿Por qué es R ¿no encontrar ar1 = 0,9 (= phi) con un error estándar abrumadoramente pequeño?

Answer 1

Tienes dos problemas y uno de ellos es interesante.

Sin un término de ruido, la serie ya no es estacionaria. Su valor aumenta asintóticamente, pero definitivamente, hacia $1:$

ARIMA sólo se aplica a los modelos estacionarios, y estos datos obviamente no son de un modelo estacionario. Eso no es terriblemente interesante. Lo que es Lo interesante es que el problema persiste incluso con el ruido.

¿Qué ocurre entonces cuando añadimos un poco de ruido?

Es obvio que todavía no es estacionario pero la razón es que los valores iniciales son inconsistentes con todo lo que sigue.

Hay que eliminar un "periodo de rodaje" durante el cual los valores simulados empiezan a comportarse como lo hará el resto de la serie. Esto es lo que parece cuando quitamos el primer $n_0=30$ valores:

¿Qué hace arima ¿regresar?

Coefficients:
         ar1  intercept
      0.9074     0.9872
s.e.  0.0309     0.0088

$0.9074 \pm 0.0309$ es una gran estimación de $\phi=0.9.$

Repetí este proceso 99 veces más, produciendo 100 estimaciones de $\phi$ junto con sus errores estándar. Aquí hay un gráfico de esas estimaciones y de los errores estándar $90\%$ límites de confianza (fijados en $1.645$ errores estándar por encima y por debajo de las estimaciones):

La línea gris horizontal se sitúa en $\phi=0.9$ como referencia. Los intervalos de confianza rojos son los que no se superponen a la referencia: hay $12$ de ellos, lo que indica que el nivel de confianza es de alrededor de $88\%,$ que coincide (dentro del error de muestreo) con el valor previsto de $90\%.$ La línea negra horizontal es la estimación media. Es un poco más baja que $\phi,$ tal vez porque después de incluso $200$ pasos de tiempo la serie aún no es del todo estacionaria. (Tampoco se espera que la distribución de la estimación sea simétrica: $1$ es un límite importante y hará que la distribución esté sesgada hacia los valores más pequeños).

Aquí está el mismo estudio pero con $2000$ pasos de tiempo en cada iteración:

El sesgo en la estimación casi ha desaparecido.

Otra solución es empezar a generar la serie en su valor medio asintótico (igual a cnst/(1-phi) en el R código más abajo). Pero eso requiere conocer la asíntota, que puede ser más difícil de conseguir en modelos más complejos, por lo que es bueno conocer la técnica de descartar el segmento inicial de una serie simulada.

Por cierto, aquí hay una forma razonablemente eficiente y compacta de generar estos conjuntos de datos:

phi <- 0.9    # AR(1) coefficient
n <- 200      # Total number of time steps after the initial value
cnst <- 0.1   # Intercept
sigma <- 0.01 # Error variance

Y <- Reduce(function(y, e) y * phi + e, rnorm(n, cnst, sigma), 0, accumulate=TRUE)
n0 <- which.max(abs(Y) >= quantile(abs(Y), 0.5)) # Estimate where Y levels off
Y <- Y[-seq_len(n0)]                             # Strip the initial values
plot(Y)                                          # LOOK at Y before doing anything else...

Answer 2

El código que has creado ni siquiera está generando ninguna salida (pseudo) aleatoria, y mucho menos un proceso AR(1). Si desea generar la salida de un proceso gaussiano estacionario AR(1), puede utilizar la función que se indica a continuación. Esta función genera la salida exacta del proceso, utilizando la distribución marginal estacionaria del proceso como la distribución de partida; esto significa que el cálculo no requiere ninguna eliminación de iteraciones de "quemado". Como tal, este código debería ser computacionalmente más rápido -y estadísticamente más exacto- que los métodos que anclan el proceso a un valor inicial fijo y luego descartan las iteraciones de quemado.

GENERATE_NAR1 <- function(n, phi = 0, mu = 0, sigma = 1) {
  if (abs(phi) >= 1) { stop('Error: This is not a stationary process --- |phi| >= 1') }
  EE    <- rnorm(n, mean = 0, sd = sigma);
  YY    <- rep(0, n);
  YY[1] <- mu + EE[1]/sqrt(1-phi^2);
  for (t in 2:n) {
     YY[t] <- mu + phi*(YY[t-1]-mu) + EE[t]; }
  YY; }

A continuación voy a generar una serie de $n=10^5$ valores observables y luego podemos utilizar el arima para ajustarlo a un modelo ARIMA con órdenes especificados.

set.seed(1);

phi   <- 0.9;
mu    <- 5;
sigma <- 4;
Y     <- GENERATE_NAR1(n = 10^5, phi, mu, sigma);
MODEL <- arima(Y, order = c(1,0,0), include.mean = TRUE);

MODEL;

    Call:
arima(x = Y, order = c(1, 0, 0), include.mean = TRUE)

Coefficients:
         ar1  intercept
      0.8978     4.9099
s.e.  0.0014     0.1242

sigma^2 estimated as 16.11:  log likelihood = -280874.1,  aic = 561754.2

Como puede ver, con tantos puntos de datos, el arima estima los verdaderos parámetros del modelo AR(1) dentro de un nivel razonable de precisión --- los valores verdaderos están dentro de dos errores estándar de las estimaciones. Un gráfico de la primera $n = 1,000$ en la serie temporal muestra que no requiere descartar ninguna iteración de "quemado".

plot(Y[1:1000], type = 'l', 
     main = paste0('Gaussian AR(1) time-series \n (phi = ',
                   phi, ', mu = ' , mu, ', sigma = ' , sigma, ')'),
     xlab = 'Time', ylab = 'Value');

Answer 3

Como también se ha señalado en los comentarios, ya no es un proceso AR(p). arima asume lo siguiente y ajusta los coeficientes en consecuencia: $$y_t=c+\phi y_{t-1}+\epsilon_t$$ Es muy normal que no se acerque a la correcta $\phi$ . Además, después de añadir el término de ruido, intente aumentar el tamaño de la muestra para obtener estimaciones más seguras.