LessThan2Periods

Vásquez Guerra Carlos Fernando

12/4/2020

Tomemos los datos del coronavirus de Our World in Data a la fecha de creación de este documento.

data <- read_csv("https://covid.ourworldindata.org/data/owid-covid-data.csv", col_names = TRUE)

Bajo la necesidad de crear un análisis de series temporales con esta información, se tomaron los datos para México.

(data_Mexico <- data %>% dplyr::filter(location == "Mexico") %>% dplyr::select(date, total_cases) %>% filter(!is.na(total_cases)))

Ahora vamos a crear un objeto de ts en el cual se colocará una frecuencia de 365, ya que tenemos datos diarios, además obtengamos una descomposición natural entre la tendencia, estacionalidad y el remanente.

tseries <- ts(data_Mexico$total_cases, 
              start = c(2020, as.numeric(format(data_Mexico$date[1], "%j"))), 
              frequency = 365)
# stl(tseries)

Al ejecutar el comando stl(tseries), se obtendría el siguiente error:

Error in stl(tseries) : series is not periodic or has less than two periods

Este error se debe a que no tenemos ni un ciclo (son necesarios al menos dos para obtener una estimación sobre la estacionalidad) ya que contamos con 280 días de un año. Para solucionar este inconveniente no utilizaremos la descomposión proporcionada por stl() o decompose(). Utilizaremos un estimación de la estacionalidad basada en series de Fourier, tal como se muestra en el post Seasonal decomposition of short time series. Donde se utiliza el hecho de que es posible utilizar un modelo de regresión lineal para descomponer una serie de tiempo, así como se explica en el siguiente bookdown: Forecasting: Principles and Practice.

Para estos en particular, consideraremos que los datos tienen una estacionalidad semanal, por lo que al tener 52 semanas, el parámetro \(K\) mostrado en dicho modelo debe cumplir que se \(K<=52/2 = 26\), tal como lo establece la frecuencia Nyquist. Para seleccionar el mejor \(K\), es decir la cantidad de pares de funciones \(sen\left(\frac{2\pi tk}{m}\right)\) y \(cos\left(\frac{2\pi tk}{m}\right)\), utilizaremos el modelo que obtenga el mínimo valor de la función CV() (leave-one-out CV statistic).

library(forecast)

k_min <- function(ts, kmax){
  season_f <- map(1:kmax, ~ fourier(ts, .x)) #Obtenemos las apxomiaciones de estacionalidad por coef. de Fourier.
  lineal_ajusts <- map(season_f, ~tslm(ts ~ trend + .x)) #Creamos los ajuste lineales
  coef_fourier <- map(lineal_ajusts, ~CV(.x)) #Aplicamos la función CV a los anteriores modelos
  coef_fourier <- map_dbl(coef_fourier, ~'[['(.x, 1)) #Obtenemos los resultados interesantes
  which(coef_fourier == min(coef_fourier)) #Determinamos cual es el k que minimiza el CV
}

(min_casesT <- k_min(tseries, 26))

## [1] 20

Así que tomaremos el modelo anterior para realizar la correspondiente descomposición.

fitted_ts <- tslm(tseries ~ trend + fourier(tseries, min_casesT))

trend <- coef(fitted_ts)[1] + coef(fitted_ts)['trend']*seq_along(data_Mexico$total_cases)

components <- cbind(
  data = (data_Mexico$total_cases),
  trend = trend,
  season = (data_Mexico$total_cases) - trend - residuals(fitted_ts),
  remainder = residuals(fitted_ts)
)

autoplot(components, facet=TRUE)

Aunque como siempre podemos dar un mejor formato

Mexico_series <- tibble(time = data_Mexico$date,
       Datos = (data_Mexico$total_cases),
       Tendencia = trend,
       Estacionalidad = (data_Mexico$total_cases) - trend - residuals(fitted_ts),
       Remanente = residuals(fitted_ts)) %>% 
  gather("type", "value", 2:5) %>% 
  mutate(type = factor(type, levels = c("Datos", "Tendencia", "Estacionalidad", "Remanente"))) %>% 
  ggplot(aes(x = time, y = value)) + 
  geom_line() + 
  facet_grid(type~., scales = "free_y") +
  labs(x = "Tiempo", y = "Población") + 
  ggtitle("Descomposición temporal. Casos totales para México")

