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-2022[marzo].
library(readxl)
library(tidyr)
library(dplyr)
library(forecast)
IVAE_03_22 <- read_excel("IVAE_03_22.xlsx",
col_names = FALSE, skip = 6, n_max = 10)
data.ivae<-pivot_longer(data = IVAE_03_22[1,],names_to = "vars",cols = 2:160,values_to = "indice") %>% select("indice")
serie.ivae.ts<- data.ivae %>% ts(start = c(2009,1),frequency = 12)
serie.ivae.ts %>%
autoplot(main ="IVAE ENE 2009- MAR 2022",
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-2022[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.8472467
## beta : 0
## gamma: 1
##
## Coefficients:
## [,1]
## a 118.3719784
## b 0.1600947
## s1 0.9529095
## s2 1.0192349
## s3 1.0440936
## s4 0.9962326
## s5 1.0009929
## s6 0.9838373
## s7 0.9554392
## s8 1.0145286
## s9 1.0872315
## s10 0.9748891
## s11 0.9626701
## s12 1.0059813#Generar el pronóstico:
PronosticosHW<-forecast(object = ModeloHW,h = 12,level = c(0.95))
PronosticosHW## Point Forecast Lo 95 Hi 95
## Apr 2022 112.9503 107.46492 118.4358
## May 2022 120.9752 113.57246 128.3779
## Jun 2022 124.0929 115.22236 132.9634
## Jul 2022 118.5640 108.86876 128.2592
## Aug 2022 119.2908 108.50116 130.0804
## Sep 2022 117.4038 105.81297 128.9947
## Oct 2022 114.1680 101.97016 126.3657
## Nov 2022 121.3911 107.66977 135.1125
## Dec 2022 130.2643 114.88360 145.6450
## Jan 2023 116.9603 102.34981 131.5708
## Feb 2023 115.6485 100.48405 130.8129
## Mar 2023 121.0126 83.19827 158.8270Por ejemplo en los resultados, se esperaría que el IVAE para abril 2022 sea de 107 a 118 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 2022 111.9184 105.09260 118.7442
## May 2022 119.3664 110.67455 128.0583
## Jun 2022 122.8422 112.63807 133.0464
## Jul 2022 118.4498 107.52066 129.3789
## Aug 2022 120.0189 107.93209 132.1056
## Sep 2022 118.4260 105.56932 131.2828
## Oct 2022 114.3087 101.05402 127.5633
## Nov 2022 120.2837 105.49291 135.0745
## Dec 2022 127.9776 111.38423 144.5709
## Jan 2023 114.4558 98.88074 130.0309
## Feb 2023 113.5923 97.43199 129.7527
## Mar 2023 119.3883 101.68918 137.0873#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
library(TSstudio)
library(forecast)
ts_plot(log(Yt))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.
Yt %>%
diff(lag = 12,diffences=D) %>%
diff(diffences=d) %>%
ts_cor(lag.max = 36)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.
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.
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.1022 -0.7232
## s.e. 0.1266 0.1135
##
## sigma^2 = 6.779: log likelihood = -351.22
## AIC=708.43 AICc=708.6 BIC=717.38
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.05202814 2.477766 1.68543 0.02782791 1.663587 0.4366552
## ACF1
## Training set 0.06837002modelo_estimado %>% autoplot(type="both")+theme_solarized()
modelo_estimado %>% check_res(lag.max = 36)Yt_fitted<-modelo_estimado$fitted
grafico_comparativo<-cbind(Yt,Yt_fitted)
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]}\)
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) %>% select(.model:BIC)## # A tibble: 4 × 6
## .model sigma2 log_lik AIC AICc BIC
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Orders 6.13 -347. 702. 702. 714.
## 2 Orders 6.72 -352. 707. 707. 713.
## 3 Orders 6.78 -351. 708. 709. 717.
## 4 Orders 7.99 -361. 726. 726. 731.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)