Fórmula ACF y PACF

Question

Quiero crear un código para trazar ACF y PACF a partir de datos de series temporales. Al igual que este gráfico generado a partir de Minitab (abajo).

He intentado buscar la fórmula, pero sigo sin entenderla bien. ¿Le importaría decirme la fórmula y cómo utilizarla, por favor? ¿Cuál es la línea roja horizontal en los gráficos ACF y PACF? ¿Cuál es la fórmula?

Gracias,

Answer 1

Autocorrelaciones

La correlación entre dos variables $y_1, y_2$ se define como:

$$ \rho = \frac{\hbox{E}\left[(y_1-\mu_1)(y_2-\mu_2)\right]}{\sigma_1 \sigma_2} = \frac{\hbox{Cov}(y_1, y_2)}{\sigma_1 \sigma_2} \,, $$

donde E es el operador de expectativa, $\mu_1$ y $\mu_2$ son las medias, respectivamente, de $y_1$ y $y_2$ y $\sigma_1, \sigma_2$ son sus desviaciones típicas.

En el contexto de una única variable, es decir auto -correlación, $y_1$ es la serie original y $y_2$ es una versión retardada de la misma. Sobre la definición anterior, autocorrelaciones muestrales de orden $k=0,1,2,...$ puede obtenerse calculando el siguiente expresión con la serie observada $y_t$ , $t=1,2,...,n$ :

$$ \rho(k) = \frac{\frac{1}{n-k}\sum_{t=k+1}^n (y_t - \bar{y})(y_{t-k} - \bar{y})}{ \sqrt{\frac{1}{n}\sum_{t=1}^n (y_t - \bar{y})^2}\sqrt{\frac{1}{n-k}\sum_{t=k+1}^n (y_{t-k} - \bar{y})^2}} \,, $$

donde $\bar{y}$ es la media muestral de los datos.

Autocorrelaciones parciales

Las autocorrelaciones parciales miden la dependencia lineal de una variable tras eliminar el efecto de otra(s) variable(s) que afectan a ambas variables. Por ejemplo, la autocorrelación parcial del orden mide el efecto (dependencia lineal) de $y_{t-2}$ en $y_t$ tras eliminar el efecto de $y_{t-1}$ en ambos $y_t$ y $y_{t-2}$ .

Cada autocorrelación parcial podría obtenerse como una serie de regresiones de la forma:

$$ \tilde{y}_t = \phi_{21} \tilde{y}_{t-1} + \phi_{22} \tilde{y}_{t-2} + e_t \,, $$

donde $\tilde{y}_t$ es la serie original menos la media muestral, $y_t - \bar{y}$ . La estimación de $\phi_{22}$ dará el valor de la autocorrelación parcial de orden 2. Ampliando la regresión con $k$ rezagos adicionales, la estimación del último término dará la autocorrelación parcial de orden $k$ .

Una forma alternativa de calcular las autocorrelaciones parciales de la muestra es resolviendo el siguiente sistema para cada orden $k$ :