library(plotly)
ggplotly(Mexico_series)

Aunque el caso anterior no tiene mucho sentido, ya que al estar tratando con casos acumulados sabemos que los valores siguientes seguirán en aumento, por lo que una componente no se está haciendo presente. Así que, manera análoga, se puede tener esta descomposición con la cantidad de casos nuevos diarios

data_Mexico_dailyCases <- data %>% dplyr::filter(location == "Mexico") %>% dplyr::select(date, new_cases) %>% filter(!is.na(new_cases))

tseries <- ts(data_Mexico_dailyCases$new_cases, 
              start = c(2020, as.numeric(format(data_Mexico_dailyCases$date[1], "%j"))), 
              frequency = 365)

fitted_ts <- tslm(tseries ~ trend + fourier(tseries, k_min(tseries, 26)))

trend <- coef(fitted_ts)[1] + coef(fitted_ts)['trend']*seq_along(data_Mexico_dailyCases$new_cases)

Mexico_series_daily <- tibble(time = data_Mexico_dailyCases$date,
       Datos = (data_Mexico_dailyCases$new_cases),
       Tendencia = trend,
       Estacionalidad = (data_Mexico_dailyCases$new_cases) - trend - residuals(fitted_ts),
       Remanente = residuals(fitted_ts)) %>% 
  gather("type", "value", 2:5) %>% 
  mutate(type = factor(type, levels = c("Datos", "Tendencia", "Estacionalidad", "Remanente"))) %>% 
  ggplot(aes(x = time, y = value)) + 
  geom_line() + 
  facet_grid(type~., scales = "free_y") +
  labs(x = "Tiempo", y = "Población") + 
  ggtitle("Descomposición temporal. Casos nuevos diarios para México")


ggplotly(Mexico_series_daily)

En el ajuste anterior \(K\) = 5. El problema aquí es la suposición de un comportamiento lineal en la tendencia. Vamos a considerar entonces un comportamiento no lineal en esta. Para más información véase el siguiente capítulo del anterior libro mencionado.

t <- time(tseries)
#Puntos de corte
# data_Mexico_dailyCases %>% filter(!is.na(new_cases)) %>% mutate(index = row_number()) %>% filter(date == "2020-5-5")
t.break1 <- time(tseries)[104]
# data_Mexico_dailyCases %>% filter(!is.na(new_cases)) %>% mutate(index = row_number()) %>% filter(date == "2020-8-1")
t.break2 <- time(tseries)[192]
# data_Mexico_dailyCases %>% filter(!is.na(new_cases)) %>% mutate(index = row_number()) %>% filter(date == "2020-10-4")
t.break3 <- time(tseries)[256]
# data_Mexico_dailyCases %>% filter(!is.na(new_cases)) %>% mutate(index = row_number()) %>% filter(date == "2020-11-17")
t.break4 <- time(tseries)[300]

#Función pair-wise
tb1 <- ts(pmax(0, t - t.break1), start = c(2020, as.numeric(format(data_Mexico_dailyCases$date[1], "%j"))), frequency = 365)
tb2 <- ts(pmax(0, t - t.break2), start = c(2020, as.numeric(format(data_Mexico_dailyCases$date[1], "%j"))), frequency = 365)
tb3 <- ts(pmax(0, t - t.break3), start = c(2020, as.numeric(format(data_Mexico_dailyCases$date[1], "%j"))), frequency = 365)
tb4 <- ts(pmax(0, t - t.break4), start = c(2020, as.numeric(format(data_Mexico_dailyCases$date[1], "%j"))), frequency = 365)

#Modelo lineal
fit.pw <- tslm(tseries ~ t + tb1 + tb2 + tb3 + tb4)
fit.pw <- fitted(fit.pw)

De acuerdo a lo anterior, el k_min(tseries) seguirá siendo al anterior a este nuevo modelo (\(k\) = 5). Estableciendo los puntos de corte (las fechas donde parece haber un gran cambio de tendencia) en 2020-10-4 y 2020-11-17.

