Lección 4. ACF, PACF y densidad espectral de modelos MA

# Módulos necesarios
import statsmodels.api as sm
from statsmodels.tsa.arima_process import ArmaProcess
from scipy import signal
from scipy.signal import freqz
from numpy.polynomial import Polynomial
import math

La siguiente función que calcula la PACF a partir de la ACF resolviendo las ecuaciones de Yule-Walker recursivamente con el algoritmo Levinson-Durbin; véase Pourahmadi, M. (2001, Capítulo 7) o Brockwell & Davis (1991, Capitulo 5)).

def theoretical_pacf_from_acf(acf_vals):
    """
    Calcula la PACF teórica (para lags >= 0) de forma paramétrica usando el algoritmo de Durbin-Levinson.
    Se asume que acf_vals es un array tal que acf_vals[0] corresponde a lag 0.
    
    Retorna un array donde el coeficiente en la posición k corresponde a la PACF en el lag k.
    """
    nlags = len(acf_vals)-1
    # Matriz para almacenar los coeficientes (phi)
    phi = np.zeros((nlags+1, nlags+1))
    # v guarda la varianza de error en cada paso
    v = np.zeros(nlags+1)
    pacf = np.zeros(nlags+1)
    v[0] = acf_vals[0]
    pacf[0] = 1.0  # Por convención, PACF(0)=1.
    
    for k in range(1, nlags+1):
        if k == 1:
            phi[1,1] = acf_vals[1] / acf_vals[0]
            pacf[1] = phi[1,1]
            v[1] = v[0]*(1 - phi[1,1]**2)
        else:
            sum_val = 0.0
            # Sumar phi[k-1,j]*acf_vals[k-j] para j=1,...,k-1.
            for j in range(1, k):
                sum_val += phi[k-1, j] * acf_vals[k - j]
            phi[k, k] = (acf_vals[k] - sum_val) / v[k-1]
            pacf[k] = phi[k, k]
            for j in range(1, k):
                phi[k, j] = phi[k-1, j] - phi[k, k] * phi[k-1, k - j]
            v[k] = v[k-1] * (1 - phi[k, k]**2)
    return pacf  # Incluye PACF en lag 0 con valor 1.

La siguiente función pinta las ACF, PACF y densidad espectral teóricas.

def plot_arma_parametric_diagnostics(ar_params=[1,], ma_params=[1,], sigma2=1, lags=20):
    """
    Genera y retorna una figura con tres subgráficas para un modelo ARMA:
      - ACF teórica (omitiendo retardo 0)
      - PACF teórica paramétrica (calculada vía el algoritmo de Durbin-Levinson, omitiendo retardo 0)
      - Densidad espectral teórica
    
    Parámetros:
      ar_params : array-like
          Coeficientes del polinomio AR (incluye el 1 inicial).
      ma_params : array-like
          Coeficientes del polinomio MA (incluye el 1 inicial).
      sigma2 : float
          Varianza del ruido Uₜ.
      lags : int, opcional
          Número de retardos considerados (por defecto 20).
    
    Retorna:
      fig : objeto matplotlib.figure.Figure
    """
    phi   = Polynomial(ar_params, symbol='B')
    theta = Polynomial(ma_params, symbol='B')
    # Crear el proceso ARMA
    arma_process = ArmaProcess(ar=phi.coef, ma=theta.coef)
    
    # 1. Calcular la ACF teórica
    acf_theo = arma_process.acf(lags=lags)
    lags_theo = np.arange(len(acf_theo))
    # Omitir el retardo 0 para gráficos de ACF y PACF:
    acf_plot = acf_theo[1:]
    lags_plot = lags_theo[1:]
    
    # 2. Calcular la PACF teórica de forma paramétrica a partir de la ACF
    pacf_theo = theoretical_pacf_from_acf(acf_theo)
    pacf_plot = pacf_theo[1:]
    lags_pacf = np.arange(len(pacf_theo))[1:]
    
    # 3. Calcular la densidad espectral teórica usando freqz
    w, h = freqz(theta.coef, phi.coef, worN=1024)
    freq = w / (2 * np.pi)  # Conversión de radianes/muestra a ciclos/muestra
    spectrum = (sigma2 / (2 * np.pi)) * np.abs(h)**2
    
    # 4. Graficar los tres subgráficos
    fig, axs = plt.subplots(1, 3, figsize=(18, 3))
    
    # ACF teórica
    axs[0].stem(lags_plot, acf_plot, basefmt=" ")
    axs[0].axhline(0, color='black', lw=0.5, linestyle='--')
    axs[0].set_xlabel('Retardo')
    axs[0].set_ylabel('ACF')
    
    # PACF teórica paramétrica
    axs[1].stem(lags_pacf, pacf_plot, basefmt=" ")
    axs[1].axhline(0, color='black', lw=0.5, linestyle='--')
    axs[1].set_xlabel('Retardo')
    axs[1].set_ylabel('PACF')
    
    # Densidad espectral teórica
    axs[2].plot(freq, spectrum)
    axs[2].axhline(0, color='black', lw=0.5, linestyle='--')
    axs[2].set_xlabel('Frecuencia')
    axs[2].set_ylabel('Densidad Espectral')
    
    plt.tight_layout()
    plt.close(fig)
    return fig