\begin{eqnarray} \left(\begin{array}{cccc} \rho(0) & \rho(1) & \cdots & \rho(k-1) \\ \rho(1) & \rho(0) & \cdots & \rho(k-2) \\ \vdots & \vdots & \vdots & \vdots \\ \rho(k-1) & \rho(k-2) & \cdots & \rho(0) \\ \end{array} \(derecha) \izquierda( \begin{array}{c} \phi_{k1} \\ \phi_{k2} \\ \vdots \\ \phi_{kk} \\ \end{array} \right) = \izquierda( \begin{array}{c} \rho(1) \\ \rho(2) \\ \vdots \\ \rho(k) \\ \end{array} \derecha) \fin

donde $\rho(\cdot)$ son las autocorrelaciones de la muestra. Esta correspondencia entre las autocorrelaciones muestrales y las autocorrelaciones parciales se denomina Recursividad Durbin-Levinson . Este enfoque es relativamente fácil de aplicar a efectos ilustrativos. Por ejemplo, en el software R, podemos obtener la autocorrelación parcial de orden 5 de la siguiente manera:

# sample data
x <- diff(AirPassengers)
# autocorrelations
sacf <- acf(x, lag.max = 10, plot = FALSE)$acf[,,1]
# solve the system of equations
res1 <- solve(toeplitz(sacf[1:5]), sacf[2:6])
res1
# [1]  0.29992688 -0.18784728 -0.08468517 -0.22463189  0.01008379
# benchmark result
res2 <- pacf(x, lag.max = 5, plot = FALSE)$acf[,,1]
res2
# [1]  0.30285526 -0.21344644 -0.16044680 -0.22163003  0.01008379
all.equal(res1[5], res2[5])
# [1] TRUE

Bandas de confianza

Las bandas de confianza pueden calcularse como el valor de las autocorrelaciones de la muestra $\pm \frac{z_{1-\alpha/2}}{\sqrt{n}}$ donde $z_{1-\alpha/2}$ es el cuantil $1-\alpha/2$ en la distribución gaussiana, por ejemplo 1,96 para bandas de confianza del 95%.

A veces se utilizan bandas de confianza que aumentan a medida que aumenta el orden. En este caso, las bandas pueden definirse como $\pm z_{1-\alpha/2}\sqrt{\frac{1}{n}\left(1 + 2\sum_{i=1}^k \rho(i)^2\right)}$ .

Answer 2

"Quiero crear un código para graficar ACF y PACF de series de tiempo series temporales".

Aunque el PO es un poco vago, posiblemente esté más orientado a una formulación de codificación tipo "receta" que a una formulación de modelo de álgebra lineal.

---

En ACF es bastante sencillo: tenemos una serie temporal, y básicamente hacemos múltiples "copias" (como en "copiar y pegar") de ella, entendiendo que cada copia va a estar desplazada por una entrada de la copia anterior, porque los datos iniciales contienen $t$ puntos de datos, mientras que la longitud de la serie temporal anterior (que excluye el último punto de datos) es de sólo $t-1$ . Podemos hacer prácticamente tantas copias como filas haya. Cada copia se correlaciona con la original, teniendo en cuenta que necesitamos longitudes idénticas y, para ello, tendremos que ir recortando la cola de la serie de datos inicial para que sean comparables. Por ejemplo, para correlacionar los datos iniciales con $ts_{t-3}$ tendremos que deshacernos del último $3$ puntos de datos de la serie temporal original (el primer $3$ cronológicamente).

Por ejemplo:

Crearemos una serie temporal con un patrón sinusoidal cíclico superpuesto a una línea de tendencia y ruido, y trazaremos el ACF generado por R. Obtuve este ejemplo de un puesto en línea de Christoph Scherber, y sólo le añadí el ruido:

x=seq(pi, 10 * pi, 0.1)
y = 0.1 * x + sin(x) + rnorm(x)
y = ts(y, start=1800)

Normalmente, tendríamos que comprobar la estacionariedad de los datos (o simplemente mirar el gráfico anterior), pero sabemos que existe una tendencia, así que nos saltamos esta parte y pasamos directamente al paso de eliminación de la tendencia:

model=lm(y ~ I(1801:2083))
st.y = y - predict(model)

Ahora estamos listos para tomar esta serie temporal generando primero la ACF con la función acf() en R y comparando los resultados con el bucle improvisado que he creado:

ACF = 0                  # Starting an empty vector to capture the auto-correlations.
ACF[1] = cor(st.y, st.y) # The first entry in the ACF is the correlation with itself (1).
for(i in 1:30){          # Took 30 points to parallel the output of `acf()`
  lag = st.y[-c(1:i)]    # Introducing lags in the stationary ts.
  clipped.y = st.y[1:length(lag)]    # Compensating by reducing length of ts.
  ACF[i + 1] = cor(clipped.y, lag)   # Storing each correlation.
}
acf(st.y)                            # Plotting the built-in function (left)
plot(ACF, type="h", main="ACF Manual calculation"); abline(h = 0) # and my results (right).

---

DE ACUERDO. Eso fue un éxito. A la PACF . Mucho más complicado de hackear... La idea aquí es clonar de nuevo el ts inicial un montón de veces, y luego seleccionar múltiples puntos temporales. Sin embargo, en lugar de simplemente correlacionar con la serie temporal inicial, juntamos todos los retardos intermedios y realizamos un análisis de regresión, de forma que la varianza explicada por los puntos temporales anteriores pueda excluirse (controlarse). Por ejemplo, si nos centramos en el PACF que finaliza en el momento $ts_{t-4}$ mantenemos $ts_t$ , $ts_{t-1}$ , $ts_{t-2}$ y $ts_{t-3}$ así como $ts_{t-4}$ y hacemos una regresión $ts_t \sim ts_{t-1} + ts_{t-2} + ts_{t-3}+ts_{t-4}$ a través de el origen y manteniendo sólo el coeficiente de $ts_{t-4}$ :

PACF = 0          # Starting up an empty storage vector.
for(j in 2:25){   # Picked up 25 lag points to parallel R `pacf()` output.
  cols = j        
  rows = length(st.y) - j + 1 # To end up with equal length vectors we clip.

  lag = matrix(0, rows, j)    # The storage matrix for different groups of lagged vectors.

for(i in 1:cols){
  lag[ ,i] = st.y[i : (i + rows - 1)]  #Clipping progressively to get lagged ts's.
}
  lag = as.data.frame(lag)
  fit = lm(lag$V1 ~ . - 1, data = lag) # Running an OLS for every group.
  PACF[j] = coef(fit)[j - 1]           # Getting the slope for the last lagged ts.
}

Y, por último, trazando de nuevo lado a lado, los cálculos generados por R y los manuales:

Que la idea es correcta, aparte de probables problemas computacionales, puede verse comparando PACF a pacf(st.y, plot = F) .

---

código aquí .

Answer 3

Pues bien, en la práctica encontramos error (ruido) que se representa por $ e_t $ las bandas de confianza ayudan a averiguar si un nivel puede considerarse sólo ruido (porque aproximadamente el 95% de los tiempos estarán dentro de las bandas).

Answer 4

He aquí un código python para calcular la ACF:

def shift(x,b):
    if ( b <= 0 ):
        return x
    d = np.array(x);
    d1 = d
    d1[b:] = d[:-b]
    d1[0:b] = 0
    return d1

# One way of doing it using bare bones
# - you divide by first to normalize - because corr(x,x) = 1
x = np.arange(0,10)
xo = x - x.mean()

cors = [ np.correlate(xo,shift(xo,i))[0]  for i in range(len(x1)) ]
print (cors/cors[0] )

#-- Here is another way - you divide by first to normalize
cors = np.correlate(xo,xo,'full')[n-1:]
cors/cors[0]