Validación cruzada en R

La siguiente entrada muestra una alternativa para obtener una serie histórica de errores de predicción de una serie temporal (validación cruzada) a partir de modelos ARIMA con represores (xARIMA). La entrada se segmenta en dos partes, primero muestra un ejemplo aplicado del uso de la función y segundo explica brevemente la construcción de la función de validación cruzada y de sus funciones implícitas.

a.       Ejemplo de uso de la función

El ejemplo utiliza datos del crecimiento trimestral del PIB real de la República Dominicana, publicados por el Banco Central de la República Dominicana, para el periodo 1991q1-2020q1.

library(forecast)

library(readxl)

library(tidyverse)

library(zoo)


setwd("C:/Users/Nery Ramirez/Dropbox/densidadCrecimiento/densidadpib")

data_prib <- read_excel("Data_pib_real.xlsx")

 

> data_prib

# A tibble: 117 x 2

   t      pib_cons_real_2007

   <chr>               <dbl>

 1 1991Q1               20.5

 2 1991Q2               23.9

 3 1991Q3               22.5

 4 1991Q4               24.4

 5 1992Q1               26.7

 6 1992Q2               26.6

 7 1992Q3               29.1

 8 1992Q4               26.8

 9 1993Q1               30.0

10 1993Q2               32.6

# ... with 107 more rows

Utilizando el paquete dplyr del tidyverse y la función as.yearqtr de la librería zoo, se recuperan las fechas correspondientes a cada trimestre, se calcula el crecimiento anual (utilizando cuatro rezagos en la función lag) y posteriormente se crea una variable binaria llamada d_cre_b1, la misma asume el valor de 1 en los periodos de crecimiento negativo posterior al 2002, y en los 10 trimestres donde el crecimiento fue menor antes de este año (esto porque no se registraron periodos de caídas interanual en la producción trimestral).

data_prib1 <- data_prib %>%

  mutate(q = as.yearqtr(time(pib_cons_real_2007), format = "q"),

         anio = as.numeric(substring(t, 1,4)),

         creci = (pib_cons_real_2007/dplyr::lag(pib_cons_real_2007,4)) - 1,

         d_cre_b1 = (creci %in% sort(creci[anio<=2002])[1:10] | creci<0)*1,

  ) %>%

  # Reordenar el nombre de las variables

  select(t,q,everything())

 

data_prib1

# A tibble: 117 x 6

   t      q         pib_cons_real_2007  anio  creci d_cre_b1

   <chr>  <yearqtr>              <dbl> <dbl>  <dbl>    <dbl>

 1 1991Q1 1 Q1                    20.5  1991 NA           NA

 2 1991Q2 2 Q1                    23.9  1991 NA           NA

 3 1991Q3 3 Q1                    22.5  1991 NA           NA

 4 1991Q4 4 Q1                    24.4  1991 NA           NA

 5 1992Q1 5 Q1                    26.7  1992  0.299        0

 6 1992Q2 6 Q1                    26.6  1992  0.114        0

 7 1992Q3 7 Q1                    29.1  1992  0.294        0

 8 1992Q4 8 Q1                    26.8  1992  0.101        0

 9 1993Q1 9 Q1                    30.0  1993  0.124        0

10 1993Q2 10 Q1                   32.6  1993  0.226        0

# ... with 107 more rows

Ahora, usamos la función (se explica debajo) error_historico_d1 para obtener una serie histórica de los errores de pronósticos a diversos pasos. La serie x es la serie de crecimiento anteriormente estimada, pero omitiendo los na, mientras que el regresor x del modelo corresponde a la dummy para indicar periodos de crecimiento económico “bajo” o negativos. Los argumentos de la función son (note que solo los dos primeros argumentos son obligatorios):

-          x: serie de tiempo sobre la que se desea trabajar.

-          xreg_d: regresores del modelo ARIMA.

-          n.valid: número de periodos que se desean proyectar.

-          Recursiva: indicar que queremos que la ventana de estimación sea móvil o recursiva.

-          Ventana: tamaño de la ventana, el valor predeterminado es 70% de las observaciones.

-          Medida_error: indicamos la serie de errores que deseamos guardar. Las opciones disponibles se colocan debajo.

("ME", "MSE", "RMSE", "MAE", "MPE", "MAPE","ACF1", "Theil's U")

