Modelos de frecuencia mixta en R

 1.    Simulando la data

La siguiente entrada muestra un ejemplo de como estimar un modelo de frecuencia mixta. Puntualmente simulamos una serie de consumo trimestral que queremos estimar en una regresión contra una serie de ventas mensuales. Como tal, la diferencia en el numero de observaciones limita la posibilidad de estimar directamente un modelo de regresión tradicional -aquí suponemos que las series son estacionarias, proceso obligatorio para estimar una regresión entre series temporales-. Todo el ejercicio se realiza en Package ‘midasr’ (Zemlys, 2022).

En primer lugar, simulamos la serie de consumo:

# Simulando series
library(ggplot2)
library(forecast)
library(midasr)
 
# Simulando consumo
# Definición de parámetros
set.seed(1234)
trend <- 0.5  # Tendencia de la serie
seasonality <- 0.2  # Estacionalidad trimestral
n <- 22*4  # Número de trimestres a simular (22 años)
 
# Generación de la serie trimestral
consumo <- rnorm(n, mean = trend + seasonality * rep(1:4, length.out = n))
consumo_trimestral <- ts(consumo, frequency = 4, start = 2000)
 
autoplot(consumo_trimestral)

 

Posteriormente simulamos la serie de ventas, teniendo presente que el número de observaciones de la variable mensual, debe coincidir con las necesarias para equiparar las observaciones de mi serie en frecuencia más alta. Puntualmente, si tenemos una serie trimestral de 10 años, entonces necesitaríamos 4*10 observaciones de mi serie en menor frecuencia.

# simulando ventas
# Definición de parámetros
trend <- 0.8  # Tendencia de la serie
seasonality <- c(0.5, 0.3, 0.2)  # Estacionalidad mensual (ejemplo: enero, febrero, marzo)
n2 <- (n/4)*12  # Número de meses a simular
 
# Generación de la serie mensual
ventas_mensuales <- numeric(n)
for (i in 1:n2) {
  ventas_mensuales[i] <- rnorm(1, mean = trend + seasonality[(i - 1) %% length(seasonality) + 1])
}
 
ventas_mensuales <- ts(ventas_mensuales, frequency = 12, start = 2000)
 
autoplot(ventas_mensuales)

2.     Estimando el modelo midas

Ahora que hemos simulado nuestras series, podemos estimar nuestro modelo en frecuencia mixta. Primero utilizamos una ventana de datos para poder comparar los valores proyectados con los valores observados. Luego usamos la función midas_r para estimar el modelo en frecuencia mixta. Note aquí que la clase está en la función fmls que realiza un collapse de la data en alta frecuencia para que este coincida en el total de observaciones. Dentro de esta función se coloca la variable independiente, el numero de rezagos a incluir, la cantidad de veces que esta contenida mi serie de baja frecuencia en la de alta (trimestral y mensual 3; anual y mensual 12…)…

 # train data
x_q <-  window(ventas_mensuales, end = c(2019, 12))
y_c <-  window(consumo_trimestral, end = c(2019, 4))
aa <- round(length(x_q)/length(y_c))
 
trend <- 1:length(y_c)
 
mr <- midas_r(y_c ~ trend + fmls(x_q, 11, aa, nealmon), start = list(x_q = rep(0, 3)))
 
MIDAS regression model with "ts" data:
Start = 2000(4), End = 2019(4)
 model: y_c ~ trend + fmls(x_q, 11, aa, nealmon)
(Intercept)       trend        x_q1        x_q2        x_q3
    0.50057     0.01001    -0.12472    20.67112    -0.84271
 
Function optim was used for fitting

 [Nota: note que la funciones fmls selecciona volares puntuales de una vector de datos en el caso de colocar 0 toma la observación actual. 1, toma la actual y crea un vector del tamaño de la serie de baja frecuencia que adhiere el vector anterior. 2, a los dos vectores anteriores, agrega uno con el segundo rezago. Finalmente, el último argumento le indica cada cuantas observaciones de tu serie de alta frecuencia tu tomaras un dato para crear un vector cuya longitud coincida con el de baja frecuencia.

x<-c(1:12)
x
 [1]  1  2  3  4  5  6  7  8  9 10 11 12
 
fmls(x, 0, 3)
     X.0/m
[1,]     3
[2,]     6
[3,]     9
[4,]    12
 
fmls(x, 1, 3)
     X.0/m X.1/m
[1,]     3     2
[2,]     6     5
[3,]     9     8
[4,]    12    11
]

