library(readxl)
library(forecast)
## Warning: package 'forecast' was built under R version 4.5.3
serie.imae <- read_excel("C:/Users/laura/Downloads/El Salvador IMAE por actividades.xlsx",
col_types = c("skip", "numeric"),
skip = 6)
serie.imae.ts <- ts(data = serie.imae,
start = c(2018, 1),
frequency = 12)
serie.imae.ts %>%
autoplot(main = "IMAE, El Salvador 2018-2026[marzo]",
xlab = "Años/Meses",
ylab = "Indice")
ma2_12 <- ma(serie.imae.ts, 12, centre = TRUE)
autoplot(serie.imae.ts,main = "IMAE, El Salvador 2018-2026[marzo]",
xlab = "Años/Meses",
ylab = "Indice")+
autolayer(ma2_12,series = "Tt")
## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_line()`).
library(magrittr)
Yt <- serie.imae.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(2018, 1), frequency = 12)
autoplot(St,
main = "Factores Estacionales",
xlab = "Años/Meses",
ylab = "Factor Estacional")
descomposicion <- decompose(serie.imae.ts, type = "additive")
It <- descomposicion$random
autoplot(It,
main = "Componente Irregular 2018-2026",
xlab = "Años/Meses",
ylab = "It")
descomposicion_aditiva<-decompose(serie.imae.ts,type = "additive")
autoplot(descomposicion_aditiva,main="Descomposición Aditiva",xlab="Años/Meses")
library(tsibble)
## Warning: package 'tsibble' was built under R version 4.5.3
##
## Adjuntando el paquete: 'tsibble'
## The following objects are masked from 'package:base':
##
## intersect, setdiff, union
library(feasts)
## Warning: package 'feasts' was built under R version 4.5.3
## Cargando paquete requerido: fabletools
## Warning: package 'fabletools' was built under R version 4.5.3
library(ggplot2)
Yt %>% as_tsibble() %>%
model(
classical_decomposition(value, type = "additive")
) %>%
components() %>%
autoplot() +
labs(title = "Descomposición Clásica Aditiva, IMAE")+xlab("Años/Meses")
## Warning: `autoplot.dcmp_ts()` was deprecated in fabletools 0.6.0.
## ℹ Please use `ggtime::autoplot.dcmp_ts()` instead.
## ℹ Graphics functions have been moved to the {ggtime} package. Please use
## `library(ggtime)` instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: Removed 24 rows containing missing values or values outside the scale range
## (`geom_line()`).
Tt<- ma(serie.imae.ts, 12, centre = TRUE)
autoplot(Tt,main = "Componente Tendencia [Ciclo]", xlab = "Años/Meses",ylab = "Tt")
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(2018, 1), frequency = 12)
autoplot(St,
main = "Factores Estacionales",
xlab = "Años/Meses",
ylab = "Factor Estacional")
It<-Yt/(Tt*St)
autoplot(It,
main = "Componente Irregular",
xlab = "Años/Meses",
ylab = "It")
descomposicion_multiplicatica<-decompose(serie.imae.ts,type = "multiplicative")
autoplot(descomposicion_multiplicatica,main="Descomposición Multiplicativa",xlab="Años/Meses")
library(tsibble)
library(feasts)
library(ggplot2)
Yt %>% as_tsibble() %>%
model(classical_decomposition(value, type = "multiplicative")) %>%
components() %>%
autoplot() +
labs(title = "Descomposición Clásica Multiplicativa, IMAE") + xlab("Años/Meses")
## Warning: Removed 24 rows containing missing values or values outside the scale range
## (`geom_line()`).
library(TSstudio)
## Warning: package 'TSstudio' was built under R version 4.5.3
ts_decompose(Yt, type = "additive", showline = TRUE)
ts_seasonal(Yt,type = "box",title = "Análisis de Valores Estacionales")
library(forecast)
#Estimar el modelo
ModeloHW<-HoltWinters(x = Yt,
seasonal = "multiplicative",
optim.start = c(0.9,0.9,0.9))
## Warning in HoltWinters(x = Yt, seasonal = "multiplicative", optim.start =
## c(0.9, : dificultades de optimización: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
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.9633383
## beta : 0
## gamma: 1
##
## Coefficients:
## [,1]
## a 132.2274689
## b 0.2152098
## s1 0.9770831
## s2 1.0264055
## s3 1.0200485
## s4 0.9932670
## s5 1.0166883
## s6 0.9794076
## s7 0.9611332
## s8 1.0135212
## s9 1.0809612
## s10 0.9696077
## s11 0.9553065
## s12 1.0058424
#Generar el pronóstico:
PronosticosHW<-forecast(object = ModeloHW,h = 12,level = c(0.95))
PronosticosHW
## Point Forecast Lo 95 Hi 95
## Apr 2026 129.4075 122.91843 135.8966
## May 2026 136.1608 126.92877 145.3928
## Jun 2026 135.5370 124.43331 146.6407
## Jul 2026 132.1922 119.69673 144.6877
## Aug 2026 135.5281 121.29717 149.7591
## Sep 2026 130.7693 115.69492 145.8436
## Oct 2026 128.5361 112.47300 144.5992
## Nov 2026 135.7603 117.71435 153.8062
## Dec 2026 145.0265 124.80033 165.2526
## Jan 2027 130.2955 111.09070 149.5002
## Feb 2027 128.5793 108.64966 148.5089
## Mar 2027 135.5976 98.76312 172.4321
#Gráfico de la serie original y del pronóstico.
PronosticosHW %>% autoplot()
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 2026 128.4312 119.6066 137.2558
## May 2026 136.2651 125.1483 147.3819
## Jun 2026 134.3411 121.8859 146.7963
## Jul 2026 132.4598 118.8600 146.0596
## Aug 2026 134.1745 119.1778 149.1712
## Sep 2026 130.7962 115.0737 146.5186
## Oct 2026 130.3138 113.6196 147.0081
## Nov 2026 136.5318 118.0220 155.0416
## Dec 2026 145.2733 124.5480 165.9986
## Jan 2027 130.3703 110.8878 149.8528
## Feb 2027 128.9350 108.8290 149.0409
## Mar 2027 135.4347 113.4679 157.4015
#Gráfico de la serie original y del pronóstico.
PronosticosHW2 %>% autoplot()
library(TSstudio)
library(forecast)
ts_plot(Yt,Xtitle = "Años/Meses")
library(kableExtra)
## Warning: package 'kableExtra' was built under R version 4.5.3
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")
Yt %>%
diff(lag = 12,diffences=D) %>%
diff(diffences=d) %>%
ts_cor(lag.max = 36)
Para elegir “p”, la parte autorregresiva regular, verificamos las barras Non-Seasonal que sobrepasan los intervalos de confianza en el gráfico PACF al inicio del mismo. Se observa que los lags no estacionales no sobrepasan claramente el intervalo de confianza en los primeros retardos, y las barras significativas se encuentran en órdenes superiores (lag 3), lo cual no constituye un patrón AR claro al inicio.
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 PACF. Se evidencia que la primera barra estacional (retardo 12) es altamente significativa, mientras que los retardos 24 y 36 no lo son.
Se concluye que el valor de “P” es 1
Para elegir “q”, la parte de media móvil regular, verificamos las barras Non-Seasonal que sobrepasan los intervalos de confianza en el gráfico ACF al inicio del mismo. Se observa que los lags no estacionales no sobrepasan el intervalo de confianza en los primeros retardos de forma consistente, y las barras apenas significativas se encuentran en órdenes superiores.
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 ACF. Se evidencia que la primera barra estacional (retardo 12) es altamente significativa, mientras que los retardos 24 y 36 no lo son.
Se concluye que el valor de “Q” es 1
Con esto se estimará el modelo SARIMA:
SARIMA(0,1,0)(1,1,1)[12]
Usando forecast
library(forecast)
library(ggthemes)
## Warning: package 'ggthemes' was built under R version 4.5.3
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.1817 -0.9979
## s.e. 0.1186 0.6130
##
## sigma^2 = 10.27: log likelihood = -231.8
## AIC=469.6 AICc=469.89 BIC=476.96
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.06265305 2.952445 2.011228 0.01413767 1.754354 0.3573143
## ACF1
## Training set 0.01975159
modelo_estimado %>% autoplot(type="both")+theme_solarized()
modelo_estimado %>% check_res(lag.max = 36)
## Warning: Ignoring 36 observations
## Warning: Ignoring 34 observations
## Warning: Ignoring 4 observations
Yt_Sarima<-modelo_estimado$fitted
Yt_HW<-PronosticosHW$fitted
grafico_comparativo<-cbind(Yt,Yt_Sarima,Yt_HW)
ts_plot(grafico_comparativo)
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)
## Warning: package 'fable' was built under R version 4.5.3
library(fabletools)
library(tidyr)
##
## Adjuntando el paquete: 'tidyr'
## The following object is masked from 'package:magrittr':
##
## extract
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.5.3
##
## Adjuntando el paquete: 'dplyr'
## The following object is masked from 'package:kableExtra':
##
## group_rows
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
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]>
## # ℹ 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_010_011 Orders 11.9 -232. 469. 469. 473. <cpl [0]> <cpl [12]>
## 2 arima_original Orders 10.3 -232. 470. 470. 477. <cpl [12]> <cpl [12]>
## 3 arima_automatico Orders 10.9 -231. 470. 471. 480. <cpl [1]> <cpl [12]>
## 4 arima_010_110 Orders 14.4 -237. 479. 479. 484. <cpl [12]> <cpl [0]>
library(forecast)
library(dplyr)
library(tsibble)
library(fable)
library(fabletools)
Yt<-Yt %>% as_tsibble() %>% rename(IMAE=value)
data.cross.validation<-Yt %>%
as_tsibble() %>%
stretch_tsibble(.init = 60,.step = 1)
TSCV<-data.cross.validation %>%
model(ARIMA(IMAE ~ pdq(0, 1, 0) + PDQ(1, 1, 1))) %>%
forecast(h=1) %>% accuracy(Yt)
## Warning in sqrt(diag(best$var.coef)): Se han producido NaNs
## Warning in sqrt(diag(best$var.coef)): Se han producido NaNs
## Warning in sqrt(diag(best$var.coef)): Se han producido NaNs
## Warning: The future dataset is incomplete, incomplete out-of-sample data will be treated as missing.
## 1 observation is missing at 2026 abr
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(IMAE ~ pdq(0, 1,… Test 0.164 3.40 2.75 0.101 2.21 0.489 0.437 -0.398