library(readxl)
library(forecast)
# Leemos el archivo completo omitiendo las 5 filas de títulos institucionales
serie.completa <- read_excel("C:/Users/Rebeca/Downloads/El Salvador IMAE por actividades.xlsx", skip = 5)
## New names:
## • `` -> `...1`
# Seleccionamos únicamente la segunda columna (donde están los valores del IMAE)
# Esto evita tener que adivinar los 'col_types' manualmente
serie.imae <- serie.completa[, 2]
# Creamos el objeto de serie de tiempo desde enero de 2018
serie.ivae.ts <- ts(data = serie.imae,
start = c(2018, 1),
frequency = 12)
# Graficamos la serie original
serie.ivae.ts %>%
autoplot(main = "IMAE, El Salvador 2018-2026[marzo]",
xlab = "Años/Meses",
ylab = "Indice")
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)
## Warning: package 'ggplot2' was built under R version 4.5.3
Yt <- serie.ivae.ts
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 every 8 hours.
## 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()`).
### Componente de Tendencia Tt
library(forecast)
ma2_12 <- ma(serie.ivae.ts, 12, centre = TRUE)
autoplot(serie.ivae.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)
## Warning: package 'magrittr' was built under R version 4.5.3
library(forecast) # Asegurar que autoplot funcione correctamente
Yt <- serie.ivae.ts # Serie original (2018 - 2026[marzo])
Tt <- ma2_12 # Media móvil centrada (2x12-MA) como tendencia
# 1. Obtener la diferencia (Modelo Aditivo)
SI <- Yt - Tt
# 2. Promediar los resultados de cada mes
St <- tapply(SI, cycle(SI), mean, na.rm = TRUE)
# 3. Ajustar para que los factores sumen exactamente "0"
St <- St - sum(St) / 12
# 4. CORRECCIÓN: Replicar los factores usando el inicio correcto (2018)
St <- rep(St, len = length(Yt)) %>%
ts(start = start(Yt), frequency = 12) # Copia automáticamente el inicio de tu serie original (2018, 1)
# 5. Graficar los factores estacionales ajustados
autoplot(St,
main = "Factores Estacionales - IMAE El Salvador",
xlab = "Años/Meses",
ylab = "Factor Estacional")
It <- Yt - Tt - St
autoplot(It,
main = "Componente Irregular - IMAE El Salvador",
xlab = "Años/Meses",
ylab = "Componente Irregular (It)")
### 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")
Tt<- ma(serie.ivae.ts, 12, centre = TRUE)
autoplot(Tt,main = "Componente Tendencia [Ciclo]", xlab = "Años/Meses",ylab = "Tt")
library(magrittr)
library(forecast)
# 1. Obtener la serie sin tendencia (Modelo Multiplicativo)
SI <- Yt / Tt
# 2. Promediar los resultados de cada mes
St <- tapply(SI, cycle(SI), mean, na.rm = TRUE)
# 3. Ajustar para que los factores promedien exactamente "1"
St <- St * 12 / sum(St)
# 4. CORRECCIÓN: Replicar los factores usando el inicio real de tu serie (2018)
St <- rep(St, len = length(Yt)) %>%
ts(start = start(Yt), frequency = 12) # Detecta automáticamente el inicio correcto
# 5. Graficar los factores estacionales multiplicativos
autoplot(St,
main = "Factores Estacionales (Modelo Multiplicativo) - IMAE",
xlab = "Años/Meses",
ylab = "Factor Estacional (Porcentaje)")
It<-Yt/(Tt*St)
autoplot(It,
main = "Componente Irregular",
xlab = "Años/Meses",
ylab = "It")
# Realizar la descomposición multiplicativa automática
descomposicion_multiplicativa <- decompose(serie.ivae.ts, type = "multiplicative")
# Graficar los 4 componentes multiplicativos
autoplot(descomposicion_multiplicativa,
main = "Descomposición Multiplicativa - IMAE El Salvador",
xlab = "Años/Meses")
library(tsibble)
library(feasts)
library(ggplot2)
library(magrittr) # Por si acaso para asegurar el operador %>%
# Descomposición clásica multiplicativa usando el enfoque moderno (Tidyverts)
Yt %>%
as_tsibble() %>%
model(classical_decomposition(value, type = "multiplicative")) %>%
components() %>%
autoplot() +
labs(title = "Descomposición Clásica Multiplicativa - IVAE El Salvador",
xlab = "Años/Meses")
## Ignoring unknown labels:
## • xlab : "Años/Meses"
## Warning: Removed 24 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Descomposición usando la libreria TSstudio
library(TSstudio)
## Warning: package 'TSstudio' was built under R version 4.5.3
ts_decompose(Yt, type = "additive", showline = TRUE)
library(TSstudio)
ts_seasonal(Yt,
type = "box",
title = "Análisis de Valores Estacionales - IVAE El Salvador")
library(forecast)
# 1. Estimar el modelo usando tus objetos definidos (Yt)
ModeloHW <- HoltWinters(x = Yt,
seasonal = "multiplicative",
optim.start = c(0.9, 0.9, 0.9))
# 2. Desplegar el objeto en consola
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.1133145
## beta : 0
## gamma: 0.5581176
##
## Coefficients:
## [,1]
## a 88.6094932
## b 0.4723193
## s1 0.7400912
## s2 1.0366972
## s3 0.9419354
## s4 0.8557054
## s5 0.8889486
## s6 0.7175855
## s7 0.7582607
## s8 0.9668398
## s9 1.0821564
## s10 0.7806829
## s11 0.6526374
## s12 0.8292333
#Generar el pronóstico:
PronosticosHW<-forecast(object = ModeloHW,h = 12,level = c(0.95))
PronosticosHW
## Point Forecast Lo 95 Hi 95
## May 2026 65.92866 45.46336 86.39397
## Jun 2026 92.84052 71.68972 113.99131
## Jul 2026 84.79910 63.48287 106.11534
## Aug 2026 77.44029 55.98874 98.89184
## Sep 2026 80.86862 58.97950 102.75775
## Oct 2026 65.61846 43.99444 87.24249
## Nov 2026 69.69608 47.57104 91.82113
## Dec 2026 89.32445 65.71526 112.93363
## Jan 2027 100.48944 75.77846 125.20042
## Feb 2027 72.86323 49.92380 95.80266
## Mar 2027 61.22066 38.76982 83.67149
## Apr 2027 78.17790 40.07261 116.28318
#Gráfico de la serie original y del pronóstico.
PronosticosHW %>% autoplot()
### Usando Forecast (Aproximación por Espacios de los Estados ETS )
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
## May 2026 70.06399 -2.680502 142.8085
## Jun 2026 103.54552 -4.487947 211.5790
## Jul 2026 97.33869 -4.705973 199.3834
## Aug 2026 85.99014 -4.580852 176.5611
## Sep 2026 75.62609 -4.395537 155.6477
## Oct 2026 64.52669 -4.058697 133.1121
## Nov 2026 61.29085 -4.143693 126.7254
## Dec 2026 90.63573 -6.548164 187.8196
## Jan 2027 119.29628 -9.164570 247.7571
## Feb 2027 64.95123 -5.282738 135.1852
## Mar 2027 55.96343 -4.800833 116.7277
## Apr 2027 77.59913 -6.997720 162.1960
#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 | 0 |
#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)
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.0386 -0.8560
## s.e. 0.1457 0.2286
##
## sigma^2 = 293: log likelihood = -376.57
## AIC=759.14 AICc=759.42 BIC=766.53
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 1.138996 15.78034 9.310496 -40.59759 73.79015 0.6028872 -0.4090379
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)
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)
##
## 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 293. -377. 757. 757. 762. <cpl [0]> <cpl [12]>
## 2 arima_original Orders 293. -377. 759. 759. 767. <cpl [12]> <cpl [12]>
## 3 arima_010_110 Orders 401. -385. 774. 774. 779. <cpl [12]> <cpl [0]>
## 4 arima_automatico Orders 223. -409. 833. 834. 851. <cpl [24]> <cpl [4]>
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)
## Warning: The future dataset is incomplete, incomplete out-of-sample data will be treated as missing.
## 1 observation is missing at 2026 may
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, 1,… Test 0.929 12.6 10.1 0.538 18.1 0.651 0.586 -0.569