-          Agregar: se muestra un agregado de las series de errores, puede ser un promedio, mediana o la función de interés, siempre que esté disponible en R.

x<-data_prib1$creci[-c(1:4)]

xreg_d<-data_prib1$d_cre_b1[-c(1:4)]

 

error_historico_d1(x,

                   xreg_d,

                   n.valid=8,

                   recursiva=TRUE,

                   ventana=round(length(!is.na(x))*0.7),

                   medida_error="RMSE",

                   agregar = "mean")

El modelo arroja por salida una lista con información sobre la estimación, el ultimo modelo autoarima seleccionado por el algoritmo desarrollado en el paquete forecast, la serie histórica de errores a cada paso, y los errores agregados, que para el caso actual representa el promedio de errores.

$Informacion

     [,1]                  

[1,] "Medida de errro: RMSE"

[2,] "Horizonte: 8"        

[3,] "Función agregar: mean"

 

$`Modelo Arima`

Series: x.train

Regression with ARIMA(1,0,0) errors

 

Coefficients:

         ar1  intercept     xreg

      0.3090     0.1407  -0.1949

s.e.  0.1059     0.0120   0.0209

 

sigma^2 estimated as 0.005801:  log likelihood=122.84

AIC=-237.68   AICc=-237.28   BIC=-227.06

 

$`Hitorico errores`

                      h1         h2         h3         h4         h5         h6         h7         h8

error_pronH 0.0120526909 0.01155399 0.01390647 0.01437705 0.01369995 0.03722491 0.04796603 0.04526843

error_pronH 0.0158480814 0.01798974 0.01688961 0.01571486 0.04056176 0.05165188 0.04824761 0.04778444

error_pronH 0.0138126750 0.01552025 0.01395166 0.04467589 0.05622090 0.05179590 0.05078809 0.04789154

error_pronH 0.0224250648 0.01677765 0.05223426 0.06335995 0.05719506 0.05523534 0.05154089 0.04958089

error_pronH 0.0164286782 0.06017829 0.07053814 0.06166394 0.05861779 0.05396812 0.05146070 0.05575162

error_pronH 0.0898086733 0.08977179 0.07400109 0.06776404 0.06108429 0.05734009 0.06103166 0.05711851

error_pronH 0.0597617913 0.04254667 0.04443076 0.03958809 0.03817338 0.04751249 0.04402317 0.04143310

error_pronH 0.0164400564 0.04262606 0.03739458 0.03541055 0.04663292 0.04258040 0.03962853 0.04393133

error_pronH 0.0517267059 0.03950836 0.03657621 0.04995537 0.04469983 0.04105940 0.04551537 0.05180642

error_pronH 0.0009055807 0.02340388 0.05000653 0.04339390 0.03929627 0.04430251 0.05116947 0.04800552

error_pronH 0.0334304486 0.06134581 0.05018948 0.04400640 0.04857884 0.05530050 0.05134836 0.05356872

error_pronH 0.0671412866 0.04759847 0.03959037 0.04701842 0.05568364 0.05103217 0.05349869 0.06633577

error_pronH 0.0214602097 0.01728465 0.04048314 0.05432041 0.04890278 0.05193941 0.06728899 0.06510779

error_pronH 0.0199991749 0.04827238 0.06183984 0.05388911 0.05639750 0.07222707 0.06924120 0.06718595

error_pronH 0.0729692207 0.07808759 0.06415510 0.06439812 0.08009231 0.07560880 0.07257203 0.07483408

error_pronH 0.0552389712 0.03981395 0.05033203 0.07538330 0.07084374 0.06789383 0.07156330 0.06906939

error_pronH 0.0101622099 0.04823030 0.08075583 0.07418822 0.07001982 0.07401316 0.07095847 0.07125072

error_pronH 0.0635089143 0.09744160 0.08450430 0.07734165 0.08030907 0.07596078 0.07555640 0.07396403

error_pronH 0.1471701781 0.10949240 0.09397442 0.09300588 0.08597771 0.08379617 0.08113644 0.07623440

error_pronH 0.1023847115 0.08070527 0.08495009 0.07743285 0.07689872 0.07435431 0.06939065 0.06493957

error_pronH 0.0833769891 0.08756248 0.07637803 0.07564973 0.07301224 0.06724270 0.06227617 0.06112950