Breve explicación de la anterior función plot_arma_parametric_diagnostics:

1. La función recibe los parámetros del modelo ARMA y el número de retardos.
2. Se calcula la ACF teórica y a partir de ella la PACF con theoretical_pacf_from_acf; en la gráfica se omite el retardo cero.

3. La densidad espectral teórica \( S(\omega) \) se obtiene mediante la función freqz y la fórmula \[S(\omega)=\frac{\sigma^2}{2\pi}\Biggl|\frac{\sum_{j=0}^{q}\theta_j e^{-i2\pi \omega\,j}}{\sum_{k=0}^{p}\phi_k e^{-i2\pi \omega\,k}}\Biggr|^2.\] Donde:

   • \(\sigma^2\) es la varianza del término de error (innovación) del modelo.
   • \((\theta_j \mid j=0:q)\) son los coeficientes del polinomio de medias móviles (MA), de modo que el polinomio MA es: \(\theta(z) = \sum_{j=0}^{q}\theta_j\,z^j\).
   • \((\phi_k \mid k=0:p)\) son los coeficientes del polinomio autorregresivo (AR), de modo que el polinomio AR es: \(\phi(z) = \sum_{k=0}^{p}\phi_k\,z^k\).
   • Al evaluar estos polinomios en \(z = e^{-i\,2\pi \omega}\) se obtiene la respuesta en frecuencia del modelo: \[H(e^{-i\,2\pi \omega})=\frac{\theta\bigl(e^{-i\,2\pi \omega}\bigr)}{\phi\bigl(e^{-i\,2\pi \omega}\bigr)}.\]
   • El módulo al cuadrado de \( H(e^{-i\,2\pi \omega}) \) multiplicado por \(\frac{\sigma^2}{2\pi}\) es la densidad espectral.

   Esta fórmula es fundamental porque permite conocer la distribución de la varianza del proceso a lo largo de las diferentes frecuencias, lo cual es clave para analizar las propiedades dinámicas en modelos ARMA.

4. Se crean tres subgráficos, y se añade una línea horizontal en y=0 en cada uno (axhline).
5. Por último, la función retorna el objeto figura.

Función que simula un modelo ARIMA y pinta la serie temporalm, las ACF, PACF y densidad espectral estimadas con los datos simulados.

def plot_arma_analysis(ar_params=[1,], ma_params=[1,], sigma2=1, lags=20, n=400, seed=0):

    lags=math.floor(n/2) if n/2 < lags else lags
    
    phi   = Polynomial(ar_params, symbol='B')
    theta = Polynomial(ma_params, symbol='B')
    arma_process = sm.tsa.ArmaProcess(phi.coef, theta.coef)
    
    # Simulación de la serie temporal
    if seed:
        np.random.seed(seed)
    data = arma_process.generate_sample(nsample=n)
       
    # Crear la figura
    fig = plt.figure(figsize=(18, 7))
    
    # Subgráfico de la serie temporal
    ax1 = fig.add_subplot(2, 1, 1)
    ax1.plot(data, color='blue')
    ax1.set_xlabel('t')
    
    # Crear una fila de 3 subgráficos para ACF, PACF y periodograma
    ax2 = fig.add_subplot(2, 3, 4)  # Fila 2, Columna 1
    sm.graphics.tsa.plot_acf(data, lags=lags, ax=ax2, title='ACF Estimada', zero=False, auto_ylims=True)
    ax2.set_xlabel('Retardo')

    ax3 = fig.add_subplot(2, 3, 5)  # Fila 2, Columna 2
    sm.graphics.tsa.plot_pacf(data, lags=lags, ax=ax3, title='PACF Estimada', zero=False, auto_ylims=True)
    ax3.set_xlabel('Retardo')
    
    ax4 = fig.add_subplot(2, 3, 6)  # Fila 2, Columna 3
    f, Pxx = signal.welch(data, fs=1, nperseg=256)
    ax4.plot(f, Pxx, color='purple')
    ax4.set_title('Periodograma Estimado')
    ax4.set_xlabel('Frecuencia')
    #ax4.set_ylabel('Densidad Espectral')
    ax4.set_xlim(0, 0.5)  # Limitar el eje x a 0.5

    plt.tight_layout()
    plt.close(fig)    
    return fig

