Descomposición de Series Temporales (Enfoque Tradicional)

library(readxl)
library(forecast)

#fundamental
serie.pib <- read_excel("C:/Users/Rebeca/Downloads/PIB a precios corrientes  trimestral.xls", col_types = c("skip", "numeric"),
             skip = 7)

serie.pib.ts <- ts(data = serie.pib,
                    start = c(2009, 1),
                    frequency = 4)
serie.pib.ts %>% 
  autoplot(main = "PIB, El Salvador 2009-2025[diciembre]",
                           xlab = "Años/Trimestres",
                           ylab = "Millones de USD")

Modelo Aditivo.

Metodo Manual

Componente de Tendencia Tt [Componente TCt]

# Metodo manual (Estimar el componente de Tendencia-Ciclo)
ma2_4 <- ma(serie.pib.ts, 4, centre = TRUE)
autoplot(serie.pib.ts,main = "PIB, El Salvador 2009-2025[diciembre]",
           xlab = "Años/Trimestres",
           ylab = "Millones de USD")+
  autolayer(ma2_4,series = "Tt")

Cálculo de los Factores Estacionales [Componente St]

library(magrittr)
## Warning: package 'magrittr' was built under R version 4.5.3
Yt <- serie.pib.ts #Serie original
Tt <- ma2_4 #Media móvil centrada 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) / 4
#Generar la serie de factores para cada valor de la serie original
St <-
  rep(St, len = length(Yt)) %>% ts(start = c(2009, 1), frequency = 4) 
autoplot(St,
         main = "Factores Estacionales del PIB",
         xlab = "Años/Trimestres",
         ylab = "Factor Estacional") 

#### Cálculo del Componente Irregular It

It<-Yt-Tt-St
autoplot(It,
         main = "Componente Irregular del PIB",
         xlab = "Años/Trimestres",
         ylab = "It")

Método con libreria (stats)

descomposicion_aditiva<-decompose(serie.pib.ts,type = "additive")
autoplot(descomposicion_aditiva,main="Descomposición Aditiva del PIB",xlab="Años/Trimestre")

Método con 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, PIB")+xlab("Años/Trimestre")

Método Multiplicativo

Metodo con Calculo Manual

Componente Tendencia Ciclo [Tt=TCt]

Tt<- ma(serie.pib.ts, 4, centre = TRUE)
autoplot(Tt,main = "Componente Tendencia del PIB [Ciclo]", xlab = "Años/Trimestres",ylab = "Tt")

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*4/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 = 4) 
autoplot(St,
         main = "Factores Estacionales del PIB",
         xlab = "Años/Trimestres",
         ylab = "Factor Estacional del PIB") 

Cálculo del Componente Irregular [It]

It<-Yt/(Tt*St)
autoplot(It,
         main = "Componente Irregular del PIB",
         xlab = "Años/Trimestre",
         ylab = "It")

Método con Libreria stats

descomposicion_multiplicatica<-decompose(serie.pib.ts,type = "multiplicative")
autoplot(descomposicion_multiplicatica,main="Descomposición Multiplicativa del PIB",xlab="Años/Trimestres")

Método con 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, PIB") + xlab("Años/Trimestres")

Pronóstico de series temporales, enfoque deterministico (clásico)

Pronóstico Modelo de Holt Winters

Método usando forecast y stats

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.7936087
##  beta : 0
##  gamma: 0.09650871
## 
## Coefficients:
##            [,1]
## a  9413.6758348
## b    67.1976250
## s1    0.9755610
## s2    1.0073567
## s3    0.9767203
## s4    1.0418755
#Generar el pronóstico:
PronosticosHW<-forecast(object = ModeloHW,h = 4,level = c(0.95))
PronosticosHW
##         Point Forecast    Lo 95     Hi 95
## 2025 Q4       9249.170 8866.463  9631.877
## 2026 Q1       9618.313 9075.102 10161.525
## 2026 Q2       9391.428 8745.512 10037.345
## 2026 Q3      10087.924 9404.769 10771.080
#Gráfico de la serie original y del pronóstico.
PronosticosHW %>% autoplot()

Método usando forecast (Aproximación por Espacios de los Estados ETS )

library(forecast)

#Generar el pronóstico:
PronosticosHW2<-hw(y = Yt,
                   h = 4,
                   level = c(0.95),
                   seasonal = "multiplicative",
                   initial = "optimal")
PronosticosHW2
##         Point Forecast    Lo 95     Hi 95
## 2025 Q4       9205.962 8555.651  9856.273
## 2026 Q1       9589.826 8703.576 10476.075
## 2026 Q2       9506.052 8463.353 10548.751
## 2026 Q3      10180.662 8914.629 11446.696
#Gráfico de la serie original y del pronóstico.
PronosticosHW2 %>% autoplot()

