Análisis Basico de series de tiempo univariadas
Paquete estadistico proyecto R
Alejandro García
2021-08-28
1 INTRODUCCIÓN
El motivo de este trabajo es una introducción a los modelos ARIMA en el estudio de serie temporales. Recordando que para dicho estudio es necesario que la serie temporal sea estacionaria, es decir que su media, varianza y autocovarianza independiente del tiempo.
Para realizar el mismo se utilizó proyecto R, con el entorno gráfico RSTUDIO. A través del cual se aplico las diferentes herramientas necesarias.
El trabajo comenzara describiendo las configuraciones principales para los “chunk” y las librerías utilizadas; luego pasará a describir los datos a utilizar en este trabajo con sus correspondientes transformaciones para llegar a su correcto análisis. Terminado con lo anterior se procede con la metodología de Box y Jenkins (Box, Jenkins, 1976) con el objetivo de obtener el mejor ajuste sobre los datos y los pronósticos correctos. Para terminar se presenta una conclusión sobre el análisis y resultados del trabajo; sin olvidar una pequeña descripción de los materiales utilizados durante el trabajo.
2 CONFIGURACIONES Y LIBRERÍAS IMPORTANTES
Las siguientes configuraciones corresponden a lo que se dejara por defecto sobre todos los chunk correspondientes a código en R.
knitr::opts_chunk$set(echo = T, # muestra codigo
eval = T, # ejecuta comando
message= T, # muestra mensajes
warning = T, # muestra advertencias
results = 'asis') # asistente de resultadoAdemás las librerías necesarias para este trabajo son:
library(readr) # base r desde la web
library(dplyr) #filter
library(tidyverse) #filter
library(tseries) # test adf
library(forecast) # auto.arima
library(ggplot2) #graficos
library(graphics) # para añadir lineas al grafico
library(dynlm) # paquete de regresiones3 DATOS
Para proceder con el análisis temporal, utilizaré la base de datos correspondiente a los casos diarias de coronavirus, serie histórica para argentina. A partir del repositorio creado por Sistemas Mapache con el objetivo de poder contar con datos abiertos de la información oficial proveniente de los partes diarios sobre la situación de COVID-19 en Argentina https://github.com/SistemasMapache/Covid19arData.
basetf <-read_csv('https://docs.google.com/spreadsheets/d/16-bnsDdmmgtSxdWbVMboIHo5FRuz76DBxsz_BbsEVWA/export?format=csv&id=16-bnsDdmmgtSxdWbVMboIHo5FRuz76DBxsz_BbsEVWA&gid=0')
base_salta <- basetf %>%
filter(., osm_admin_level_4 == "Salta" ) %>%
select(., "fecha", "nue_casosconf_diff") %>%
rename(., casos_d = nue_casosconf_diff)
base_salta$fecha = as.Date(base_salta$fecha, format="%d/%m/%Y")
base_salta <- base_salta %>%
mutate(.,
casos_a = cumsum(0 + casos_d),
logcasos_a = log(casos_a)) %>%
filter(., fecha<="2021-06-01" & fecha >= "2021-01-01") Luego de cargar y filtrar los datos de manera de retener solamente los correspondientes a “Salta”, se genero la variable stock de casos acumulados, sobre la cual se aplicará logaritmos para mejorar la forma exponencial de la función en nivel y evitar problemas en la dispersión de datos.
El periodo seleccionado para el análisis de los contagios corresponde a la primera mitad de año, es decir desde el primero de enero del 2021 hasta el 1 de junio del mismo. Para en este primer momento realizar pronósticos a corto plazo y evitar tener en cuenta el efecto de las restricciones.
Los siguientes gráficos corresponden a la función a la función casos acumulados y casos acumulados en logaritmos.
ggplot(data = base_salta,
aes(x = fecha,
y = casos_a)) +
geom_line()
Dada esta forma exponencial es que se procede a observar que sucede con su transformación logarítmica.
ggplot(data = base_salta,
aes(x = fecha,
y = logcasos_a)) +
geom_line()
A pesar de tener los datos en logaritmos este tramos de la función mantiene una forma creciente que no tiende a aplanarse.
3.1 CONFIGURANDO SERIE TEMPORAL
Configurando las objetos correspondiente a las series temporales para su análisis. Dado que se quiere estudiar una serie temporal, en este caso univariada, es necesario que esta sea estacionaria.
# Objeto con logaritmos
casos_ats <- ts(base_salta$casos_a,
start = c(1,1),
frequency = 7)
plot(casos_ats)
plot(log(casos_ats))
log(casos_ats) == ts(base_salta$logcasos_a,
start = c(1,1),
frequency = 7)Time Series: Start = c(1, 1) End = c(22, 5) Frequency = 7 [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE [16] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE [31] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE [46] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE [61] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE [76] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE [91] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE [106] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE [121] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE [136] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE [151] TRUE TRUE
Se graficó dentro del modelo temporal nuevamente ambas funciones y se configuro como una serie acumulada diaria con frecuencia semanal para captar los ciclos estacionales de tamaño semanal. Dado que los datos son seleccionados desde el 1 de enero del 2021, día viernes, día uno sera viernes, día dos sábado hasta llegar a día siete en jueves. .Claramente los gráficos no son estables con el tiempo en su media. También se comparo que los datos logarítmicos obtenidos de formas diferentes son idénticos.
Dado que el principal objetivo es contar con una serie estacionaria se procede a observar los test de raíces unitarias sobre los datos y las correspondientes funciones de autocorrelación para la función de nivel casos acumulados en logaritmos sobre el total. Ademas se enseñara el modelo mas adecuados para los datos según el procedimiento auto.arima como una guía para continuar el presente estudio.
Acf es la función que captura todos los retardos posibles (valor q de un ARIMA), mientras Pacf describe la correlación de todos los puntos que son k pasos retardados o separados (valor p de un ARIMA).
# Primera inspección
adf.test(log(casos_ats))## Warning in adf.test(log(casos_ats)): p-value greater than printed p-valueAugmented Dickey-Fuller Testdata: log(casos_ats) Dickey-Fuller = 4.3659, Lag order = 5, p-value = 0.99 alternative hypothesis: stationary
modelog <- auto.arima(log(casos_ats))
modelogSeries: log(casos_ats) ARIMA(1,1,1)(0,1,1)[7]
Coefficients: ar1 ma1 sma1 0.9853 -0.7303 -0.7552 s.e. 0.0206 0.0844 0.0609
sigma^2 estimated as 8.377e-06: log likelihood=717.76 AIC=-1427.53 AICc=-1427.24 BIC=-1415.65
# Funciones de autocorrelación
par(mfrow=c(1,2))
Acf(log(casos_ats), main="")
Pacf(log(casos_ats), main="")
# Descomposición
par(mfrow=c(1,1))
plot(decompose(log(casos_ats)), xlab='Semana')
boxplot(log(casos_ats) ~ cycle(log(casos_ats)),
main="Box - Plot sobre Serie logaritmica",
ylab="Log casos acumulados",
xlab="Ciclo")
cycle(log(casos_ats))Time Series: Start = c(1, 1) End = c(22, 5) Frequency = 7 [1] 1 2 3 4 5 6 7 1 2 3 4 5 6 7 1 2 3 4 5 6 7 1 2 3 4 5 6 7 1 2 3 4 5 6 7 1 2 [38] 3 4 5 6 7 1 2 3 4 5 6 7 1 2 3 4 5 6 7 1 2 3 4 5 6 7 1 2 3 4 5 6 7 1 2 3 4 [75] 5 6 7 1 2 3 4 5 6 7 1 2 3 4 5 6 7 1 2 3 4 5 6 7 1 2 3 4 5 6 7 1 2 3 4 5 6 [112] 7 1 2 3 4 5 6 7 1 2 3 4 5 6 7 1 2 3 4 5 6 7 1 2 3 4 5 6 7 1 2 3 4 5 6 7 1 [149] 2 3 4 5
El estadístico de Dickey y Fuller (Dickey, Fuller, 1979), se calculo a través del comando adf.test, para prueba de raíz unitaria o mas conocido como test ADF. Presenta un valor-p mayor al nivel de significancia tradicional, \(\alpha = 0,05\). Esto me permite quedar en la zona de no rechazo de hipótesis nula y por tanto indicaría la presencia de al menos una raíz unitaria. Además la función ACF esta descendiendo lentamente no llegando a cero en la gráfica, lo que significa que no existe asociación lineal entre las observaciones separada por k retardos; por lo tanto se advierte presencia de una raíz unitaria.
La descomposición muestra una tendencia similar al comportamiento de la serie, todavía no es posible observar un comportamiento estacional sobre la serie original en el modelo de box - plot.
Paso a calcular la primera diferencias de la serie en logaritmos y sus correspondientes funciones de autocorrelación.
# Primera diferencia
head(diff(log(casos_ats)))Time Series: Start = c(1, 2) End = c(1, 7) Frequency = 7 [1] 0.0005358100 0.0003570472 0.0009366429 0.0011139084 0.0022241013 [6] 0.0030611590
par(mfrow=c(1,1))
plot(diff(log(casos_ats)),
main="Primera diferencia de la serie log",
ylab="Log casos acumulados",
xlab="Tiempo")
# Prueba estacionariedad
adf.test(diff(log(casos_ats), lag =1))Augmented Dickey-Fuller Testdata: diff(log(casos_ats), lag = 1) Dickey-Fuller = -2.0594, Lag order = 5, p-value = 0.5516 alternative hypothesis: stationary
par(mfrow=c(1,2))
Acf(diff(log(casos_ats), lag =1), main="")
Pacf(diff(log(casos_ats), lag =1), main="")
Nuevamente se realizó la prueba de raíz unitaria a través del test ADF, obteniendo en este caso valor-p mayor a el alfa de significancia, por lo tanto aún no puedo rechazar la hipótesis nula de raíz unitaria.
Considerando los valores estacionales se observa sobre la función ACF desciende lentamente con lo que estoy advirtiendo una posibles raíz unitaria estacional.
Se procede a observar el comportamiento estacional en la serie, y en caso de considerarse necesario, aplicar la desestacionalización.
# Descomposición de la serie en 1º dif
casos_logts.desc = decompose(diff(log(casos_ats), lag =1))
plot(casos_logts.desc, xlab='Semana')
# Box-plot y ciclo sobre la serie en 1º dif
boxplot(diff(log(casos_ats), lag =1) ~ cycle(diff(log(casos_ats), lag =1)),
main="Box - Plot sobre Serie log con una diferencia",
ylab="Log casos acumulados",
xlab="Ciclo")
cycle(diff(log(casos_ats), lag =1))Time Series: Start = c(1, 2) End = c(22, 5) Frequency = 7 [1] 2 3 4 5 6 7 1 2 3 4 5 6 7 1 2 3 4 5 6 7 1 2 3 4 5 6 7 1 2 3 4 5 6 7 1 2 3 [38] 4 5 6 7 1 2 3 4 5 6 7 1 2 3 4 5 6 7 1 2 3 4 5 6 7 1 2 3 4 5 6 7 1 2 3 4 5 [75] 6 7 1 2 3 4 5 6 7 1 2 3 4 5 6 7 1 2 3 4 5 6 7 1 2 3 4 5 6 7 1 2 3 4 5 6 7 [112] 1 2 3 4 5 6 7 1 2 3 4 5 6 7 1 2 3 4 5 6 7 1 2 3 4 5 6 7 1 2 3 4 5 6 7 1 2 [149] 3 4 5
La descomposición de la serie muestra un fuerte comportamiento estacional, además notar que la tendencia se muestra como una función escalonada y suavizada de los ciclos de los datos.
Sobre la función box-plot planteada para observar el comportamiento estacional, se observa que los valores mas bajos son una constante del día domingo, mientras que lo datos mas grandes se encuentran presentes entre martes y viernes.
Se calcula una diferenciación estacional sobre la serie logarítmica en primera diferencia:
# Diferencia estacional
head(diff(diff(log(casos_ats) ,lag =1),lag=7))Time Series: Start = c(2, 2) End = c(2, 7) Frequency = 7 [1] 0.0026872220 0.0005682056 0.0012189570 0.0005107160 0.0007110858 [6] -0.0005708241
par(mfrow=c(1,1))
plot(diff(diff(log(casos_ats) ,lag =1),lag=7),
main="Serie log con diferencia primera y estacional",
ylab="Log casos acumulados",
xlab="Tiempo")
# Prueba raíz unitaria
adf.test(diff(diff(log(casos_ats) ,lag =1),lag=7))## Warning in adf.test(diff(diff(log(casos_ats), lag = 1), lag = 7)): p-value
## smaller than printed p-valueAugmented Dickey-Fuller Testdata: diff(diff(log(casos_ats), lag = 1), lag = 7) Dickey-Fuller = -4.7095, Lag order = 5, p-value = 0.01 alternative hypothesis: stationary
par(mfrow=c(1,2))
Acf(diff(diff(log(casos_ats) ,lag =1),lag=7), main = "")
Pacf(diff(diff(log(casos_ats) ,lag =1),lag=7), main = "")
En primer lugar notar que por primera vez el test ADF, permite rechazar la hipótesis nula de raíz unitaria. Por lo tanto se logro quitar los coeficientes que estaban fuera del circulo unitario.
Teniendo en cuenta que la serie ya no cuenta con raíces unitaria, paso a plantear una descomposición de la nueva serie, y observar su comportamiento.
# Descomposición de la serie "estacionaria"
casos_logts.descf = decompose(diff(diff(log(casos_ats) ,lag =1),lag=7))
plot(casos_logts.descf, xlab='Semana')
# Box-plot y ciclo sobre la serie con dif_1 + dif_S
boxplot(diff(diff(log(casos_ats) ,lag =1),lag=7) ~ cycle(diff(diff(log(casos_ats) ,lag =1),lag=7)),
main="Box - Plot sobre Serie log con Dif_1 y Dif_S",
ylab="Log casos acumulados",
xlab="Ciclo")
cycle(diff(diff(log(casos_ats) ,lag =1),lag=7))Time Series: Start = c(2, 2) End = c(22, 5) Frequency = 7 [1] 2 3 4 5 6 7 1 2 3 4 5 6 7 1 2 3 4 5 6 7 1 2 3 4 5 6 7 1 2 3 4 5 6 7 1 2 3 [38] 4 5 6 7 1 2 3 4 5 6 7 1 2 3 4 5 6 7 1 2 3 4 5 6 7 1 2 3 4 5 6 7 1 2 3 4 5 [75] 6 7 1 2 3 4 5 6 7 1 2 3 4 5 6 7 1 2 3 4 5 6 7 1 2 3 4 5 6 7 1 2 3 4 5 6 7 [112] 1 2 3 4 5 6 7 1 2 3 4 5 6 7 1 2 3 4 5 6 7 1 2 3 4 5 6 7 1 2 3 4 5
Se observa que persiste un comportamiento “estacional” en la serie, y un comportamiento irregular variable con t.
Sin embargo no es claro en el gráfico de box-plot, dado que parece haberse removido correctamente el comportamiento estacional.
4 PROCEDIMIENTO DE BOX - JENKINS
Selección del modelo adecuado con el procedimiento de BOX Y JENKINS (Box, Jenkins, 1976).
4.0.1 SUPUESTOS DEL MODELO ARIMA
Los datos deben ser estacionarios; por estacionario, significa que las propiedades de la serie no dependen del momento en que se capturan. Una serie de ruido blanco y series con comportamiento cíclico también pueden considerarse como series estacionarias.
4.0.2 PASOS A SEGUIR PARA EL MODELADO ARIMA
Ahora este nuevo procedimiento consta de pasos en la selección adecuada del modelo:
Identificación del Modelo.
Ajustar el modelo.
Diagnostico del modelo.
Pronósticos.
Es decir, este procedimiento requiere que previamente la serie cumpla las propiedades de estacionariedad, objetivo que se consiguió y diagnosticó en la unidad anterior.
4.1 IDENTIFICACIÓN DEL MODELO
La IDENTIFICACIÓN consisten en la inspección de los datos y de manera formal platear las funciones logarítmica original y sus transformaciones, es decir las respectivas funciones de autocorrelación simple y parcial.
# Funciones de autocorrelación
Acf(diff(diff(log(casos_ats) ,lag =1),lag=7),
main="",sub="Figura: ACF" ,lag=40)
Pacf(diff(diff(log(casos_ats) ,lag =1),lag=7),
main = '',sub="Figura: PACF", lag=40)
Notar que la función de ACF, muestra como posibles rezados en 1 y en 7 y una forma descendiente oscilatoria partir del rezago 7. Mientras que la función de PACF muestra como posibles retardo en 1, 7, 8 y 14.
Sin embargo dada la dificultad a la hora de identificar el modelos solo por la inspección de la serie se procederá a realizar también un ajuste automático.
4.2 AJUSTE DEL MODELO
Aplicando la función auto.arima se calculó anteriormente el mejor modelo para los datos. Ahora también teniendo en cuenta los datos observados en las funciones de autocorrelación se ajustaran estos distintos modelos.
# Modelando
modelog <- auto.arima(log(casos_ats),
stepwise = F,
method = "ML",
# trace = T,
seasonal = T)
modelogSeries: log(casos_ats) ARIMA(1,1,1)(0,1,1)[7]
Coefficients: ar1 ma1 sma1 0.9854 -0.7373 -0.7588 s.e. 0.0192 0.0798 0.0597
sigma^2 estimated as 8.377e-06: log likelihood=717.76 AIC=-1427.51 AICc=-1427.23 BIC=-1415.63
modelog1 = auto.arima(diff(diff(log(casos_ats) ,lag =7),lag=1),
stepwise = F,
method = "ML",
approximation=F,
seasonal = T,
allowdrift= F,
allowmean= F,
stationary= T,
# trace = T,
max.p = 5,
max.q = 5,
max.P = 3,
max.Q = 3)
modelog1Series: diff(diff(log(casos_ats), lag = 7), lag = 1) ARIMA(1,0,1)(0,0,1)[7] with zero mean
Coefficients: ar1 ma1 sma1 0.9854 -0.7295 -0.7548 s.e. 0.0205 0.0844 0.0609
sigma^2 estimated as 2.692e-06: log likelihood=717.82 AIC=-1427.65 AICc=-1427.36 BIC=-1415.77
Se observa que los dos procesos señalan un modelo similar, sin embargo los coeficiente presentan un pequeña diferencia. Se aplico al modelo diferenciado las mismas características que se observo en el modelo original. Notar: los valores ajustados para los datos originales no diferenciados no son equivalentes a los valores ajustados de los datos diferenciados. Para ver esto, tenga en cuenta que los valores ajustados en los datos originales están dados por \(\hat{X}_t = X_{t-1} + \phi(X_{t-1}-X_{t-2})\) mientras que para el modelo diferenciado \(\hat{Y}_t = \phi (X_{t-1}-X_{t-2})\), donde sombrero indica valor estimado, pero sucede que \(\hat{X}_t - \hat{X}_{t-1} \ne \hat{Y}_t\).
Dado que deseo estimaciones y pronostico sobre el modelo sin diferenciar dejare la tarea de diferenciar dentro del proceso ARIMA.
# Modelando continuación
modelog <- auto.arima(log(casos_ats),
stepwise = F,
method = "ML",
trace = T,
seasonal = T,
max.p = 5,
max.q = 5,
max.P = 3,
max.Q = 3,
nmodels = 200)Fitting models using approximations to speed things up…
ARIMA(0,1,0)(0,1,0)[7] : -1364.572 ARIMA(0,1,0)(0,1,1)[7] : -1376.077 ARIMA(0,1,0)(0,1,2)[7] : -1374.467 ARIMA(0,1,0)(0,1,3)[7] : -1374.812 ARIMA(0,1,0)(1,1,0)[7] : -1377.229 ARIMA(0,1,0)(1,1,1)[7] : -1375.349 ARIMA(0,1,0)(1,1,2)[7] : -1373.301 ARIMA(0,1,0)(1,1,3)[7] : -1372.67 ARIMA(0,1,0)(2,1,0)[7] : -1375.282 ARIMA(0,1,0)(2,1,1)[7] : -1373.249 ARIMA(0,1,0)(2,1,2)[7] : -1371.312 ARIMA(0,1,0)(2,1,3)[7] : Inf ARIMA(0,1,0)(3,1,0)[7] : -1373.812 ARIMA(0,1,0)(3,1,1)[7] : -1371.944 ARIMA(0,1,0)(3,1,2)[7] : Inf ARIMA(0,1,1)(0,1,0)[7] : -1374.83 ARIMA(0,1,1)(0,1,1)[7] : -1397.128 ARIMA(0,1,1)(0,1,2)[7] : -1395.452 ARIMA(0,1,1)(0,1,3)[7] : -1393.651 ARIMA(0,1,1)(1,1,0)[7] : -1396.712 ARIMA(0,1,1)(1,1,1)[7] : -1395.559 ARIMA(0,1,1)(1,1,2)[7] : -1393.533 ARIMA(0,1,1)(1,1,3)[7] : -1391.524 ARIMA(0,1,1)(2,1,0)[7] : -1395.266 ARIMA(0,1,1)(2,1,1)[7] : -1393.53 ARIMA(0,1,1)(2,1,2)[7] : -1391.396 ARIMA(0,1,1)(3,1,0)[7] : -1393.382 ARIMA(0,1,1)(3,1,1)[7] : -1391.53 ARIMA(0,1,2)(0,1,0)[7] : -1374.238 ARIMA(0,1,2)(0,1,1)[7] : -1404.467 ARIMA(0,1,2)(0,1,2)[7] : -1403.461 ARIMA(0,1,2)(0,1,3)[7] : -1401.381 ARIMA(0,1,2)(1,1,0)[7] : -1402.819 ARIMA(0,1,2)(1,1,1)[7] : -1403.576 ARIMA(0,1,2)(1,1,2)[7] : -1401.46 ARIMA(0,1,2)(2,1,0)[7] : -1402.764 ARIMA(0,1,2)(2,1,1)[7] : -1401.457 ARIMA(0,1,2)(3,1,0)[7] : -1400.924 ARIMA(0,1,3)(0,1,0)[7] : -1372.283 ARIMA(0,1,3)(0,1,1)[7] : -1405.775 ARIMA(0,1,3)(0,1,2)[7] : -1405.146 ARIMA(0,1,3)(1,1,0)[7] : -1403.483 ARIMA(0,1,3)(1,1,1)[7] : -1405.315 ARIMA(0,1,3)(2,1,0)[7] : -1403.892 ARIMA(0,1,4)(0,1,0)[7] : -1374.278 ARIMA(0,1,4)(0,1,1)[7] : -1404.033 ARIMA(0,1,4)(1,1,0)[7] : -1401.844 ARIMA(0,1,5)(0,1,0)[7] : -1377.1 ARIMA(1,1,0)(0,1,0)[7] : -1376.692 ARIMA(1,1,0)(0,1,1)[7] : -1412.476 ARIMA(1,1,0)(0,1,2)[7] : -1411.286 ARIMA(1,1,0)(0,1,3)[7] : -1409.189 ARIMA(1,1,0)(1,1,0)[7] : -1407.293 ARIMA(1,1,0)(1,1,1)[7] : -1411.34 ARIMA(1,1,0)(1,1,2)[7] : -1409.253 ARIMA(1,1,0)(1,1,3)[7] : -1407.09 ARIMA(1,1,0)(2,1,0)[7] : -1408.833 ARIMA(1,1,0)(2,1,1)[7] : -1409.241 ARIMA(1,1,0)(2,1,2)[7] : -1407.079 ARIMA(1,1,0)(3,1,0)[7] : -1407.609 ARIMA(1,1,0)(3,1,1)[7] : -1407.242 ARIMA(1,1,1)(0,1,0)[7] : -1375 ARIMA(1,1,1)(0,1,1)[7] : -1427.226 ARIMA(1,1,1)(0,1,2)[7] : -1425.676 ARIMA(1,1,1)(0,1,3)[7] : -1423.542 ARIMA(1,1,1)(1,1,0)[7] : -1412.282 ARIMA(1,1,1)(1,1,1)[7] : -1425.71 ARIMA(1,1,1)(1,1,2)[7] : -1423.976 ARIMA(1,1,1)(2,1,0)[7] : -1415.194 ARIMA(1,1,1)(2,1,1)[7] : -1423.663 ARIMA(1,1,1)(3,1,0)[7] : -1417.448 ARIMA(1,1,2)(0,1,0)[7] : -1373.449 ARIMA(1,1,2)(0,1,1)[7] : Inf ARIMA(1,1,2)(0,1,2)[7] : Inf ARIMA(1,1,2)(1,1,0)[7] : -1412.065 ARIMA(1,1,2)(1,1,1)[7] : Inf ARIMA(1,1,2)(2,1,0)[7] : -1416.103 ARIMA(1,1,3)(0,1,0)[7] : -1371.314 ARIMA(1,1,3)(0,1,1)[7] : Inf ARIMA(1,1,3)(1,1,0)[7] : -1410.607 ARIMA(1,1,4)(0,1,0)[7] : -1372.374 ARIMA(2,1,0)(0,1,0)[7] : -1374.83 ARIMA(2,1,0)(0,1,1)[7] : -1416.767 ARIMA(2,1,0)(0,1,2)[7] : -1416.113 ARIMA(2,1,0)(0,1,3)[7] : -1414.091 ARIMA(2,1,0)(1,1,0)[7] : -1409.665 ARIMA(2,1,0)(1,1,1)[7] : -1416.26 ARIMA(2,1,0)(1,1,2)[7] : -1414.343 ARIMA(2,1,0)(2,1,0)[7] : -1411.671 ARIMA(2,1,0)(2,1,1)[7] : -1414.275 ARIMA(2,1,0)(3,1,0)[7] : -1411.355 ARIMA(2,1,1)(0,1,0)[7] : Inf ARIMA(2,1,1)(0,1,1)[7] : Inf ARIMA(2,1,1)(0,1,2)[7] : Inf ARIMA(2,1,1)(1,1,0)[7] : -1412.436 ARIMA(2,1,1)(1,1,1)[7] : Inf ARIMA(2,1,1)(2,1,0)[7] : -1416.548 ARIMA(2,1,2)(0,1,0)[7] : -1371.302 ARIMA(2,1,2)(0,1,1)[7] : Inf ARIMA(2,1,2)(1,1,0)[7] : -1410.32 ARIMA(2,1,3)(0,1,0)[7] : Inf ARIMA(3,1,0)(0,1,0)[7] : -1373.119 ARIMA(3,1,0)(0,1,1)[7] : -1417.409 ARIMA(3,1,0)(0,1,2)[7] : -1416.129 ARIMA(3,1,0)(1,1,0)[7] : -1408.394 ARIMA(3,1,0)(1,1,1)[7] : -1416.212 ARIMA(3,1,0)(2,1,0)[7] : -1410.397 ARIMA(3,1,1)(0,1,0)[7] : -1371.462 ARIMA(3,1,1)(0,1,1)[7] : Inf ARIMA(3,1,1)(1,1,0)[7] : -1410.363 ARIMA(3,1,2)(0,1,0)[7] : Inf ARIMA(4,1,0)(0,1,0)[7] : -1372.228 ARIMA(4,1,0)(0,1,1)[7] : -1418.753 ARIMA(4,1,0)(1,1,0)[7] : -1407.085 ARIMA(4,1,1)(0,1,0)[7] : -1370.108 ARIMA(5,1,0)(0,1,0)[7] : -1370.259
Now re-fitting the best model(s) without approximations…
Best model: ARIMA(1,1,1)(0,1,1)[7]
modelogSeries: log(casos_ats) ARIMA(1,1,1)(0,1,1)[7]
Coefficients: ar1 ma1 sma1 0.9854 -0.7373 -0.7588 s.e. 0.0192 0.0798 0.0597
sigma^2 estimated as 8.377e-06: log likelihood=717.76 AIC=-1427.51 AICc=-1427.23 BIC=-1415.63
modelog2 = arima(diff(diff(log(casos_ats) ,lag =1),lag= 7),
order = c(1,0,1),
seasonal = list(order = c(0,0,1), period = 7 ),
include.mean = F,
method="ML")
modelog2Call: arima(x = diff(diff(log(casos_ats), lag = 1), lag = 7), order = c(1, 0, 1), seasonal = list(order = c(0, 0, 1), period = 7), include.mean = F, method = “ML”)
Coefficients: ar1 ma1 sma1 0.9854 -0.7295 -0.7548 s.e. 0.0205 0.0844 0.0609
sigma^2 estimated as 2.636e-06: log likelihood = 717.82, aic = -1427.65
modelog3 = Arima(log(casos_ats),
order = c(1,1,1),
seasonal = list(order = c(0,1,1), period = 7 ),
method="ML")
modelog3Series: log(casos_ats) ARIMA(1,1,1)(0,1,1)[7]
Coefficients: ar1 ma1 sma1 0.9854 -0.7373 -0.7588 s.e. 0.0192 0.0798 0.0597
sigma^2 estimated as 8.377e-06: log likelihood=717.76 AIC=-1427.51 AICc=-1427.23 BIC=-1415.63
modelog4 = Arima(log(casos_ats),
order = c(1,1,1),
seasonal = list(order = c(1,1,1), period = 7 ),
method="ML")
modelog4Series: log(casos_ats) ARIMA(1,1,1)(1,1,1)[7]
Coefficients: ar1 ma1 sar1 sma1 0.9831 -0.7168 -0.0882 -0.7164 s.e. 0.0229 0.0882 0.1119 0.0861
sigma^2 estimated as 8.424e-06: log likelihood=718.07 AIC=-1426.14 AICc=-1425.71 BIC=-1411.3
modelog5 = Arima(log(casos_ats),
order = c(1,1,1),
seasonal = list(order = c(2,1,1), period = 7 ),
method="ML")
modelog5Series: log(casos_ats) ARIMA(1,1,1)(2,1,1)[7]
Coefficients: ar1 ma1 sar1 sar2 sma1 0.9830 -0.7167 -0.0647 0.0435 -0.7376 s.e. 0.0232 0.0883 0.1283 0.1203 0.1008
sigma^2 estimated as 8.484e-06: log likelihood=718.14 AIC=-1424.28 AICc=-1423.66 BIC=-1406.46
Notar que se apelo a distintos modelos, guiados por el procedimiento de inspección del de las funciones de autocorrelación y el comando auto.arima.
Se concluye que el modelo que mejor ajustó a los datos, basados en los criterios AIC y SBC, corresponde a un modelo ARIMA(1,1,1)(0,1,1)[periodo estacional = 7].
4.3 PRUEBAS DE DIAGNÓSTICO
Box y Jenkins sugieren realizar test para verificar que el modelo seleccionado se ajuste correctamente a los datos. Si el modelo aproxima correctamente la serie original, los residuos deben comportarse como un ruido blanco, motivo por el cual se deberán observar las funciones de autocorrelación de estos.
# Pruebas diagnósticos
tsdisplay(residuals(modelog), lag.max=10, main='(1,1,1)x(0,1,1)_s Model Residuals')
qqnorm(modelog$residuals,main="", sub="Figura: Gráfico Q para evaluar normalidad"); qqline(modelog$residuals)
tsdiag(modelog)
checkresiduals(modelog)
Ljung-Box test
data: Residuals from ARIMA(1,1,1)(0,1,1)[7] Q* = 9.2035, df = 11, p-value = 0.6031
Model df: 3. Total lags used: 14
autoplot(modelog)
Se comenzó con el comando tsdisplay, con el cual se observa las funciones de autocorrelación simple y parcial de los residuos. Las funciones se encuentra dentro de los intervalos de confianza, mostrando así que no hay autocorrelaciones significativas en esta primera inspección.
Aplicado con qqnorm se obtuvo el diagrama Q-Q de los residuos del modelo ARIMA(1,1,1)(0,1,1)[7] estimado para la serie. Los puntos parecen seguir la línea recta y permanecer próximos a esta con lo que no puedo rechazar la hipotiposis de normalidad.
Ademas se procede a realizar un segundo test sobre los residuos, es este caso se llamara al comando tsdiag. Observando el gráfico del estadístico de Ljung-Box (Ljung, Box, 1978), todos los Q observados están fuera de la región de rechazo, por lo tanto se acepta la hipótesis nula de ruido blanco. Por su parte el test LM encontrado en el comando checkresiduals reporta resultados similares. Así también se observan que la concentración de los residuos esta en cero, con una cola larga a izquierda por el residuo aislado de gran tamaño que se observa.
Por ultimo se procede a observar gráficamente las raíces del modelo, las cuales esta cercanas al circulo unitario, indicando que este modelos pueda ser poco estable. Sin embargo dada que proporciona el mejor ajuste temporal sobre los datos se continuará con el mismo.
4.3.1 EFECTOS GARCH
Aunque el modelo de residuos aparenta comportarse como un ruido blanco, se realizara a continuación una prueba sobre los residuos al cuadrada, donde espero corroborar los resultados anteriores de no presencia de un modelo autorregresivo sobre los residuos.
# Calcular los residuales al cuadrado de nuestro modelo ARIMA (hay Heterocedastica)
Errores_cuadrado = resid(modelog)^2
plot(Errores_cuadrado) #La varianza es constante, no existe Heterocedacistidad
La varianza es constante, es decir salvo por el residuo atípico registrado en los apartados anteriores, no se observa presencia de heterocedaticidad.
A continuación se procede con una regresión de los residuos al cuadrado para terminar de “rechazar” un comportamiento GARCH sobre los datos.
# Regresión con los residuales al cuadrado rezagados
Regresion1 = dynlm(Errores_cuadrado ~ L(Errores_cuadrado))
summary(Regresion1)Time series regression with “ts” data: Start = 1(2), End = 22(5)
Call: dynlm(formula = Errores_cuadrado ~ L(Errores_cuadrado))
Residuals: Min 1Q Median 3Q Max -1.460e-05 -7.340e-06 -6.750e-06 -4.620e-06 7.298e-04
Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 7.476e-06 4.940e-06 1.513 0.132 L(Errores_cuadrado) 1.598e-02 8.186e-02 0.195 0.845
Residual standard error: 6.019e-05 on 149 degrees of freedom Multiple R-squared: 0.0002558, Adjusted R-squared: -0.006454 F-statistic: 0.03812 on 1 and 149 DF, p-value: 0.8455
Al parecer no nos enfrentamos a un modelo GARCH. Dado que no se rechaza la hipótesis nula de no efectos ARCH. Que terminaremos de corroborar en las siguientes funcione de autocorrelación
autoplot(acf(Errores_cuadrado, lag.max = 5394, ylim = c(-0.5,1))) +
labs(title = "Autocorrelación de los errores al cuadrado") +
xlab("Rezagos") +
ylab("Autocorrelación")
autoplot(pacf(Errores_cuadrado, lag.max = 5394, ylim = c(-0.5,1))) +
labs(title = "Autocorrelación parcial de los errores al cuadrado") +
xlab("Rezagos") +
ylab("Autocorrelación parcial")
Como se esperaba no hay presencia de autocorrelaciones, dado que ninguna sale de los limites de criterio, con lo que podemos decir que no hay heterocedasicidad en la varianza.
Dada la fiabilidad de los resultados diagnósticos observados se procede.
4.4 PRONÓSTICOS
Uno de los principales intereses sobre el análisis temporal, es proporcionar predicciones y pronóstico óptimos, aquella con errores cuadrático medios de predicción mínimos. Esta viene dada por la esperanza condicionada a t, es decir \(E[Y_{t+1}/Y_t]\) que se calcula de forma recursiva sobre el modelo ARIMA, estimando este modelos ARIMA como la esperanza condicionada para luego remplazar en el modelo de pronostico.
Dado que se considero a través de todas la pruebas que el modelo es el correcto para los datos, podemos especificar un horizonte de pronóstico “h” periodos por delante para realizar perdiciones sobre el mejor modelo estimado.
pronostico <- forecast(modelog, h = 30)
pronostico$x <- exp(pronostico$x)
pronostico$mean <- exp(pronostico$mean)
pronostico$lower <- exp(pronostico$lower)
pronostico$upper <- exp(pronostico$upper)
plot(pronostico,
main="Serie casos acumulado y pronóstico",
ylab="Cantidad Total",
xlab="Numero de Semana")
Los pronósticos ajustados por el modelo, muestran que la curva sigue un ascenso pronunciado con intervalos de confianza afectados por la diferenciación.
A continuación se procederá a adicionar al gráfico de pronósticos los valores reales, que no fueron tenidos en cuenta para este análisis. Para observar gráficamente que tan bueno es el modelo.
# Nueva sub-muestra de la serie original
base_salta2 <- basetf %>%
filter(., osm_admin_level_4 == "Salta" ) %>%
select(., "fecha", "nue_casosconf_diff") %>%
rename(., casos_d = nue_casosconf_diff)
base_salta2$fecha = as.Date(base_salta2$fecha, format="%d/%m/%Y")
base_salta2 <- base_salta2 %>%
mutate(.,
casos_a = cumsum(0 + casos_d),
logcasos_a = log(casos_a)) %>%
filter(., fecha<="2021-07-01" & fecha >= "2021-01-01")
casos_ats2 <- ts(base_salta2$casos_a,
start = c(1,1),
frequency = 7)
# Gráfico de comparación final
plot(pronostico,
main="Serie casos acumulado y pronóstico",
ylab="Cantidad Total",
xlab="Numero de Semana")
lines(casos_ats2, col="green")
legend("topleft",col=c("blue","green"),legend =c("Pronosticos","Serie Real"), lwd=3)
La linea azul representa los valores pronosticado, mientras la linea de color verde muestra los valores reales. Se observa que para la primera semana y media de pronósticos los valores obtenido se comportan de manera similar que los datos reales, mientras que pasado este tiempo los valores reales quedan un poco por debajo del pronósticos pero siempre dentro de las bandas de confianza proporcionadas por el modelo.
5 CONCLUSIONES
La series casos acumulados de coronavirus para la ciudad de salta durante la primera mitad del 2021, parece tener un buen ajuste de corto plazo por parte del modelo seleccionado, indicando así que la serie debería en el plazo pronosticado continuar ascendiendo.