Prueba de independencia de los residuos

Pruebas de Autocorrelación Individual (por Lag)

Función que arroja en una tabla la prueba de independencia (autocorrelación y no autocorrelación) de los residuos para cada lag de la afc en conjunto con el test de Ljung-Box o Box -Pierce.

  1. LjBox.Pierce_Inde(model, max.lag = 25, type = c("Ljung-Box", "Box-Pierce"))

  • model : Se ingresa el modelo Arima() generado con library(forecast).

  • max.lag : Cantidad de lag o retrasos que se desean considerar en la tabla, por defecto arroja 25.

  • type : Tipo de test que deses que evalúen a los lag de afc de los residuos, se puede elegir entre "Ljung-Box", "Box-Pierce"

uso los datos AirPassengers para probar las funciones.

modelo_arima<-forecast::auto.arima(AirPassengers)
library(forecast)
library(tseries)
library(ggplot2)
library(knitr)


LjBox.Pierce_Inde <- function(model, max.lag = 25, type = c("Ljung-Box", "Box-Pierce")) {
  type <- match.arg(type)
  
  # Aqui extraigo los residuos
  residuos <- resid(model)
  
  #  ACF y PACF
  acf_resid <- acf(residuos, type = "correlation", plot = FALSE, lag.max = max.lag)
  pacf_resid <- pacf(residuos, plot = FALSE, lag.max = max.lag)
  
  # Genero las tabla de ACF y PACF
  Lag <- 1:max.lag
  acf_table <- data.frame(ACF = round(acf_resid$acf[-1], 3))  # Excluyo el lag 0 y redonde
  pacf_table <- data.frame(PACF = round(pacf_resid$acf, 3))  # Redondeo
  
  # Prueba de Ljung-Box o Box-Pierce
  p.value <- rep(NA, max.lag)
  Qm <- rep(NA, max.lag)
  fit <- sum(arimaorder(model)[c("p", "q", "P", "Q")], na.rm = TRUE)
  
  for (i in 1:max.lag) {
    if (i > fit) {
      test <- Box.test(residuos, type = type, fitdf = fit, lag = i)
      p.value[i] <- test$p.value
      Qm[i] <- test$statistic
    } else {
      # Relleno a Qm con el valor estadístico para los lags <= fit
      test <- Box.test(residuos, type = type, fitdf = 0, lag = i)
      Qm[i] <- test$statistic
    }
  }
  
  # Redondeo a 3 decimales
  Qm <- round(Qm, 3)
  p.value <- round(p.value, 3)
  
  # data frame final
  autocorres <- data.frame(Lag, ACF = round(acf_table$ACF, 3), PACF = pacf_table$PACF, Qm, pvalue = p.value)
  # Convertir NA a cadena vacía debido alos terminos de cada modelo
  autocorres_print <- autocorres
  autocorres_print[is.na(autocorres_print)] <- " "
  
  return(autocorres)
}

Probando la función

# pruebo la funcióm 
LJ <- LjBox.Pierce_Inde(modelo_arima, max.lag = 30, type = "Ljung-Box")

# para tabla en consola
LJ_print <- LJ
#LJ_print$pvalue <- ifelse(is.na(LJ_print$pvalue), " ", LJ_print$pvalue)
#print(LJ_print, row.names = FALSE)
library(kableExtra)
# para gnrar una tabla en formato HTML
LJ_print %>%
  kable(format = "html", table.attr = "class='table table-striped'", align = "c") %>%
  kable_styling(full_width = FALSE, position = "center")
