Clase 1 - tercera unidad
JACQUELINE VANESSA CLAROS RUIZ
2022-05-31
Análisis de Series Temporales
I. Definiciones Básicas
Serie Temporal [Yt]: Conjunto de observaciones de una variable, ordenadas cronológicamente y con la misma frecuencia (m), durante el año { m=1 anual, m=2 semestral, m=3 Cuatrimestral, m=4 trimestral, m=6 Bimensual, m=12 Mensual}.
Componente de tendencia [Tt]: indica el movimiento de largo plazo de la serie, ya sea ascendente o descendente, con relativa identificabilidad en el patrón gráfico de la serie temporal.
Componente Cíclico [Ct]: recoge los movimientos derivados de la situación de auge o decadencia del nivel de la variable respecto al valor de la tendencia.
Componente Tendencia-Ciclo [TCt]: En la práctica los componentes Tt y Ct, son difíciles de separar, ya que la premisa sobre la cual se realiza la sepación, arroja resultados diferentes en función de la técnica usada. De forma práctica las rutinas de descomposición de una serie temporal, mostrará el Componente TCt que deberá interpretarse como una parte de la serie temporal de contiene los movimientos de mediano y largo plazo de la serie.
Componente Estacional [St]: este se encuentra presente, en series con alta frecuencia (m>1), y muestra los movimientos de la serie propios de la estación (como por ejemplo la navidad, u otro factor) que tenga regularidad dentro del patrón de la serie temporal.
Componente Irregular [It]: denota los movimientos de la serie que no pueden ser atribuidos al resto de componentes, se caracteriza por no tener un patrón especifico, lo que redunda en una mayor impredictibilidad de su comportamiento, normalmente esta asociado a factores o ‘shocks’ no previstos.
II. Descomposición de Series Temporales (Enfoque Tradicional)
Se usará en este apartado como datos de ejemplo, la serie del Indice de Volumen de la Actividad Económica [IVAE], para el periodo 2009-2021[marzo].
library(readxl)
library(forecast)
serie.ivae <-
read_excel("IVAE_SLV.xlsx",
col_types = c("skip", "numeric"),
skip = 5)
serie.ivae.ts <- ts(data = serie.ivae,
start = c(2009, 1),
frequency = 12)
serie.ivae.ts %>%
autoplot(main = "IVAE, El Salvador 2009-2021[marzo]",
xlab = "Años/Meses",
ylab = "Indice")
2.1. Modelo Aditivo.
Modelo Aditivo: La serie temporal es el resultado de la suma de los componentes teóricos.
2.1.1. Componente de Tendencia Tt [Componente TCt]
Se procede a estimar el componente de Tendencia-Ciclo a través de medias móviles:
ma2_12 <- ma(serie.ivae.ts, 12, centre = TRUE)
autoplot(serie.ivae.ts,main = "IVAE, El Salvador 2009-2021[marzo]",
xlab = "Años/Meses",
ylab = "Indice")+
autolayer(ma2_12,series = "Tt")
2.1.2. Cálculo de los Factores Estacionales [Componente St]
library(magrittr)
Yt <- serie.ivae.ts #Serie original
Tt <- ma2_12 #Media móvil centrada (2x12-MA) como componente de Tendencia Ciclo
SI <- Yt - Tt #Diferencia que contiene componentes Estacional e Irregular
St <- tapply(SI, cycle(SI), mean, na.rm = TRUE) #Promediando los resultados de cada mes
#Los factores estacionales deben sumar "0" en el modelo aditivo
St <- St - sum(St) / 12
#Generar la serie de factores para cada valor de la serie original
St <-
rep(St, len = length(Yt)) %>% ts(start = c(2009, 1), frequency = 12)
autoplot(St,
main = "Factores Estacionales",
xlab = "Años/Meses",
ylab = "Factor Estacional") 
2.1.3. Cálculo del Componente Irregular It.
It<-Yt-Tt-St
autoplot(It,
main = "Componente Irregular",
xlab = "Años/Meses",
ylab = "It")
2.1.4. Descomposición Aditiva (usando la libreria stats):
descomposicion_aditiva<-decompose(serie.ivae.ts,type = "additive")
autoplot(descomposicion_aditiva,main="Descomposición Aditiva",xlab="Años/Meses")
2.1.5. Descomposición Aditiva usando libreria feasts
library(tsibble)
library(feasts)
library(ggplot2)
Yt %>% as_tsibble() %>%
model(
classical_decomposition(value, type = "additive")
) %>%
components() %>%
autoplot() +
labs(title = "Descomposición Clásica Aditiva, IVAE")+xlab("Años/Meses")
2.2. Modelo Multiplicativo.
Modelo Multiplicativo: La serie temporal es el resultado de la amplificación/atenuación de la tendencia a causa del resto de los componentes.
Un modelo Multiplicativo puede expresarse como un un modelo aditivo, a través de una transformación logaritmica, (las minusculas indican el logaritmo natural de la variable/componente):
2.2.1. Componente Tendencia Ciclo [Tt=TCt]
Tt<- ma(serie.ivae.ts, 12, centre = TRUE)
autoplot(Tt,main = "Componente Tendencia [Ciclo]", xlab = "Años/Meses",ylab = "Tt")
2.2.2 Cálculo de Factores Estacionales [St]
SI<-Yt/Tt #Serie sin tendencia.
St <- tapply(SI, cycle(SI), mean, na.rm = TRUE) #Promediando los resultados de cada mes
#Los factores estacionales deben promediar "1" en el modelo multiplicativo
St <- St*12/sum(St)
#Generar la serie de factores para cada valor de la serie original
St <-
rep(St, len = length(Yt)) %>% ts(start = c(2009, 1), frequency = 12)
autoplot(St,
main = "Factores Estacionales",
xlab = "Años/Meses",
ylab = "Factor Estacional") 
2.2.3.Cálculo del Componente Irregular [It]
It<-Yt/(Tt*St)
autoplot(It,
main = "Componente Irregular",
xlab = "Años/Meses",
ylab = "It")
2.2.4 Descomposición Multiplicativa (usando la libreria stats):
descomposicion_multiplicatica<-decompose(serie.ivae.ts,type = "multiplicative")
autoplot(descomposicion_multiplicatica,main="Descomposición Multiplicativa",xlab="Años/Meses")
2.2.5 Descomposición Multiplicativa usando libreria feasts
library(tsibble)
library(feasts)
library(ggplot2)
Yt %>% as_tsibble() %>%
model(classical_decomposition(value, type = "multiplicative")) %>%
components() %>%
autoplot() +
labs(title = "Descomposición Clásica Multiplicativa, IVAE") + xlab("Años/Meses")
2.3. Descomposición usando la libreria TSstudio
library(TSstudio)
ts_decompose(Yt, type = "additive", showline = TRUE)ts_seasonal(Yt,type = "box",title = "Análisis de Valores Estacionales")III. Pronóstico de series temporales, enfoque deterministico (clásico)
3.1. Pronóstico Modelo de Holt Winters
Si la serie resulta tener un patrón estacional, está técnica permite realizar predicciones extrapolando los componentes de la serie temporal observada.
Este es el modelo más importante.
La técnica se fundamenta en la extrapolación, hay dos versiones, el holtwinter aditivo y multiplicativo. En optim.start es opcional, pero permite establecer los valores iniciales para la busqueda de los parámetros que utiliza hotwinters, uno para estimar el nivel, otro para estimar el cambio de nivel y el otro para el componente estacional.
3.1.1 Usando Stats y forecast
library(forecast)
#Estimar el modelo
ModeloHW<-HoltWinters(x = Yt,
seasonal = "multiplicative",
optim.start = c(0.9,0.9,0.9))
ModeloHW## Holt-Winters exponential smoothing with trend and multiplicative seasonal component.
##
## Call:
## HoltWinters(x = Yt, seasonal = "multiplicative", optim.start = c(0.9, 0.9, 0.9))
##
## Smoothing parameters:
## alpha: 0.8408163
## beta : 0
## gamma: 1
##
## Coefficients:
## [,1]
## a 117.0799442
## b 0.1600306
## s1 0.9502255
## s2 1.0233274
## s3 1.0518154
## s4 0.9900276
## s5 1.0007537
## s6 0.9807088
## s7 0.9600206
## s8 1.0149628
## s9 1.0915423
## s10 0.9796752
## s11 0.9584676
## s12 0.9994880#Generar el pronóstico:
PronosticosHW<-forecast(object = ModeloHW,h = 12,level = c(0.95))
PronosticosHW## Point Forecast Lo 95 Hi 95
## Apr 2021 111.4044 105.94952 116.8593
## May 2021 120.1386 112.77972 127.4976
## Jun 2021 123.6515 114.83356 132.4694
## Jul 2021 116.5461 107.01096 126.0813
## Aug 2021 117.9689 107.30373 128.6341
## Sep 2021 115.7630 104.33417 127.1918
## Oct 2021 113.4746 101.36814 125.5810
## Nov 2021 120.1312 106.57272 133.6897
## Dec 2021 129.3698 114.12877 144.6109
## Jan 2022 116.2681 101.78191 130.7543
## Feb 2022 113.9046 98.99575 128.8134
## Mar 2022 118.9394 81.61739 156.2614Por ejemplo en los resultados, se esperaría que el IVAE para abril 2021 sea de 105 a 116 y, la lectura es la misma para el resto de las observaciones.
Entre más busquemos los pronósticos del futuro más alejado, la incertidumbre crece más.
Nunca se pronóstica más allá del doble de la frecuencia, por ejemplo, si nuestra frecuencia es mensual, no deberíamos pronósticar más de dos años, se puede pero no es recomendado.
#Gráfico de la serie original y del pronóstico.
PronosticosHW %>% autoplot()
3.1.2 Usando Forecast (Aproximación por Espacios de los Estados ETS)
Forecast utiliza una versión de hotwinters, implementada en el espacio de los estados.En level se explica como un vector atómico porque se pueden escribir varios niveles de confianza. Es la misma lectura que la forma anterior.
library(forecast)
#Generar el pronóstico:
PronosticosHW2<-hw(y = Yt,
h = 12,
level = c(0.95),
seasonal = "multiplicative",
initial = "optimal")
PronosticosHW2## Point Forecast Lo 95 Hi 95
## Apr 2021 107.4928 99.59377 115.3918
## May 2021 115.3370 106.85846 123.8155
## Jun 2021 115.6946 107.18632 124.2029
## Jul 2021 109.6301 101.56399 117.6962
## Aug 2021 112.2047 103.94477 120.4646
## Sep 2021 110.2284 102.10923 118.3476
## Oct 2021 107.9393 99.98348 115.8951
## Nov 2021 114.4306 105.99022 122.8710
## Dec 2021 122.6459 113.59234 131.6994
## Jan 2022 108.5830 100.56060 116.6054
## Feb 2022 107.1239 99.20182 115.0460
## Mar 2022 113.1678 104.79016 121.5455#Gráfico de la serie original y del pronóstico.
PronosticosHW2 %>% autoplot()
IV. Pronóstico de series temporales, enfoque moderno (estocástico)
Metodología Box-Jenkins
Acá se le da un planteamiento estadístico, se trata de identificar cuál es el proceso o distribución del modelo probabilístico que generó la serie temporal. Podemos generar pronósticos.
Primero es de definir que la varianza es estable, luego se debe establecer si es estacionaria.
La estacionaria alude a esa situación donde la serie tiene media y varianza constante.
La estacional hace referencia a los cambios que puede tener la serie de acuerdo del periodo de la estación.
La presencia de un factor estacional puede coincidir en que nuestra serie no sea estacionaria. Entonces se calcula la diferencia estacional.
nsdiffs: verificar cuantas diferencias de estaciones es necesario aplicar para que nuestra diferencia sea estacionaria. Tambien debido que la serie crece puede ocurrir que la media sea diferente o vaya creciendo en el tiempo, otra fuente es la dependencia se puede poner la diferencia ordinaria y se usará para esto ndiffs.
Luego se procede al proceso autorregresivo - AR, el componente “A” hace referencia a que la serie está relacionada a valores pasados.La parte de media movil - MA hace referencia a errores de pronósticos, las diferencias grandes en el pasado que afectan en el presente.
Autorregresivos: depende de valores pasados.
Media movil: depende de errores de medición en el pasado.
library(TSstudio)
library(forecast)
ts_plot(Yt,Xtitle = "Años/Meses")Una serie se considera estacionaria, en sentido aproximado, si su media es constante a lo largo del tiempo, y su varianza es finita.
Normalmente una serie no es estacionaria, y debe ser sometida a transformaciones para lograr que lo sea. Dichas transformaciones pretenden estabilizar la media y la varianza de la serie, de tal manera que sea relativamente fácil identificar los componentes del proceso generador de los datos. > Las transformaciones para estabilizar la varianza están representadas por la familia de transformaciones de Box & Cox, la más popular es la transformación de logaritmo natural.
Las trasformaciones para estabilizar la media se consiguen a través de lo que se denomina diferencias regulares \(\Delta^d{Y_t}\) y de las diferencias estacionales \(\Delta_m^D{Y_t}\), puede ser necesario aplicar ambas d>0,D>0, o ninguna dependiendo lo que indiquen las pruebas de raices unitarias, apropiadas.
La serie se observa con poca volatilidad, a excepción del periodo de pandemia, lo cual indicaría que habrá que intervenir el modelo, para captar el efecto de la cuarentena sobre la actividad económica. Por lo que no habrá necesidad de aplicar ninguna transformación para estabilizar la varianza.
Identificación
Orden de Integración
Verificar el orden de integración ordinario y estacional (d & D) Se deberán aplicar pruebas de raices unitarias, para efectos de aplicar la técnica se usaran los comandos ndiffs & nsdiffs de la librería forecast para obtener dichos ordenes de integración.
library(kableExtra)
library(magrittr)
d<-ndiffs(Yt)
D<-nsdiffs(Yt)
ordenes_integracion<-c(d,D)
names(ordenes_integracion)<-c("d","D")
ordenes_integracion %>% kable(caption = "Ordenes de Integración") %>% kable_material()| x | |
|---|---|
| d | 1 |
| D | 1 |
#Gráfico de la serie diferenciada
Yt %>%
diff(lag = 12,diffences=D) %>%
diff(diffences=d) %>%
ts_plot(title = "Yt estacionaria")Verificar los valores para (p,q) & (P,Q)
Se construyen las funciones de autocorrelación simple y parcial.
Lleva 36 en meses, donde nos referimos a 3 años.
Yt %>%
diff(lag = 12,diffences=D) %>%
diff(diffences=d) %>%
ts_cor(lag.max = 36)- Parte ordinaria: (Azul) (p,q)
Cada una de las barras azules permite comprobar la hipótesis nula, si ese coeficiente de correlación es igual a 0 o no. En el primer ejemplo no hay valor significativo, p = 0, no hay componente autorregresivo en ACF.
En la segunda es igual, no hay valor significativo, tampoco hay componente de media movil ordinal.A nivel muestral es 0.
En esta parte ordinaria, el IVAE no depende del valor inmediato anterior ni de los errores de medición.
- Parte estacional: (Rojo) (P,Q)
Son correlaciones a un año, dos años, etc. Hay barras que salen significativamente, es significativa al primer año, al segundo ya no, al tercero tampoco. Entonces el componente autorregresivo estacional es 1, tiene mucho que ver la situación vista en el mismo mes con el año inmediato anterior. Hay un componente autorregresivo estacional.
En la media movil Q, resulta para el segundo año bastante significativa. Se puede concluir que depende de esa diferencia del error que se cometió al generar el pronóstico del año anterior. Efectivamente hay un componente de media movil estacional.
Componentes Autorregresivos
Para elegir “p”, es decir la parte autorregresiva regular, verificamos las barras “Non-Seasonal”, que sobrepasan los intervalos de confianza en el gráfico PAFC al inicio del mismo, se observa bien en los retardos 4 y 13 apenas y sobrepasa el intervalo de confianza, y se encuentran en los retardos de ordenes superiores a 3.
Se concluye entonces que el valor de “p” es 0
Para la elección de “P”, la parte autorregresiva estacional, verificamos las barras “Seasonal Lag” que sobrepasan los intervalos de confianza en el gráfico PAFC se evidencia que el la primera barra (retardo 12) es altamente significativa, y luego los retardos 24 y 36 no lo son.
Se concluye que el valor de “P” es 1
Componentes de Media Móvil
Para elegir “q”, es decir la parte de media móvil regular, verificamos las barras “Non-Seasonal”, que sobrepasan los intervalos de confianza en el gráfico AFC al inicio del mismo, se observa bien en los retardos 4 y 14 apenas y sobrepasa el intervalo de confianza, y se encuentran en los retardos de ordenes superiores a 3.
Se concluye entonces que el valor de “q” es 0
Para la elección de “Q”, la parte de media móvil estacional, verificamos las barras “Seasonal Lag” que sobrepasan los intervalos de confianza en el gráfico AFC se evidencia que el la primera barra (retardo 12) es altamente significativa, y luego los retardos 24 y 36 no lo son.
Se concluye que el valor de “Q” es 1
Con esto se estimará el modelo SARIMA:
Se concluye que el valor de “P” es 1
Estimación del modelo.
Parte ordinaria: (0,1,0) = (Componente p, diferencia regular, componente de media movil regular)
Parte estacional: (1,1,1) = (Componente autorregresivo, diferencia estacional, componente de media movil estacional)
Usando forecast
library(forecast)
library(ggthemes)
modelo_estimado <- Yt %>%
Arima(order = c(0, 1, 0),
seasonal = c(1, 1, 1))
summary(modelo_estimado)## Series: .
## ARIMA(0,1,0)(1,1,1)[12]
##
## Coefficients:
## sar1 sma1
## -0.2040 -0.6422
## s.e. 0.1586 0.1554
##
## sigma^2 = 6.897: log likelihood = -323.43
## AIC=652.85 AICc=653.04 BIC=661.55
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.04712358 2.488576 1.666188 0.02203353 1.658867 0.4981404
## ACF1
## Training set 0.08970896En la práctica nuestro modelo depende del valor del mismo mes en el año anterior y del error de medición que se cometió en el año anterior.
El MAPE es de 1.65%, quiere decir que de cada 100 valores pronosticados, la distancia o el error cometido es de 1.65%.
modelo_estimado %>% autoplot(type="both")+theme_solarized()
Se puede comprobar que el modelo es estable, o se cumple el teorema de invertibilidad si los puntos se encuentran dentro del círculo.
modelo_estimado %>% check_res(lag.max = 36)Yt_Sarima<-modelo_estimado$fitted
Yt_HW<-PronosticosHW$fitted
grafico_comparativo<-cbind(Yt,Yt_Sarima,Yt_HW)
ts_plot(grafico_comparativo)Verificación de sobre ajuste/sub ajuste
Se estimarán los modelos: Partiendo del modelo original \(SARIMA(0,1,0)(1,1,1)_{[12]}\), se estima un nuevo modelo con P-1 \(SARIMA(0,1,0)(0,1,1)_{[12]}\), y otro con Q-1: \(SARIMA(0,1,0)(1,1,0)_{[12]}\)
Nos va a permitir generar varios modelos al mismo tiempo. A través de model se pueden agrupar distintos modelos…
library(tsibble)
library(feasts)
library(fable)
library(fabletools)
library(tidyr)
library(dplyr)
a<-Yt %>% as_tsibble() %>%
model(arima_original=ARIMA(value ~ pdq(0, 1, 0) + PDQ(1, 1, 1)),
arima_010_011 = ARIMA(value ~ pdq(0, 1, 0) + PDQ(0, 1, 1)),
arima_010_110 = ARIMA(value ~ pdq(0, 1, 0) + PDQ(1, 1, 0)),
arima_automatico=ARIMA(value,ic="bic",stepwise = FALSE)
)
print(a)## # A mable: 1 x 4
## arima_original arima_010_011 arima_010_110
## <model> <model> <model>
## 1 <ARIMA(0,1,0)(1,1,1)[12]> <ARIMA(0,1,0)(0,1,1)[12]> <ARIMA(0,1,0)(1,1,0)[12]>
## # … with 1 more variable: arima_automatico <model>a %>% pivot_longer(everything(), names_to = "Model name",
values_to = "Orders") %>% glance() %>%
arrange(AICc) -> tabla
tabla## # A tibble: 4 × 9
## `Model name` .model sigma2 log_lik AIC AICc BIC ar_roots ma_roots
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <list> <list>
## 1 arima_automatico Orders 6.30 -320. 649. 649. 660. <cpl [1]> <cpl [12]>
## 2 arima_010_011 Orders 6.88 -324. 652. 653. 658. <cpl [0]> <cpl [12]>
## 3 arima_original Orders 6.90 -323. 653. 653. 662. <cpl [12]> <cpl [12]>
## 4 arima_010_110 Orders 7.60 -328. 660. 660. 666. <cpl [12]> <cpl [0]>La unidad de tiempo puede ser cualquiera, y la frecuencia, cualquier sub partición de la unidad de tiempo (por ejemplo días-horas, en cuyo caso la frecuencia sería m=24)
Validación Cruzada (Cross Validated)
Complemento del tema de simulación de los datos de corte transversal pero para series temporales. Cada línea ilustra la misma serie temporal, la data roja es la sub muestra de la serie original, se estima un modelo y se genera un pronóstrico (verde). Entonces podemos saber qué tan apropiado es nuestro modelo estimado para reproducir la información.
library(forecast)
library(dplyr)
library(tsibble)
library(fable)
library(fabletools)
Yt<-Yt %>% as_tsibble() %>% rename(IVAE=value)
data.cross.validation<-Yt %>%
as_tsibble() %>%
stretch_tsibble(.init = 60,.step = 1)
TSCV<-data.cross.validation %>%
model(ARIMA(IVAE ~ pdq(0, 1, 0) + PDQ(1, 1, 1))) %>%
forecast(h=1) %>% accuracy(Yt)
print(TSCV)## # A tibble: 1 × 10
## .model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ARIMA(IVAE ~ pdq(0, … Test 0.0234 2.98 2.08 -0.0260 2.01 0.623 0.612 0.184Podemos observar que el MAPE es de 2.01%, a diferencia del anterior que era de 1.65% en la estimación del modelo, podríamos decir entonces que nuestro modelo sigue teniendo un gran poder predictivo. Este ha sido calculado sobre la base de la diferencias de toda la data disponible con el pronóstico.