fitted_ts <- tslm(tseries ~ fit.pw + fourier(tseries, k_min(tseries, 26)))
Mexico_series_daily <- tibble(time = data_Mexico_dailyCases$date,
       Datos = (data_Mexico_dailyCases$new_cases),
       Tendencia = as.vector(fit.pw),
       Estacionalidad = (data_Mexico_dailyCases$new_cases) - as.vector(fit.pw) - residuals(fitted_ts),
       Remanente = residuals(fitted_ts)) %>% 
  gather("type", "value", 2:5) %>% 
  mutate(type = factor(type, levels = c("Datos", "Tendencia", "Estacionalidad", "Remanente"))) %>% 
  ggplot(aes(x = time, y = value)) + 
  geom_line() + 
  facet_grid(type~., scales = "free_y") +
  labs(x = "Tiempo", y = "Población") + 
  ggtitle("Descomposición temporal. Casos nuevos diarios para México")

ggplotly(Mexico_series_daily)

Vamos a quitar los datos en 0s.

data_Mexico_dailyCases_N0 <- data_Mexico_dailyCases %>% filter(new_cases != 0)
tseries <- ts(data_Mexico_dailyCases_N0$new_cases, 
              start = c(2020, as.numeric(format(data_Mexico_dailyCases_N0$date[1], "%j"))), 
              frequency = 365)

t <- time(tseries)
#Puntos de corte
# data_Mexico_dailyCases_N0 %>% filter(!is.na(new_cases)) %>% mutate(index = row_number()) %>% filter(date == "2020-5-5")
t.break1 <- time(tseries)[61]
# data_Mexico_dailyCases_N0 %>% filter(!is.na(new_cases)) %>% mutate(index = row_number()) %>% filter(date == "2020-8-1")
t.break2 <- time(tseries)[149]
# data_Mexico_dailyCases_N0 %>% filter(!is.na(new_cases)) %>% mutate(index = row_number()) %>% filter(date == "2020-10-4")
t.break3 <- time(tseries)[213]
# data_Mexico_dailyCases_N0 %>% filter(!is.na(new_cases)) %>% mutate(index = row_number()) %>% filter(date == "2020-11-17")
t.break4 <- time(tseries)[257]

#Función pair-wise
tb1 <- ts(pmax(0, t - t.break1), start = c(2020, as.numeric(format(data_Mexico_dailyCases_N0$date[1], "%j"))), frequency = 365)
tb2 <- ts(pmax(0, t - t.break2), start = c(2020, as.numeric(format(data_Mexico_dailyCases_N0$date[1], "%j"))), frequency = 365)
tb3 <- ts(pmax(0, t - t.break3), start = c(2020, as.numeric(format(data_Mexico_dailyCases_N0$date[1], "%j"))), frequency = 365)
tb4 <- ts(pmax(0, t - t.break4), start = c(2020, as.numeric(format(data_Mexico_dailyCases_N0$date[1], "%j"))), frequency = 365)

#Modelo lineal
fit.pw <- tslm(tseries ~ t + tb1 + tb2 + tb3 + tb4)
fit.pw <- fitted(fit.pw)

fitted_ts <- tslm(tseries ~ fit.pw + fourier(tseries, k_min(tseries, 26)))
Mexico_series_daily <- tibble(time = data_Mexico_dailyCases_N0$date,
       Datos = (data_Mexico_dailyCases_N0$new_cases),
       Tendencia = as.vector(fit.pw),
       Estacionalidad = (data_Mexico_dailyCases_N0$new_cases) - as.vector(fit.pw) - residuals(fitted_ts),
       Remanente = residuals(fitted_ts)) %>% 
  gather("type", "value", 2:5) %>% 
  mutate(type = factor(type, levels = c("Datos", "Tendencia", "Estacionalidad", "Remanente"))) %>% 
  ggplot(aes(x = time, y = value)) + 
  geom_line() + 
  facet_grid(type~., scales = "free_y") +
  labs(x = "Tiempo", y = "Población") + 
  ggtitle("Descomposición temporal. Casos nuevos diarios para México")

ggplotly(Mexico_series_daily)