Lag ACF PACF Qm pvalue
1 -0.001 -0.001 0.000 NA
2 0.038 0.038 0.211 NA
3 -0.077 -0.077 1.084 NA
4 -0.096 -0.099 2.480 0.115
5 0.054 0.061 2.924 0.232
6 0.012 0.014 2.946 0.400
7 -0.070 -0.092 3.703 0.448
8 -0.032 -0.034 3.858 0.570
9 0.149 0.175 7.294 0.295
10 -0.094 -0.112 8.688 0.276
11 -0.004 -0.048 8.691 0.369
12 -0.120 -0.080 10.973 0.278
13 0.049 0.086 11.358 0.330
14 0.028 -0.019 11.482 0.404
15 0.025 -0.004 11.582 0.480
16 -0.126 -0.116 14.186 0.361
17 -0.042 -0.013 14.479 0.415
18 -0.005 -0.031 14.483 0.489
19 -0.165 -0.186 19.052 0.266
20 -0.143 -0.199 22.524 0.165
21 -0.038 0.027 22.773 0.199
22 -0.107 -0.176 24.744 0.169
23 0.250 0.192 35.620 0.017
24 0.111 0.085 37.784 0.014
25 -0.025 0.015 37.891 0.019
26 0.042 -0.006 38.213 0.024
27 -0.100 -0.064 40.000 0.021
28 -0.041 -0.070 40.311 0.027
29 -0.023 -0.040 40.405 0.036
30 -0.080 -0.132 41.582 0.036

Bibliotecas y Funciones Utilizadas:

  • forecast: Utilizas la función resid() para extraer los residuos de los modelos ARIMA o SARIMA.

  • tseries: Utilizas Box.test() para realizar las pruebas de Ljung-Box y Box-Pierce.

  • stats (cargada implícitamente con R): Utilizas acf() y pacf() para calcular la función de autocorrelación y la autocorrelación parcial de los residuos.

Visualización de diagnostico de los residuos

Función que arroja una matriz de gráficos residuales conteniendo los residuos como tal, el correlograma afc, un histograma, y el gráfico de la evaluación de “Ljung-Box” o “Box-Pierce”(según se elija), para cada lag (afc) de los residuos.

check_residuals(model, lag.max = 25, title = "nombre de modelo")

# aqui incluyo los graficos que forman parte de la funcion de autocorrelación para las pruebas de Ljung-Box y Box-Pierce, en conjunto con las autocorrelciones. Ademas de la grafica de residuos y su distribución

library(ggplot2)
library(gridExtra)
library(forecast)

