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")
# 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")
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")
descomposicion_aditiva<-decompose(serie.pib.ts,type = "additive")
autoplot(descomposicion_aditiva,main="Descomposición Aditiva del PIB",xlab="Años/Trimestre")
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")
Tt<- ma(serie.pib.ts, 4, centre = TRUE)
autoplot(Tt,main = "Componente Tendencia del PIB [Ciclo]", xlab = "Años/Trimestres",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*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")
It<-Yt/(Tt*St)
autoplot(It,
main = "Componente Irregular del PIB",
xlab = "Años/Trimestre",
ylab = "It")
descomposicion_multiplicatica<-decompose(serie.pib.ts,type = "multiplicative")
autoplot(descomposicion_multiplicatica,main="Descomposición Multiplicativa del PIB",xlab="Años/Trimestres")
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")
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()
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()
library(TSstudio)
## Warning: package 'TSstudio' was built under R version 4.5.3
library(forecast)
ts_plot(Yt,Xtitle = "Años/Trimestre")
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()
| 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")
Yt %>%
diff(lag = 4,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(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)
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>
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