error_pronH 0.1048222591 0.08473302 0.08124879 0.07688083 0.06939527 0.06336327 0.06178965 0.05841032

error_pronH 0.0266498858 0.04711935 0.05482717 0.04827173 0.04317618 0.04416090 0.04167408 0.03904759

error_pronH 0.0530020671 0.06221638 0.05169243 0.04476714 0.04557993 0.04250907 0.03942900 0.04917571

error_pronH 0.0887500842 0.06311191 0.05156113 0.05037165 0.04590345 0.04195031 0.05175833 0.05205893

error_pronH 0.0388291674 0.02796296 0.03736565 0.03424014 0.03078301 0.04705677 0.04843461 0.04533792

error_pronH 0.0049771555 0.03339035 0.02975767 0.02594643 0.04714780 0.04861765 0.04505684 0.04292588

 

$`Errores agregados`

        h1         h2         h3         h4         h5         h6         h7         h8

0.04789196 0.05163872 0.05346425 0.05381539 0.05499571 0.05628511 0.05645721 0.05626474

a.       Funciones

La función error_historico_d1 utiliza la función EvaluacionP, que es en realidad una modificación de la función acurracy del paquete forcast. La misma obtiene una serie de medidas de error (me, mse, sqrt(mse), mae, mpe, mape, r1, theil) de pronóstico para el periodo de validación de las series (test data) (ver referencia en https://otexts.com/fpp2/accuracy.html), al comparar el valor observado (x) con el proyectado (f).

# Solo se reescribe la original para hacerla menos restrictiva

EvaluacionP <- function(f, x){

  # versión 6/12/2019

  n <- length(x)

  test <- 1:n

  error <- (x - f)[test]

  pe <- error / x[test] * 100

 

  me <- mean(error, na.rm = TRUE)

  mse <- mean(error ^ 2, na.rm = TRUE)

  mae <- mean(abs(error), na.rm = TRUE)

  mape <- mean(abs(pe), na.rm = TRUE)

  mpe <- mean(pe, na.rm = TRUE)

 

  fpe <- (c(f[2:n]) / c(x[1:(n - 1)]) - 1)[test - 1]

  ape <- (c(x[2:n]) / c(x[1:(n - 1)]) - 1)[test - 1]

  theil <- sqrt(sum((fpe - ape) ^ 2, na.rm = TRUE) / sum(ape ^ 2, na.rm = TRUE))

 

  # Se coloca este Condicional para poder evaluar pronosticos a 1 paso

  if (length(f)==1){

    r1 <- NA

  } else {

    r1 <- acf(error, plot = FALSE, lag.max = 2, na.action = na.pass)$acf[2, 1, 1]

  }

 

  out <- c(me, mse, sqrt(mse), mae, mpe, mape, r1, theil)

  names(out) <- c("ME", "MSE", "RMSE", "MAE", "MPE", "MAPE","ACF1", "Theil's U")

 

  return(out)

 

}

Para correr la función para cada uno de los periodos seleccionados, se utiliza un pequeño bucle, que modifica los vectores del set de entrenamiento (vistos en la función anterior), para a partir de su longitud evaluar los errores a distintos pasos. Es decir, especificando en el argumento horizonte, los pasos que deseamos evaluar, podemos obtener la evaluación del error a dichos pasos. De manera predeterminada están los periodos 1-6, y el 12.

EvaluacionPh <- function(f, x, horizontes=c(1:6,12)){

  #

  vh <- horizontes

  Resultado<-c()

 

  for(i in vh){

    performance<-as.vector(EvaluacionP(f[1:i], x[1:i]))

    Resultado<-rbind(Resultado, performance)

  }

 

  # Resultados de la función

  row.names(Resultado) <- as.character(horizontes)

  colnames(Resultado) <- c("ME", "MSE", "RMSE", "MAE", "MPE", "MAPE","ACF1", "Theil's U")

 

  return(Resultado)

} 

Finalmente, llegamos a la función de validación cruzada que usamos anteriormente a partir de una ventada recursiva, es decir, que la ventana de estimación se va expandiendo en la medida en que nos desplazamos en el tiempo. La función trabaja de forma muy sencilla:

a.       Segmenta la data en una sección de estimación y otra de validación (test data). Note esta se crean en función de los pasos que queremos evaluar. En el caso de una serie con 100 observaciones en la que deseamos proyectar tres pasos, las primeras 97 serán del set de estimación y las últimas tres del set de validación. Note que esto se logra mediante indexación numérica, por tanto, data estimación: observaciones ength(data.ejercicio) - n.valid; data entrenamiento: test [-c(1:n.train)].

b.       Luego se utiliza la data de estimación o entrenamiento para estimar el modelo ARIMA, incluyendo la dummy, así como la proyección correspondiente.

auto.arima(x.train, xreg=d.train, seasonal=T)

c.       Posteriormente se usa la función EvaluacionPh, para evaluar el pronóstico obtenido en la secuencia de pasos seleccionados (n.valid) y recuperar la medida de error indicada, mediante indexación del vector obtenido ([,medida_error]) en la función de evaluación, llamado resultados.

EvaluacionPh(f=x.hat, x=x.test, h=1:n.valid)[,medida_error]

d.       Note que la selección del set de entrenamiento puede ser móvil o recursiva. La clave de esto es el contador incluido en el bucle (c<-1), dado que las observaciones del set de entrenamiento es x[c:obs] y la evolución de c esta condicionada a al argumento recursiva de la función if (recursiva == FALSE){c <- c+1}, en el caso verdadero no se cumple la condición y c vale 1 en todas las vueltas que de el bucle, así la ventana de estimación va desde 1 hasta obs (que note va creciendo en el bucle que crea los set de entrenamiento (for (obs in (ventana+n.valid):t_obs)). En el caso sea falso, c iria sumando 1 en cada vuelta, dejando una ventana del mismo tamaño en cada vuelta, en vez de ir creciendo.

error_historico_d1 <- function(x,

                               xreg_d,

                               ventana=round(length(!is.na(x))*0.7),

                               n.valid,

                               recursiva=TRUE,

                               medida_error="RMSE",

                               agregar = "mean"){

  # ********************************************************************

  # x = time serie vector

  # ventana = longitud del vector de estimación. Default= 70% de la muestra

  # n.valid = periodos proyectados

  # medida_error:  ME;MSE;RMSE;MAE;MPE;MAPE;ACF1;Theil's U

  # agregar: qué función deseas agregar de los errores.

 

  # Ejemplo: error_historico(AirPassengers, ,ventana=70, n.valid=5, recursiva=TRUE, medida_error="RMSE"                              #    ,agregar = "median")

 

  t_obs <- length(x[!is.na(x)])

  Errores_histo <- c()

  c<-1

 

  for (obs in (ventana+n.valid):t_obs){

   

    # Extrae las observaciones con que se van a trabajar

    data.ejercicio <- x[c:obs]

    dummy.ejercicio <- xreg_d[c:obs]

   

    # Prepara los datos:

    n.train <- length(data.ejercicio) - n.valid  # intervalo de estimación

    x.train <- data.ejercicio[1:n.train]        # Datos del modelo

    x.test <-  data.ejercicio[-c(1:n.train)]

   

    d.train <- dummy.ejercicio[1:n.train]

    d.test <-  dummy.ejercicio[-c(1:n.train)]

   

    # Modelo Autoarima Móvil

    modelo <- auto.arima(x.train, xreg=d.train, seasonal=T)

    x.hat <- forecast::forecast(modelo, xreg=d.test, h=n.valid)$mean

    error_pronH <- EvaluacionPh(f=x.hat, x=x.test, h=1:n.valid)[,medida_error]

   

    Errores_histo <- rbind(Errores_histo,error_pronH)

   

    if (recursiva == FALSE){c <- c+1}   #"Movil"

   

  } #Cierrra el bucle

 

  colnames(Errores_histo) <- paste("h",1:n.valid, sep="")

 

  errores_h <- apply(Errores_histo, MARGIN=2, FUN=get(agregar), na.rm = TRUE)

 

  info <- rbind(paste("Medida de errro: ", medida_error, sep=""),

                paste("Horizonte: ", n.valid, sep=""),

                paste("Función agregar: ", agregar, sep="") )

 

  resultados<-list("Informacion"=info,

                   "Modelo Arima"=modelo,

                   "Hitorico errores"=Errores_histo,

                   "Errores agregados"=errores_h)

  return(resultados)

} #Cierra la función

Entradas populares de este blog

Ejemplos de tablas en Stata

¿Cómo hacer pruebas de hipótesis en STATA?

Métodos y técnicas de Investigación Económica