Si un proceso estocástico es estacionario, la covarianza entre \(X_t\) y \(X_{t+k}\) no depende de \(t\); tan solo depende de la distancia temporal \(k\) entre ambas variables aleatorias: \(\;Cov(X_t,X_{t-k})=\gamma_k.\;\)

• La secuencia \(\boldsymbol{\gamma}=(\gamma_k\mid k\in\mathbb{Z})\) se denomina función de autocovarianzas

• Dividiendo \(\boldsymbol{\gamma}\) por \(\gamma_0\) obtenemos la secuencia \(\boldsymbol{\rho}=(\rho_k\mid k\in\mathbb{Z})\) donde \[\rho_k=\frac{Cov(X_t,X_{t-k})}{\sqrt{Var(X_t)Var(X_{t-k})}}=\frac{\gamma_k}{\gamma_0},\]

  que se denomina función de autocorrelación (ACF).

Veamos otras transformaciones de \(\boldsymbol{\gamma}\) que subrayan diferentes características de la interdependencia temporal en un proceso estocástico.

La correlación parcial entre dos variables \(X_t\) y \(X_{t-k}\) de un proceso estacionario \(\boldsymbol{X}\) mide su correlación una vez descontado el efecto de las variables \(X_{t-1},\ldots X_{t-(k-1)}\) que median entre ambas. \[X_t,\ \overbrace{X_{t-1},\ldots X_{t-(k-1)},}\ X_{t-k}\]

Si denotamos con \(\widehat{X_t}\) y \(\widetilde{X_{t-k}}\) los ajustes de las respectivas regresiones de \(X_t\) y \(X_{t-k}\) sobre \(X_{t-1},\ldots X_{t-(k-1)};\;\) la correlación parcial \(\pi_k\) entre \(X_t\) y \(X_{t+k}\) es la correlación entre los residuos de sendas regresiones: \[\pi_k=Corr\Big((X_t-\widehat{X_t}),\ (X_{t-k}-\widetilde{X_{t-k}})\Big)\]

En un proceso débilmente estacionario las correlaciones parciales solo dependen de la distancia \(k\), lo que permite definir la siguiente secuencia.

La función de autocorrelación parcial (PACF) de un proceso estocástico estacionario \(\boldsymbol{X}\) es la siguiente secuencia \(\boldsymbol{\pi}\) simétrica (\(\pi_{-k}=\pi_k\)):

\[\boldsymbol{\pi}=(\pi_k\mid k\in\mathbb{Z})= \begin{cases} \pi_0 & = \rho_0 ={\color{blue}{1}}\\\\ \pi_1 &=Corr(X_t,\ X_{t-1})=\rho_1=\frac{\gamma_1}{\gamma_0}\\\\ \pi_k &=Corr\Big((X_t-\widehat{X_t}),\ (X_{t-k}-\widetilde{X_{t-k}})\Big) \end{cases},\] donde \(\widehat{X_t}\) y \(\widetilde{X_{t-k}}\) son los ajustes de las respectivas regresiones de \(X_t\) y \(X_{t-k}\) sobre \(X_{t-1},\ldots X_{t-(k-1)};\;\) i.e., las proyecciones ortogonales sobre \(\bar{sp}(1,X_{t-1},\ldots X_{t-(k-1)})\).