check_residuals <- function(model, lag.max = 25, title = "Modelo SARIMA", type = "Ljung-Box") {
  
  # Extraer residuos
  residuos <- resid(model)
  
  # Gráfico de los residuos
  P5 <- autoplot(residuos, series = "Residuos", color = "#8B8B00") +
    geom_point(aes(y = residuos), size = 1, color = "#09622A") + 
    labs(y = "Residuos", x = "Índice") +
    ggtitle(paste("Residuos del", title)) +
    theme(
      plot.title = element_text(hjust = 0.5),
      axis.text.x = element_text(angle = 0, hjust = 1)
    )
  
  # Función para graficar ACF de residuos
  afc.re <- function(data, lag.max) {
    acf_result <- acf(data, lag.max = lag.max, plot = FALSE)
    significant <- abs(acf_result$acf) > 1.96/sqrt(length(data))
    
    p <- ggplot() +
      geom_blank() +
      xlim(0, lag.max) +
      labs(x = "Lag", y = "ACF") +
      geom_rect(aes(xmin = 0, xmax = lag.max, ymin = -1.96/sqrt(length(data)), ymax = 1.96/sqrt(length(data))), 
                fill = rgb(139/255, 69/255, 19/255, alpha = 0.2), color = NA) +
      geom_hline(yintercept = c(-1.96/sqrt(length(data)), 0, 1.96/sqrt(length(data))), 
                 linetype = "dashed", color = "brown") +
      geom_hline(yintercept = 0, linetype = "solid", color = "black") +
      geom_segment(aes(x = 1:lag.max, y = 0, xend = 1:lag.max, yend = acf_result$acf[2:(lag.max + 1)]), 
                   color = ifelse(significant[2:(lag.max + 1)], "black", "#005F52"), size = 1, linetype = "solid") +
      geom_point(aes(x = 1:lag.max, y = acf_result$acf[2:(lag.max + 1)]), 
                 shape = 21, size = 3, fill = ifelse(significant[2:(lag.max + 1)], "black", "#005F52")) +
      expand_limits(y = c(-1, 1)) +  
      coord_cartesian(ylim = c(-0.25, 0.25))
    
    return(p)
  }
  
  P6 <- afc.re(residuos, lag.max = lag.max)
  
  # Histograma de los residuos
  histores <- as.vector(residuos)
  residh <- as.data.frame(histores)
  
  P8 <- ggplot(residh, aes(x = histores)) +
    geom_histogram(bins = 60, aes(y = after_stat(density)),
                   colour = 1, fill = "#BD5DA6") +
    labs(y = 'Frecuencia', x = "Residuos del modelo") +
    geom_rug() +
    geom_density(lwd = 0.3, colour = 4, fill = 4, alpha = 0.29) +
    stat_function(fun = dnorm, args = list(mean = mean(residh$histores), sd = sd(residh$histores)), 
                  geom = "area", fill = "black", alpha = 0.2, linetype = "dashed")
  
  # Función Ljung-Box
  LjBox.Pierce_Inde <- function(model, max.lag = 25, type = type) {
    type <- match.arg(type)
    residuos <- resid(model)
    p.value <- rep(NA, max.lag)
    Qm <- rep(NA, max.lag)
    fit <- sum(arimaorder(model)[c("p", "q", "P", "Q")], na.rm = TRUE)
    
    for (i in 1:max.lag) {
      if (i > fit) {
        test <- Box.test(residuos, type = type, fitdf = fit, lag = i)
        p.value[i] <- test$p.value
        Qm[i] <- test$statistic
      } else {
        test <- Box.test(residuos, type = type, fitdf = 0, lag = i)
        Qm[i] <- test$statistic
      }
    }
    
    Qm <- round(Qm, 3)
    p.value <- round(p.value, 3)
    autocorres <- data.frame(Lag = 1:max.lag, Qm, pvalue = p.value)
    
    return(autocorres)
  }
  
  LJ <- LjBox.Pierce_Inde(model, max.lag = lag.max, type = type)
  LJ_valid <- LJ[!is.na(LJ$pvalue), ]
  
  P9 <- ggplot(LJ_valid, aes(x = Lag, y = pvalue)) +
    geom_point(stat = "identity", size = 2, shape = 16, color = "#4CAF50", fill = "#4CAF50") +
    geom_hline(yintercept = 0.05, linetype = "dashed", color = "red", size = 1) +
    geom_hline(yintercept = 0.01, linetype = "dashed", color = "blue", size = 0.5) +  
    labs(title = "P_values for Ljung-Box statistic", x = "Lag (H)", y = "p-value") +
    scale_x_continuous(breaks = 1:lag.max) +
    theme(plot.title = element_text(hjust = 0.5, size = 10))
  
  
  layout <- rbind(
    c(5, 5),
    c(6, 8),
    c(9, 9)
  )
  
  grid.arrange(
    P5, P6, P8, P9,
    layout_matrix = layout
  )
}

Probando la función:

check_residuals(modelo_arima, lag.max = 30, title = "ARIMA(2,1,1)(0,1,0)[12] ")

Bibliotecas y Funciones Utilizadas:

Test de los 3 supuestos de los residuos