Pronóstico de series temporales, enfoque moderno (estocástico)

Metodología Box-Jenkins

library(TSstudio)
## Warning: package 'TSstudio' was built under R version 4.5.3
library(forecast)

ts_plot(Yt,Xtitle = "Años/Trimestre")

Identificación

Orden de Integración

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 para el PIB") %>% kable_material()
Ordenes de Integración para el PIB
x
d 1
D 1
#Gráfico de la serie diferenciada
Yt %>% 
  diff(lag = 4,diffences=D) %>% 
  diff(diffences=d) %>% 
  ts_plot(title = "Yt estacionaria")

Verificar los valores para (p,q) & (P,Q)

Yt %>% 
  diff(lag = 4,diffences=D) %>% 
  diff(diffences=d) %>% 
  ts_cor(lag.max = 36)

Estimación del modelo.

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(2, 1, 1))

summary(modelo_estimado)
## Series: . 
## ARIMA(0,1,0)(2,1,1)[4] 
## 
## Coefficients:
##         sar1    sar2     sma1
##       0.0964  0.0306  -0.8881
## s.e.  0.1661  0.1499   0.1697
## 
## sigma^2 = 66258:  log likelihood = -433.21
## AIC=874.42   AICc=875.12   BIC=882.92
## 
## Training set error measures:
##                    ME     RMSE      MAE        MPE     MAPE      MASE
## Training set 13.40819 241.5511 128.6044 0.08026422 1.983535 0.3522423
##                    ACF1
## Training set -0.1710487
library(forecast)
library(TSstudio)

# Ajustado a tu objeto real del PIB trimestral
pib.arima.automatico <- auto.arima(y = serie.pib.ts, 
                                   stepwise = FALSE, 
                                   approximation = FALSE)

# Mostrar el resumen con el modelo SARIMA óptimo y sus coeficientes
summary(pib.arima.automatico)
## Series: serie.pib.ts 
## ARIMA(1,0,0)(0,1,1)[4] with drift 
## 
## Coefficients:
##          ar1     sma1    drift
##       0.8163  -0.7773  72.7765
## s.e.  0.0871   0.1493  11.1742
## 
## sigma^2 = 61796:  log likelihood = -437.33
## AIC=882.66   AICc=883.35   BIC=891.23
## 
## Training set error measures:
##                     ME     RMSE      MAE        MPE     MAPE      MASE
## Training set -2.254541 235.2437 125.0171 -0.2393929 1.926357 0.3424168
##                     ACF1
## Training set -0.07401954
modelo_estimado %>% autoplot(type="both")+theme_solarized()

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

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(2,1,1)),
        arima_010_011 = ARIMA(value ~ pdq(0, 1, 0) + PDQ(1, 1, 1)),
        arima_010_110 = ARIMA(value ~ pdq(0, 1, 0) + PDQ(2, 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)(2,1,1)[4]> <ARIMA(0,1,0)(1,1,1)[4]> <ARIMA(0,1,0)(2,1,0)[4]>
## # ℹ 1 more variable: arima_automatico <model>

Validación Cruzada (Cross Validated)

library(forecast)
library(dplyr)
library(tsibble)
library(fable)
library(fabletools)

# ── PIB (Yt) ─────────────────────────────────────────────────────────────────
# Convertimos a tsibble y renombramos automáticamente la columna de datos a "PIB"
Yt_tb <- Yt %>% 
  as_tsibble() %>% 
  rename_with(~ "PIB", .cols = 2) # Cambia el nombre de la columna de datos de forma segura

data.cv.pib3 <- Yt_tb %>%
  stretch_tsibble(.init = 36, .step = 1)

TSCV_PIB3 <- data.cv.pib3 %>%
  model(
    ARIMA3   = ARIMA(PIB ~ pdq(0,1,0) + PDQ(2,1,1)),
    HoltWin3 = ETS(PIB ~ error("M") + trend("A") + season("M"))
  ) %>%
  forecast(h = 3) %>%
  accuracy(Yt_tb)

print(TSCV_PIB3)
## # A tibble: 2 × 10
##   .model   .type    ME  RMSE   MAE   MPE  MAPE  MASE RMSSE  ACF1
##   <chr>    <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ARIMA3   Test   125.  601.  351.  1.33  5.07 0.961  1.26 0.610
## 2 HoltWin3 Test   160.  741.  346.  1.96  5.11 0.947  1.55 0.656