BoxJenkins
Ana Karina Haique
2025-06-10
library(fpp3)
library(ggplot2)
library(urca)Datos
datos <- read.csv2("Table 21_1.csv", sep = ";", dec = ",")
datos<- datos |>
mutate(Tiempo = yearquarter(paste(Year, "Q", Quarter)))
datos_ts <- datos |>
select(-Year, -Quarter) |>
as_tsibble(index = Tiempo)
#logaritmos:
datos_ts <- datos_ts|>
mutate(L_DPI = log(DPI),
L_PCE = log(PCE),
L_GDP = log(GDP)) Estacionalidad
datos_ts |>
gg_season(L_GDP, labels = "right") +
labs(y = "$ (millions)",
title = "Seasonal plot")
datos_ts |>
gg_subseries(L_GDP) +
labs(
y = "$ (millions)",
title = "seasonal plot"
)
Benchmark
- Naive: el pronostico es la proyeccion de la ultima observacion
modelo <- datos_ts |>
model(Naive = NAIVE(L_GDP))
L_GDP_forcast <- modelo |> #pronostico
forecast(h = 8)
L_GDP_forcast|>
autoplot(datos_ts, level = NULL) + #Null: no IC
labs(title = "Naive")
- Seasonal Naive: proyectar el ultimo ciclo estacional
modelo <- datos_ts |>
model(Seasonal_Naive =
SNAIVE(L_GDP ~ lag("year")))
#"year" Quiero que el valor pronosticado para cada período sea igual al valor del mismo período del AÑO anterior".
L_GDP_forcast <- modelo |>
forecast(h = 8)
L_GDP_forcast|>
autoplot(datos_ts, level = NULL) +
labs(title = "Pronostico")
- Drift : proyeccion de una recta entre ultima y la primer observacion
modelo <- datos_ts |>
model(Drift = RW(L_GDP ~ drift()))
L_GDP_forcast <- modelo |>
forecast(h = 8)
L_GDP_forcast|>
autoplot(datos_ts, level = NULL) +
labs(title = "Pronostico")
- Mean
modelo <- datos_ts |>
model(Mean = MEAN(L_GDP))
L_GDP_forcast <- modelo |>
forecast(h = 8)
L_GDP_forcast|>
autoplot(datos_ts, level = NULL) +
labs(title = "Pronostico")
TODOS JUNTOS
modelos <- datos_ts |>
model(
Naive = NAIVE(L_GDP),
Seasonal_Naive = SNAIVE(L_GDP ~ lag("year")),
Drift = RW(L_GDP ~ drift()),
Mean = MEAN(L_GDP)
)
L_GDP_forcast <- modelos |>
forecast(h = 8)
L_GDP_forcast|>
autoplot(datos_ts, level = NULL) +
labs(title = "Pronostico")
###
autoplot(L_GDP_forcast, level = NULL) +
labs(title = "Pronostico")
###
recientes <- datos_ts |>
slice_tail(n = 24)
autoplot(recientes, L_GDP) +
autolayer(L_GDP_forcast, series = "fit", level = NULL)## Warning in geom_line(mapping = without(mapping, "shape"), data =
## unpack_data(object[single_row[["FALSE"]], : Ignoring unknown parameters:
## `series`
labs(title = "Pronostico")## $title
## [1] "Pronostico"
##
## attr(,"class")
## [1] "labels"Box Jenkins
Paso 1: identificar los valores de p, d, q
La serie: L_GDP
datos_ts |>
autoplot(L_GDP)
Vimos que era un proceso integrado en diferencia: de modo que integramos la serie.
datos_ts <- datos_ts |>
mutate(dL_GDP = difference(L_GDP))
datos_ts|>
autoplot(dL_GDP)## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).
## Test DFLuego, \(d=1\)
Para determinar p y q, revisamos los diagramas de ACF y PACF
datos_ts|>
ACF(dL_GDP, lag_max = 60) |>
autoplot()+
ggtitle("ACF")
Los resagos 1 y 2 son significativos de modo que podrÃa tratarse de un proceso MA(2)
datos_ts|>
ACF(dL_GDP, type = "partial", lag_max = 60) |>
autoplot()+
ggtitle("PACF")
El resago 1 son significativos de modo que podrÃa tratarse de un proceso AR(1)
datos_ts |>
gg_tsdisplay(dL_GDP, plot_type = "partial")## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).
En conclusion, Los posibles modelos son:
- MA(2)
- AR(1)
- ARMA(1,2)
⚠️️ Aclaraciones:
En una serie larga, algunos rezagos pueden aparecer significativos por pura casualidad, aunque en realidad no exista una relación estructural fuerte. Con 100 rezagos y bandas de confianza al 95%, es esperable que al menos 5% de ellos (es decir, unos 5) aparezcan falsamente como significativos simplemente por azar.
Si el rezago 12 es significativo, puede estar indicando estacionalidad anual (si fueran datos mensuales).
Lo mismo con el rezago 30, que podría estar captando un ciclo más largo. En ese caso, habría que considerar un modelo SARIMA.
Hacemos la estimacion:
mod_p0_d1_q2 <- datos_ts|>
model(ARIMA(L_GDP ~ pdq(0,1,2)))
report(mod_p0_d1_q2)## Series: L_GDP
## Model: ARIMA(0,1,2) w/ drift
##
## Coefficients:
## ma1 ma2 constant
## 0.2908 0.2010 0.0082
## s.e. 0.0626 0.0621 0.0009
##
## sigma^2 estimated as 8.501e-05: log likelihood=795.75
## AIC=-1583.5 AICc=-1583.34 BIC=-1569.53mod_p1_d1_q2 <- datos_ts|>
model(ARIMA(L_GDP ~ pdq(1,1,2)))
report(mod_p1_d1_q2)## Series: L_GDP
## Model: ARIMA(1,1,2) w/ drift
##
## Coefficients:
## ar1 ma1 ma2 constant
## 0.2016 0.1008 0.1677 0.0066
## s.e. 0.2095 0.2049 0.0809 0.0007
##
## sigma^2 estimated as 8.51e-05: log likelihood=796.12
## AIC=-1582.24 AICc=-1581.98 BIC=-1564.77mod_arima <- datos_ts|>
model(ARIMA(L_GDP))
report(mod_arima) #autoarima: ideal para arrancar## Series: L_GDP
## Model: ARIMA(2,1,2) w/ drift
##
## Coefficients:
## ar1 ar2 ma1 ma2 constant
## 1.3626 -0.7791 -1.1092 0.6208 0.0034
## s.e. 0.1476 0.1781 0.2203 0.2304 0.0003
##
## sigma^2 estimated as 8.341e-05: log likelihood=798.98
## AIC=-1585.95 AICc=-1585.6 BIC=-1564.99# o #####
modelos <- datos_ts |>
model(
p0_d1_q2 = ARIMA(L_GDP ~ pdq(0,1,2)),
p1_d1_q2 = ARIMA(L_GDP ~ pdq(1,1,2)),
arima = ARIMA(L_GDP)
)
L_GDP_forcast <- modelos |>
forecast(h = 8)
L_GDP_forcast## # A fable: 24 x 4 [1Q]
## # Key: .model [3]
## .model Tiempo L_GDP .mean
## <chr> <qtr> <dist> <dbl>
## 1 p0_d1_q2 2008 Q1 N(9.4, 8.5e-05) 9.37
## 2 p0_d1_q2 2008 Q2 N(9.4, 0.00023) 9.38
## 3 p0_d1_q2 2008 Q3 N(9.4, 0.00042) 9.39
## 4 p0_d1_q2 2008 Q4 N(9.4, 6e-04) 9.39
## 5 p0_d1_q2 2009 Q1 N(9.4, 0.00079) 9.40
## 6 p0_d1_q2 2009 Q2 N(9.4, 0.00098) 9.41
## 7 p0_d1_q2 2009 Q3 N(9.4, 0.0012) 9.42
## 8 p0_d1_q2 2009 Q4 N(9.4, 0.0014) 9.43
## 9 p1_d1_q2 2008 Q1 N(9.4, 8.5e-05) 9.37
## 10 p1_d1_q2 2008 Q2 N(9.4, 0.00023) 9.38
## # ℹ 14 more rowsrecientes <- datos_ts |>
slice_tail(n = 24)
autoplot(recientes, L_GDP) +
autolayer(L_GDP_forcast) +
labs(title = "Pronostico")
autoplot(recientes, L_GDP) +
autolayer(L_GDP_forcast, level = NULL) +
labs(title = "Pronostico")
accuracy(modelos) |>
select(.model, RMSE, MAE)## # A tibble: 3 × 3
## .model RMSE MAE
## <chr> <dbl> <dbl>
## 1 p0_d1_q2 0.00914 0.00677
## 2 p1_d1_q2 0.00913 0.00676
## 3 arima 0.00902 0.00666Pronostico en escala original: GDP
mod_p0_d1_q2 <- datos_ts|>
model(ARIMA(L_GDP ~ pdq(0,1,2)))
report(mod_p0_d1_q2)## Series: L_GDP
## Model: ARIMA(0,1,2) w/ drift
##
## Coefficients:
## ma1 ma2 constant
## 0.2908 0.2010 0.0082
## s.e. 0.0626 0.0621 0.0009
##
## sigma^2 estimated as 8.501e-05: log likelihood=795.75
## AIC=-1583.5 AICc=-1583.34 BIC=-1569.53L_GDP_forcast <- mod_p0_d1_q2 |>
forecast(h = 8)
L_GDP_forcast## # A fable: 8 x 4 [1Q]
## # Key: .model [1]
## .model Tiempo L_GDP .mean
## <chr> <qtr> <dist> <dbl>
## 1 ARIMA(L_GDP ~ pdq(0, 1, 2)) 2008 Q1 N(9.4, 8.5e-05) 9.37
## 2 ARIMA(L_GDP ~ pdq(0, 1, 2)) 2008 Q2 N(9.4, 0.00023) 9.38
## 3 ARIMA(L_GDP ~ pdq(0, 1, 2)) 2008 Q3 N(9.4, 0.00042) 9.39
## 4 ARIMA(L_GDP ~ pdq(0, 1, 2)) 2008 Q4 N(9.4, 6e-04) 9.39
## 5 ARIMA(L_GDP ~ pdq(0, 1, 2)) 2009 Q1 N(9.4, 0.00079) 9.40
## 6 ARIMA(L_GDP ~ pdq(0, 1, 2)) 2009 Q2 N(9.4, 0.00098) 9.41
## 7 ARIMA(L_GDP ~ pdq(0, 1, 2)) 2009 Q3 N(9.4, 0.0012) 9.42
## 8 ARIMA(L_GDP ~ pdq(0, 1, 2)) 2009 Q4 N(9.4, 0.0014) 9.43fc_nivel <- L_GDP_forcast |>
mutate(.mean = exp(.mean))
fc_nivel## # A fable: 8 x 4 [1Q]
## # Key: .model [1]
## .model Tiempo L_GDP .mean
## <chr> <qtr> <dist> <dbl>
## 1 ARIMA(L_GDP ~ pdq(0, 1, 2)) 2008 Q1 N(9.4, 8.5e-05) 11752.
## 2 ARIMA(L_GDP ~ pdq(0, 1, 2)) 2008 Q2 N(9.4, 0.00023) 11829.
## 3 ARIMA(L_GDP ~ pdq(0, 1, 2)) 2008 Q3 N(9.4, 0.00042) 11927.
## 4 ARIMA(L_GDP ~ pdq(0, 1, 2)) 2008 Q4 N(9.4, 6e-04) 12025.
## 5 ARIMA(L_GDP ~ pdq(0, 1, 2)) 2009 Q1 N(9.4, 0.00079) 12124.
## 6 ARIMA(L_GDP ~ pdq(0, 1, 2)) 2009 Q2 N(9.4, 0.00098) 12224.
## 7 ARIMA(L_GDP ~ pdq(0, 1, 2)) 2009 Q3 N(9.4, 0.0012) 12325.
## 8 ARIMA(L_GDP ~ pdq(0, 1, 2)) 2009 Q4 N(9.4, 0.0014) 12427.# Gráfico combinado
ggplot() +
geom_line(data = datos_ts, aes(x = Tiempo, y = GDP), color = "black") + # Serie histórica
geom_line(data = fc_nivel, aes(x = Tiempo, y = .mean), color = "blue") + # Pronóstico
labs(
title = "Serie histórica y pronóstico del PIB en niveles",
x = "Tiempo",
y = "PIB (niveles)"
)
Solo pronóstico
fc_nivel |>
ggplot(aes(x = Tiempo, y = .mean)) +
geom_line(color = "blue") +
labs(
title = "Pronóstico del PIB en niveles",
x = "Tiempo",
y = "PIB (niveles)"
)