La magnitud de la correlación parcial \(\;\pi_k\;\) refleja la mejora en la predicción de \(\widehat{X_t}\) si en lugar de usar una combinación lineal con solo los \(k-1\) primeros retardos, \(\widehat{X_t}=\widehat{\alpha_0}+\sum_{j=1}^{k-1}\widehat{\alpha_j} X_{t-j}\), empleamos \(k\) retardos (i.e., un retardo más).

Así, la PACF nos ayudará a identificar el orden de procesos autoregresivos.

Hay una correspondencia uno-a-uno entre la función de autocovarianzas \(\boldsymbol{\gamma}\) y la PACF \(\boldsymbol{\pi}\).

Es decir, es posible reconstruir una de las secuencias a partir de la otra (por ejemplo resolviendo las ecuaciones de Yule-Walker recursivamente con el algoritmo Levinson-Durbin; véase Pourahmadi, M. (2001, Capítulo 7) o Brockwell & Davis (1991, Capitulo 5))

Así, la PACF \(\boldsymbol{\pi}\) puede verse como una reparametrización de la función de autocovarianzas \(\boldsymbol{\gamma}\).

Veamos otra transformación de \(\boldsymbol{\gamma}\) que también arrojará luz sobre las propiedades de un proceso estocástico.

Si \(\boldsymbol{X}\) es un proceso estocástico con función de autocovarianzas \(\boldsymbol{\gamma}\in\ell^1\), es decir, con función de autocovarianzas absolutamente sumable, \[\sum_{k\in\mathbb{Z}}|\gamma_k|<\infty,\] entonces definimos la densidad espectral de \(\boldsymbol{X}\) como \[f(\omega) \quad=\quad\frac{1}{2\pi}\sum_{k\in\mathbb{Z}}\gamma_k e^{-i k\omega} \quad=\quad{\color{blue}{\frac{1}{2\pi}\sum_{k\geq0}\gamma_k \cos(k\omega)}},\qquad \omega\in[-\pi,\pi],\] (donde \(w\) son las frecuencias de oscilación de la serie).

La densidad espectral \(f(\omega)\) satisface las siguientes propiedades:

1. \(f(\omega)=f(-\omega)\)
2. \(f(\omega)\geq0\)
3. \(\int_{-\pi}^\pi f(\omega)d\omega<\infty\)

Es más, partiendo de la densidad espectral se pueden calcular las covarianzas \[\gamma_k=\int_{-\pi}^\pi f(\omega)e^{-i k\omega}d\omega;\quad k\in\mathbb{Z}\]

Consecuentemente, para \(k=0\) \[\sigma^2=\gamma_0=\int_{-\pi}^\pi f(\omega)d\omega.\]

Por tanto podemos interpretar la densidad espectral como una descomposición de la varianza en sumas de oscilaciones con distintas frecuencias (o periodos).

(\(Periodo=\frac{1}{Frecuencia}\), donde la Frecuencia = Nº de ciclos por unidad de tiempo).

La ACF, la PACF y la densidad espectral son funciones que dependen únicamente de los dos primeros momentos de la distribución.

Su estimación y posterior análisis son la herramienta fundamental para elegir un modelo ARMA para una serie temporal ``estacionaria''.

Para entenderlo debemos ver cómo son estas funciones en algunos modelos lineales concretos.

Siempre estacionario.

Para ser invertible raíces \(\boldsymbol{\theta}\) fuera del círculo unidad \((\boldsymbol{\theta}^{-\triangleright}=\boldsymbol{\theta}^{-1}\in\ell^1)\)

Tipos de representación del proceso

Como suma ponderada finita

(número finito de parámetros) \[\boldsymbol{X}=\boldsymbol{\theta}*\boldsymbol{U} \quad\Rightarrow\quad X_t=\boldsymbol{\theta}(\mathsf{B})U_t\]

Como suma ponderada infinita

(solo existe si es invertible): \[\frac{1}{\boldsymbol{\theta}}*\boldsymbol{X}=\boldsymbol{U} \quad\Rightarrow\quad \frac{1}{\boldsymbol{\theta}}(\mathsf{B})X_t=U_t\]

ACF (\(\boldsymbol{\rho}\)): Grado (\(q\)) y cogrado (\(-q\))

PACF (\(\boldsymbol{\pi}\)): Grado (\(\infty\)) y cogrado (\(-\infty\)): exponenciales y/o sinusoidales amortiguadas.

Densidad espectral Suma de dos cosenos más una constante.