Ahora, solo témenos que considerar cuales son los valores futuros de las variables independientes. Esto con el fin de poder crear pronósticos condicionados.

##New trend values
n.valid <- 8
xn <- forecast::snaive(x_q,aa*n.valid)$mean
trendn <- length(y_c) + 1:n.valid

Finalmente, usamos la función forecast para obtener los valores predichos para mi variable independiente.

fh_midas_1 <- forecast(mr,  list(trend = trendn, x_q = xn), method = "static")
fh_midas_1
         Qtr1     Qtr2     Qtr3     Qtr4
2020 1.213594 1.255759 1.216617 1.158818
2021 1.253628 1.295793 1.256651 1.198851

 3.     Estimando un modelo dinámico (inluyendo un ar de y)

 En el siguiente ejemplo usamos mls para ingresar en el modelo valores rezagados de la variable dependiente.

 mr.dyn <- midas_r(y_c ~ trend +
                    mls(y_c, 1:2, 1, "*") +
                    fmls(x_q, 12, aa, nealmon),
                    start = list(x_q = rep(0, 3)))
 
MIDAS regression model with "ts" data:
Start = 2001(3), End = 2019(4)
 model: y_c ~ trend + mls(y_c, 1:2, 1, "*") + fmls(x_q, 18, aa, nealmon)
(Intercept)       trend        y_c1        y_c2        x_q1        x_q2        x_q3
   0.531114    0.009247    0.105423   -0.017154   -0.214118    2.527028   -0.111902
 
Function optim was used for fitting
 
fh_midas_2 <- forecast(mr.dyn,  list(trend = trendn, x_q = xn), method = "dynamic")
fh_midas_2
         Qtr1     Qtr2     Qtr3     Qtr4
2020 1.223314 1.274915 1.186098 1.196951
2021 1.278432 1.315320 1.226402 1.237496

 4.     Comparando diversas estrategias de pronósticos

 Finalmente armamos una base de datos, un poco tosca porque seguro hay formas de hacer esto de manera mas eficiente y elegante. Adicionalmente, agregamos una alternativa de estimación mediante el modelo auto.arima. es importante verificar que mayor complejidad no siempre significa mejor precisión en la estimación de pronósticos.

 # AutoArima
fh_arima <- forecast::forecast(auto.arima(y_c, stepwise = FALSE))
 
### Graficando
obs_data <- window(consumo_trimestral, start = c(2020, 1))
 
data_fores <- data.frame(
  y_obs = c(y_c, rep(NA, n.valid)),
  f.midas1 = c(rep(NA, length(y_c)-1), tail(y_c,1), fh_midas_1$mean),
  f.midas2 = c(rep(NA, length(y_c)-1), tail(y_c,1), fh_midas_2$mean),
  fh_arima = c(rep(NA, length(y_c)-1), tail(y_c,1), fh_arima$mean),
  obs = c(rep(NA, length(y_c)-1), tail(y_c,1), obs_data),
  fecha = time(consumo_trimestral)
)
 
data_fores |> tail(8*4) |>
ggplot() +
  geom_line(aes(fecha,y_obs,colour="Obs"), size=1) +
  geom_line(aes(fecha, obs,colour="Obs1"),linetype = 2, size=1) +
  geom_line(aes(fecha, f.midas1,colour="xMidas"),linetype = 3, size=1) +
  geom_line(aes(fecha, f.midas2,colour="xArMidas"), linetype = 3, size=1) +
  geom_line(aes(fecha, fh_arima,colour="Arima"),linetype = 3, size=1) +
  theme_classic() +
  scale_color_manual(name = "Serie", values = c("Obs" = "black", "Obs1" = "black", "xMidas" = "darkred", "xArMidas" = "darkblue", "Arima"="darkgreen")) +
  theme(legend.position = "bottom")

Referencias

Zemlys, V. (2022). Mixed Data Sampling Regression. Package ‘midasr’. R CRAN.