Función que crea una tabla conteniendo test de los 3 supuestos que deben cumplir los residuos de modelos de la familia ARIMA, en la primera columna imprime el tipo de evaluación, en la segunda el tipo de test, en la tercera el estadístico de prueba, en la cuarta el valor de significancia y finalmente la decisión en base al p-valor.

  1. test.residuals(modelo_arima, test_normal = c("Jarque-Bera", "Shapiro-Wilk", "Kolmogorov-Smirnov")

    • modelo_arima : Se introduce el modelo

    • test_normal:Se puede elegir entre Jarque-Bera, Shapiro-Wilk y Kolmogorov-Smirnov para probar la normalidad de los residuos por defecto arroja el gráfico con Shapiro-Wilk.

library(lmtest)
library(kableExtra)

test.residuals <- function(modelo_arima, test_normal = "Shapiro-Wilk", nrepl = 2000) {
  residuos <- residuals(modelo_arima)
  
  # Seleccionar prueba de normalidad
  if (test_normal == "Shapiro-Wilk") {
    test_normalidad <- shapiro.test(residuos)
    test_name <- "Shapiro-Wilk"
  } else if (test_normal == "Jarque-Bera") {
    # Prueba de Jarque-Bera personalizada
    NAME <- deparse(substitute(residuos))
    n <- length(residuos)
    x_mean <- mean(residuos)
    skewness <- sqrt(n) * sum((residuos - x_mean)^3) / (sum((residuos - x_mean)^2)^(3/2))
    kurtosis <- n * sum((residuos - x_mean)^4) / (sum((residuos - x_mean)^2)^2)
    jb_statistic <- (n / 6) * (skewness^2 + ((kurtosis - 3)^2) / 4)
    
    count <- 0
    for (i in 1:nrepl) {
      z <- rnorm(n)
      z_mean <- mean(z)
      a1 <- sqrt(n) * sum((z - z_mean)^3) / (sum((z - z_mean)^2)^(3/2))
      a2 <- n * sum((z - z_mean)^4) / (sum((z - z_mean)^2)^2)
      T <- (n / 6) * (a1^2 + ((a2 - 3)^2) / 4)
      if (T > jb_statistic) {
        count <- count + 1
      }
    }
    
    p_value <- count / nrepl
    test_normalidad <- list(statistic = jb_statistic, p.value = p_value, method = "Jarque-Bera", data.name = NAME)
    class(test_normalidad) <- "htest"
    test_name <- "Jarque-Bera"
  } else if (test_normal == "Kolmogorov-Smirnov") {
    test_normalidad <- ks.test(residuos, "pnorm", mean(residuos), sd(residuos))
    test_name <- "Kolmogorov-Smirnov"
  }
  
  # Test de Ljung-Box para autocorrelación
  test_acf <- Box.test(residuos, lag = 20, type = "Ljung-Box")
  
  # Test de Breusch-Pagan para heterocedasticidad
  fitted_values <- fitted(modelo_arima)
  modelo_bp <- lm(I(residuos^2) ~ fitted_values)
  test_bp <- bptest(modelo_bp)
  
  # Compilado de resultados
  resultados <- data.frame(
    Evaluacion = c("Normalidad", "Autocorrelación", "Heterocedasticidad"),
    Test = c(test_name, "Ljung-Box", "Breusch-Pagan"),
    Estadistico = c(as.numeric(test_normalidad$statistic), as.numeric(test_acf$statistic), as.numeric(test_bp$statistic)),
    P_valor = c(test_normalidad$p.value, test_acf$p.value, test_bp$p.value),
    Decisión = c(ifelse(test_normalidad$p.value < 0.05, "No Normal", "Normal"),
                 ifelse(test_acf$p.value < 0.05, "Autocorrelación", "Sin Autocorrelación"),
                 ifelse(test_bp$p.value < 0.05, "Heterocedasticidad", "Homocedasticidad"))
  )
  
  return(resultados)
}

Pruebo la función

Modeloair <- test.residuals(modelo_arima, test_normal = "Jarque-Bera")
Modeloair %>%
  kable(format = "html", table.attr = "class='table table-striped'", align = "c") %>%
  kable_styling(full_width = FALSE, position = "center")
Evaluacion Test Estadistico P_valor Decisión
Normalidad Jarque-Bera 15.130694 0.0065000 No Normal
Autocorrelación Ljung-Box 22.523605 0.3127870 Sin Autocorrelación
Heterocedasticidad Breusch-Pagan 4.699511 0.0301712 Heterocedasticidad

Bibliotecas y Funciones Utilizadas: