Pronosticar PIB Trimestral en Honduras
Modelos Univariados
Este trabajo tiene como objetivo desarrollar un ejercicio utilizando el PIB Trimestral como referencia para evaluar diferentes modelos ARIMA y los resultados de pronóstico dentro y fuera de muestra que se obtienen, además de todas las pruebas estadísticas pertinentes. Se observa que, este tipo de ejercicios solamente son útiles para pronosticar resultados de las series en el corto plazo, puesto que los mismos presentan una reversión a la media; por tanto, en caso de requerir un análisis de mediano plazo, se recomienda el uso de modelos determinísticos.
Conviene resaltar que se programó todo el ejercicio en R, incluyendo el presente informe; esto facilita la implementación para otras series económicas de utilidad en el Banco Central de Honduras. La referencia a los códigos básicos se tomó del curso de Macroeconomic Forecasting recibido en junio de 2019 en el Study Centre Gerzensee (Suiza). Los resultados no representan ninguna postura institucional, solamente se publican con fines de mostrar una aplicación de pronóstico con modelos ARIMA.
knitr::opts_chunk$set(fig.width=6, fig.height=3.5)
options(kableExtra.auto_format = FALSE)
rm(list=ls())
library(CADFtest)
library(forecast)
library(fredr)
library(ggfortify)
library(kableExtra)
library(Metrics)
library(plotly)
library(reshape2)
library(tidyverse)
library(tsbox)
library(x12)
library(xts)
source("R/functions.R")
set.seed(123)1) Importar los datos ajustados estacionalmente
Se utiliza el PIB trimestral para Honduras y el PIB de Estados Unidos de América.
##### Definir fechas #####
inicio <- "2000-01-01"
fin <- "2019-10-01"
dates <- seq(as.Date(inicio), as.Date(fin), by = "quarter")
##### PIB Honduras, serie de tiempo #####
Data_HN <- readRDS("RData/PIB_Trimestral.rds")
Data_HN <- Data_HN %>%
dplyr::filter(Enfoque_PIB == "Gasto",
Tipo_Serie == "Desestacionalizada",
Tipo_Valor == "Constantes",
Nombre_Serie == "Producto Interno Bruto a precios de mercado") %>%
dplyr::select(Fecha, Valor)
PIB_HN = xts::xts(Data_HN[, 2], order.by = dates)
PIB_HN = ts_ts(ts_span(PIB_HN, inicio, fin))
##### PIB USA #####
# https://fred.stlouisfed.org/series/GDPC1
fredr_set_key("8828b5fd4b7cf9795f597aca7954c957")
GDP_USA = fredr::fredr(series_id = "GDPC1",
observation_start = as.Date(inicio)
)
GDP_USA = xts::xts(GDP_USA[, 3], order.by = dates)
GDP_USA = ts_ts(ts_span(GDP_USA, inicio, fin))2) Inspección de los datos: PIB_HN
if(getOutputFormat() == 'html_document') {
plot_ly(Data_HN,
x = ~Fecha,
y = ~log(Valor),
type = 'scatter', mode = 'lines') %>%
layout(title = "PIB Honduras (en logs): Valores Constantes, Serie Desestacionalizada",
xaxis = list(title = ""),
yaxis = list (title = ""))
} else {
ts_plot(
`PIB_HN`= log(PIB_HN),
title = "PIB Honduras (en logs)",
subtitle = "Valores Constantes, Serie Desestacionalizada"
)
}if(getOutputFormat() == 'html_document') {
plot_ly(Data_HN, x = ~Fecha, y = ~c(rep(NA,1),diff(log(Valor),1))*100,
type = 'scatter', mode = 'lines') %>%
layout(title = "PIB Honduras: Tasa de crecimiento trimestral",
xaxis = list(title = ""),
yaxis = list (title = "%"))
} else {
ts_plot(
`PIB_HN`= tsbox::ts_pc(PIB_HN),
title = "PIB Honduras",
subtitle = "Tasa de crecimiento respecto al trimestre anterior"
)
}if(getOutputFormat() == 'html_document') {
plot_ly(Data_HN, x = ~Fecha, y = ~c(rep(NA,4),diff(log(Valor),4))*100,
type = 'scatter', mode = 'lines') %>%
layout(title = "PIB Honduras: Tasa de crecimiento interanual",
xaxis = list(title = ""),
yaxis = list (title = "%"))
} else {
ts_plot(
`PIB_HN`= tsbox::ts_pcy(PIB_HN),
title = "PIB Honduras",
subtitle = "Tasa de crecimiento interanual"
)
}Conclusiones preliminares: podemos observar que el PIB es no estacionario. Debido a esto, estimaremos el modelo en primeras diferencias. La tasa de crecimiento trimestral (que se aproxima a través de la primera diferencia logarítmica) aparenta ser estacionaria. Pueden observarse tasas de crecimiento atípicas (outliers) en diciembre 2008 y marzo 2009, siendo mayor esta última fecha.
Regularmente, el PIB se reporta en tasas de crecimiento interanual. Al finalizar el ejercicio, se reportará el pronóstico en variaciones interanuales, calculadas a través de la estimación del modelo en primeras diferencias.
3) Pruebas de raíz unitaria
Una explicación respecto a esta prueba, basada en Pindyck & Rubinfeld (2001), se encuentra en los Anexos 1 y 2. Se realizó una prueba de Dickey Fuller Aumentada (ADF, por sus siglas en inglés)1 para verificar si las series contienen una raíz unitaria. Debido a que las series tienen una tendencia creciente, se agrega una tendencia como un componente determinístico2. Los rezagos de la regresión se seleccionan mediante el criterio de información bayesiano (BIC, por sus siglas en inglés).
Las pruebas sugieren que tenemos que utilizar primeras diferencias para obtener series estacionarias; esto puede verse con los resultados en dlogs (diferenciales logarítmicos), donde se verifica si las primeras diferencias de las series contienen una raíz unitaria. Con la serie transformada en dlogs, se rechaza la hipótesis nula de raíz unitaria. El PIB en niveles es claramente no estacionario, mientras que su tasa de crecimiento (dlogs) es estacionaria. La misma situación se presenta para el PIB de Estados Unidos de América (GDP).
# Criterio en una sola tabla
df <- list(log_PIB_HN = as.matrix(log(PIB_HN)),
dlog_PIB_HN = as.matrix(ts_diff(log(PIB_HN))),
log_GDP_USA = as.matrix(log(GDP_USA)),
dlog_GDP_USA = as.matrix(ts_diff(log(GDP_USA))))
tabla_adf(df, maxLags = 10, lagSelection = "BIC") %>%
kable_styling(latex_options = "hold_position")| transf | Nombre | Stat | Drift | None | Trend |
|---|---|---|---|---|---|
| log | GDP | est.t | -0.4468 | 3.5255 | -1.8367 |
| log | GDP | p.value | 0.8943 | 0.9998 | 0.6760 |
| log | GDP | sel.lag | 1.0000 | 1.0000 | 1.0000 |
| log | PIB | est.t | -1.1843 | 7.1234 | -2.0534 |
| log | PIB | p.value | 0.6768 | 1.0000 | 0.5619 |
| log | PIB | sel.lag | 0.0000 | 0.0000 | 0.0000 |
| dlog | GDP | est.t | -5.2652 | -3.5798 | -5.2276 |
| dlog | GDP | p.value | 0.0000 | 0.0005 | 0.0003 |
| dlog | GDP | sel.lag | 0.0000 | 0.0000 | 0.0000 |
| dlog | PIB | est.t | -8.3194 | -2.0236 | -8.4682 |
| dlog | PIB | p.value | 0.0000 | 0.0420 | 0.0000 |
| dlog | PIB | sel.lag | 0.0000 | 2.0000 | 0.0000 |
4) Modelos de pronóstico ARIMA
Debido a que la serie de PIB se incorpora en primeras diferencias, con base en las pruebas de raíz unitaria, los modelos ARIMA que se especificarán en este trabajo pueden ser definidos como un proceso ARMA(p,q)3.
\[\begin{equation} Y_t=c+\sum_{i=1}^p\phi_i Y_{t-i}+\sum_{j=0}^q\theta_j \varepsilon_{t-i} \end{equation}\]
utilizando el operador de rezagos:
\[\begin{equation} \begin{split}\phi(L) Y_{t}&=c+\theta(L)\varepsilon_t \\ \phi(L)&=1-\phi_1L-\dots-\phi_pL^p \\ \theta(L)&=\theta_0+\theta_1L+\dots+\theta_qL^q\end{split} \end{equation}\]
Los coeficientes \(\phi\) corresponden al proceso autorregresivo y \(\theta\) a los de media móvil.
a) Modelo inicial (ARIMA(1,0,0))
Nota: se modela el PIB en diferenciales logarítmicos (dlog) y se incluye una constante. En la función Arima(), el primer elemento es la serie que queremos pronosticar/modelar, el orden es un vector con el número de rezagos autorregresivos (p), el número de veces que se aplican diferenciales a la serie (d) y el número de rezagos de media móvil (q). Debido a que la serie ya está expresada en diferencias, se incluye d = 0 en todas las aplicaciones.
ARIMA(1,0,0): Evaluación de los coeficientes
dPIB_HN = ts_diff(log(PIB_HN))
Model_ar1 = forecast::Arima(dPIB_HN, order = c(1,0,0), include.constant= TRUE)
# Funcion para estadisticos de coeficientes, modelos ARIMA
table_coef(Model_ar1, caption = "Coeficientes del ARIMA(1,0,0)") %>%
kable_styling(latex_options = "hold_position")| Betas | Error_Estandar | t_Estadistic | p_value | |
|---|---|---|---|---|
| ar1 | -0.0458490 | 0.1117331 | -0.4103439 | 0.6826787 |
| intercept | 0.0095568 | 0.0011719 | 8.1551247 | 0.0000000 |
coefs_ar1 <- coefs
table_stats(Model_ar1, caption = "Estad\u00EDsticos del ARIMA(1,0,0)") %>%
kable_styling(latex_options = "hold_position")| sigma_2 | logLik | AIC | AICc | BIC |
|---|---|---|---|---|
| 0.0001205 | 245.3616 | -484.7232 | -484.4032 | -477.6149 |
stats_ar1 <- cbind(data.frame(Model = "ar1"), stats)Puede verse que solamente la constante es estadísticamente significativa.4.
ARIMA(1,0,0): Análisis de los residuales
table_resids(Model_ar1, caption = "Estad\u00EDsticos de residuales ARIMA(1,0,0): Training set error measures") %>%
kable_styling(latex_options = "hold_position")| ME | RMSE | MAE | MPE | MAPE | MASE | ACF1 |
|---|---|---|---|---|---|---|
| 1e-06 | 0.0108369 | 0.0086114 | -881.9949 | 1303.383 | 0.6860476 | -0.0028669 |
evalResids_ar1 <- cbind(data.frame(Model = "ar1"), eval_resids)
LjungBox_ar1 <- checkresiduals(Model_ar1)
Ljung-Box test
data: Residuals from ARIMA(1,0,0) with non-zero mean
Q* = 6.8691, df = 6, p-value = 0.3331
Model df: 2. Total lags used: 8autoplot(Model_ar1) +
ggtitle("Evaluaci\u00F3n del modelo ARIMA(1,0,0)") +
labs(caption = "rojo = pron\u00F3stico, negro = observado",
x = "", y = "Variaci\u00F3n trimestral") +
theme_minimal()
BoxLjung_ar1 <- Box.test(Model_ar1$residuals, lag=5, fitdf=1, type="Lj")
BoxLjung_ar1
Box-Ljung test
data: Model_ar1$residuals
X-squared = 3.087, df = 4, p-value = 0.5434jarque.bera.test(Model_ar1$residuals[!is.na(Model_ar1$residuals)])
Jarque Bera Test
data: Model_ar1$residuals[!is.na(Model_ar1$residuals)]
X-squared = 1.5448, df = 2, p-value = 0.4619Tanto la prueba de Ljung-Box como la de Box-Ljung muestran que los residuos del modelo ARIMA(1,0,0) no tienen autocorrelación. Adicionalmente, la prueba de Jarque-Bera implica normalidad en los residuos.
ARIMA(1,0,0): Pronóstico fuera de muestra
h corresponde al número de períodos a ser estimados (dos años).
For_ar1 = forecast(Model_ar1, h = 8)#, level = c(25,80,95))
df <- fortify(For_ar1)
plot_ar1 <- ggplot(data = df,aes(x = Index, y = Data*100)) +
geom_line(col = 'black') +
geom_line(aes(y=Fitted*100), col='blue', linetype = "dashed") +
geom_line(aes(y=`Point Forecast`*100, col='red')) +
ggtitle("PIB Honduras: Pron\u00F3stico ARIMA(1,0,0)") +
labs(caption = "rojo = pron\u00F3stico, azul = estimado, negro = observado",
x = "", y = "Variaci\u00F3n trimestral (%)") +
theme_minimal() +
theme(legend.position="none")
if(getOutputFormat() == 'html_document') {
ggplotly(plot_ar1) %>%
layout(annotations =
list(x = 1, y = -0.1, text = "rojo = pron\u00F3stico, azul = estimado, negro = observado,
abanico = intervalos de confianza al 80% y 95%",
showarrow = F, xref='paper', yref='paper',
xanchor='right', yanchor='auto', xshift=0, yshift=0,
font=list(size = 12, color="black"))
)
} else {
plot_ar1
}Algunas veces se necesita presentar el pronóstico en variaciones interanuales. Una manera de hacerlo es transformando las series históricas y el pronóstico de las tasas de crecimiento y graficar el resultado. El modelo se estima en dlogs, luego tenemos que transformar la serie para volverla a niveles antes de calcular las tasas de crecimiento.
# Paso 1: Unir datos historicos y de pronostico
temp1 = tsbox::ts_bind(For_ar1$x, For_ar1$mean)
temp1[1] = 0
# Paso 2: Calcular la suma acumulada de dlogs y utilizar exp para obtener el nivel
temp2 = exp(ts_cumsum(temp1))
# Paso 3: Calcular tasa de crecimiento interanual
temp3 = tsbox::ts_pcy(temp2)
# Paso 4: Separar datos historicos y de pronostico para graficar
fcstStart_ar1 = tsbox::ts_start(For_ar1$mean)
histEnd = tsbox::ts_end(For_ar1$x)
HistYoY = tsbox::ts_span(temp3, NULL, histEnd)
FcstYoY_ar1 = tsbox::ts_span(temp3, fcstStart_ar1, NULL)
if(getOutputFormat() == 'html_document') {
dfHist <- data.frame(date = as.Date(time(HistYoY)), HistYoY = as.matrix(HistYoY))
dfFor_ar1 <- data.frame(date = as.Date(time(FcstYoY_ar1)),
FcstYoY_ar1 = as.matrix(FcstYoY_ar1))
df <- merge(dfHist, dfFor_ar1, all.x = TRUE, all.y = TRUE, by="date")
plotly::plot_ly(df, x = ~date, y = ~HistYoY, name = 'Observado', type = 'scatter', mode = 'lines') %>%
add_trace(x = ~date, y = ~FcstYoY_ar1, name = 'Pron\u00F3stico ARIMA(1,0,0)', type = 'scatter', mode = 'lines') %>%
layout(title = list(text = paste0('PIB Honduras: Pron\u00F3stico ARIMA(1,0,0)')),
xaxis = list(title = ""),
yaxis = list (title = "Variaci\u00F3n interanual (%)"))
} else {
tsbox::ts_plot(
`History`= HistYoY,
`Forecast` = FcstYoY_ar1,
title = "Honduras: PIB y Pron\u00F3stico ARIMA(1,0,0)",
subtitle = "Variaci\u00F3n interanual (%)"
)
}Los resultados sugieren que este modelo sencillo puede mejorarse. Los residuales muestran valores atípicos, además que el pronóstico converge rápidamente al promedio de los datos observados y no puede incorporarse fenómenos como choques aleatorios (ej: efectos del COVID-19 o una crisis interna similar a la de 2008-2009).
En los ejercicios siguientes, se trata de seleccionar el orden de los rezagos con un criterio de información y se incluyen variables exógenas adicionales que pueden ser importantes para explicar el crecimiento del PIB en Honduras.
b) Procedimiento automático para seleccionar el modelo de pronóstico
Una estrategia para seleccionar un modelo de pronóstico es mediante el uso de la función auto.arima. Esta función selecciona automáticamente el modelo de acuerdo con algún criterio. En este ejemplo se usa el AIC.5
ARIMA automático: Evaluación de los coeficientes
Model_Auto = auto.arima(dPIB_HN,
max.p = 4, max.q = 4, d = 0, max.P = 4, max.Q = 4,
ic = c("aic"), allowdrift = T, allowmean = TRUE,
seasonal = T, stepwise = FALSE)
table_coef(Model_Auto, caption = "Coeficientes del ARIMA autom\u00E1tico") %>%
kable_styling(latex_options = "hold_position")| Betas | Error_Estandar | t_Estadistic | p_value | |
|---|---|---|---|---|
| ar1 | -0.9066397 | 0.1134672 | -7.990322 | 0.0000000 |
| ar2 | -0.1725819 | 0.1111626 | -1.552518 | 0.1245866 |
| ma1 | 0.9683044 | 0.0437542 | 22.130568 | 0.0000000 |
| intercept | 0.0095793 | 0.0010923 | 8.769501 | 0.0000000 |
coefs_Auto <- coefs
table_stats(Model_Auto, caption = "Estad\u00EDsticos del ARIMA autom\u00E1tico") %>%
kable_styling(latex_options = "hold_position")| sigma_2 | logLik | AIC | AICc | BIC |
|---|---|---|---|---|
| 0.0001093 | 249.6487 | -489.2974 | -488.4755 | -477.4502 |
stats_Auto <- cbind(data.frame(Model = "autom\u00E1tico"), stats)Todos los coeficientes (excepto el correspondiente al AR(2)) son estadísticamente significativos al 5%.
ARIMA automático: Análisis de los residuales
table_resids(Model_Auto,
caption = "Estad\u00EDsticos de residuales ARIMA autom\u00E1tico: Training set error measures") %>%
kable_styling(latex_options = "hold_position")| ME | RMSE | MAE | MPE | MAPE | MASE | ACF1 |
|---|---|---|---|---|---|---|
| -1.41e-05 | 0.0101855 | 0.0081016 | -1115.759 | 1450.788 | 0.6454329 | 0.0162983 |
evalResids_Auto <- cbind(data.frame(Model = "autom\u00E1tico"), eval_resids)
LjungBox_Auto <- checkresiduals(Model_Auto)
Ljung-Box test
data: Residuals from ARIMA(2,0,1) with non-zero mean
Q* = 3.9544, df = 4, p-value = 0.4122
Model df: 4. Total lags used: 8autoplot(Model_Auto) +
ggtitle("Evaluaci\u00F3n del modelo ARIMA autom\u00E1tico") +
labs(caption = "rojo = pron\u00F3stico, negro = observado",
x = "", y = "Variaci\u00F3n trimestral") +
theme_minimal()
BoxLjung_Auto <- Box.test(Model_Auto$residuals, lag=5, fitdf=1, type="Lj")
BoxLjung_Auto
Box-Ljung test
data: Model_Auto$residuals
X-squared = 0.95808, df = 4, p-value = 0.9161jarque.bera.test(Model_Auto$residuals[!is.na(Model_Auto$residuals)])
Jarque Bera Test
data: Model_Auto$residuals[!is.na(Model_Auto$residuals)]
X-squared = 2.6446, df = 2, p-value = 0.2665## OJO MPE y MAPE son diferentes al summaryAl igual que con el modelo ARIMA(1,0,0), los residuos no presentan autocorrelación y se ajustan a una distribución normal, con valores atípicos significativos. Nuevamente, se transforman los resultados a tasas de crecimiento interanual.
ARIMA automático: Pronóstico fuera de muestra
For_Auto = forecast(Model_Auto, h = 8)
df <- fortify(For_Auto)
plot_Auto <- ggplot(data = df,aes(x = Index, y = Data*100)) +
geom_line(col = 'black') +
geom_line(aes(y=Fitted*100), col='blue', linetype = "dashed") +
geom_line(aes(y=`Point Forecast`*100, col='red')) +
ggtitle("PIB Honduras: Pron\u00F3stico ARIMA autom\u00E1tico") +
labs(caption = "rojo = pron\u00F3stico, azul = estimado, negro = observado",
x = "", y = "Variaci\u00F3n trimestral (%)") +
theme_minimal() +
theme(legend.position="none")
if(getOutputFormat() == 'html_document') {
ggplotly(plot_Auto) %>%
layout(annotations =
list(x = 1, y = -0.1, text = "rojo = pron\u00F3stico, azul = estimado, negro = observado",
showarrow = F, xref='paper', yref='paper',
xanchor='right', yanchor='auto', xshift=0, yshift=0,
font=list(size=12, color="black"))
)
} else {
plot_Auto
}Puede observarse que el pronóstico converge con menor velocidad respecto al modelo ARIMA(1,0,0).
Variaciones interanuales:
# Paso 1: Unir datos historicos y de pronostico
temp1 = ts_bind(For_Auto$x, For_Auto$mean)
temp1[1] = 0
# Paso 2: Calcular la suma acumulada de los dlogs y aplicar exp para obtener el nivel
temp2 <- exp(ts_cumsum(temp1))
# Paso 3: Calcular la tasa de crecimiento interanual
temp3 = ts_pcy(temp2)
# Paso 4: Separar series historicas y de pronostico para el grafico
fcstStart_Auto = ts_start(For_Auto$mean)
FcstYoY_Auto = ts_span(temp3, fcstStart_Auto, NULL)
dfFor_Auto <- data.frame(date = as.Date(time(FcstYoY_Auto)),
FcstYoY_ar1 = as.matrix(FcstYoY_ar1),
FcstYoY_Auto = as.matrix(FcstYoY_Auto))
if(getOutputFormat() == 'html_document') {
df <- merge(dfHist, dfFor_Auto, all.x = TRUE, all.y = TRUE, by="date")
plot_ly(df, x = ~date, y = ~HistYoY, name = 'Observado', type = 'scatter', mode = 'lines') %>%
add_trace(x = ~date, y = ~FcstYoY_ar1, name = 'ARIMA(1,0,0)', type = 'scatter', mode = 'lines') %>%
add_trace(x = ~date, y = ~FcstYoY_Auto, name = 'ARIMA Autom\u00E1tico', type = 'scatter', mode = 'lines') %>%
layout(title = list(text = paste0('PIB Honduras: Pron\u00F3sticos ARIMA')),
xaxis = list(title = ""),
yaxis = list (title = "Variaci\u00F3n interanual (%)"))
} else {
ts_plot(
`Hist.`= ts_span(HistYoY, "2000-01-01"),
`ARIMA(1,0,0)` = FcstYoY_ar1,
`Autom.` = FcstYoY_Auto,
title = "Honduras: PIB y Pron\u00F3stico ARIMA",
subtitle = "Variaci\u00F3n interanual (%)"
)
}# Tabla de resultados
names(dfFor_Auto) <- c("Fecha","ARIMA(1,0,0)","Autom\u00E1tico")
kable(dfFor_Auto,
caption = "Comparaci\u00F3n de resultados de pron\u00F3stico: Variaci0nes Interanuales") %>%
kable_styling(c("striped", "bordered")) %>%
kable_styling(latex_options = "hold_position")| Fecha | ARIMA(1,0,0) | Automático |
|---|---|---|
| 2020-01-01 | 3.203224 | 3.367647 |
| 2020-04-01 | 4.530621 | 4.525679 |
| 2020-07-01 | 4.077984 | 4.205252 |
| 2020-10-01 | 3.888411 | 3.930016 |
| 2021-01-01 | 3.897095 | 3.832820 |
| 2021-04-01 | 3.896697 | 3.968383 |
| 2021-07-01 | 3.896715 | 3.862239 |
| 2021-10-01 | 3.896714 | 3.935063 |
Los resultados del modelo automático son bastante similares a los del modelo ARIMA(1,0,0). Puede observarse que, en todos los casos, los resultados de pronóstico fuera de muestra para un período mayor a seis meses convergen hacia valores similares a la media histórica, por tanto los mismos no son eficientes si se toman en cuenta los últimos datos observados y también si se muestra que la serie puede ser explicada por un proceso más determinístico, evaluando variables explicativas como posibles regresores.
5) Extensión del modelo
a) Control por dummies
Con el fin de que la media de los errores en el procedimiento automático no se vea afectada por valores extremos (variaciones trimestrales mayores a 1.5%), se incluyen dummy en dicho modelo para la siguiente fecha:
Primer trimestre 2009 (D09a);
Segundo trimestre 2009 (D09b);
D <- ts_xts(dPIB_HN)
D[D != 0] = 0
D09a <- D
D09a["2009-01-01"] <- 1
D09a <- ts_xts(D09a)
D09b <- D
D09b["2009-04-01"] = 1
D09b <- ts_xts(D09b)
Model_ArimaDummy = Arima(dPIB_HN, order = c(2, 0, 1),
include.constant= TRUE, xreg = ts_c(D09a, D09b))
table_coef(Model_ArimaDummy, caption = "Coeficientes del ARIMA, incluyendo dummies") %>%
kable_styling(latex_options = "hold_position")| Betas | Error_Estandar | t_Estadistic | p_value | |
|---|---|---|---|---|
| ar1 | -1.0496039 | 0.1048954 | -10.006199 | 0.0000000 |
| ar2 | -0.3696799 | 0.1036767 | -3.565700 | 0.0006236 |
| ma1 | 0.9999998 | 0.0433748 | 23.054870 | 0.0000000 |
| intercept | 0.0102271 | 0.0008426 | 12.137621 | 0.0000000 |
| D09a | -0.0355877 | 0.0067707 | -5.256154 | 0.0000012 |
| D09b | -0.0170904 | 0.0067404 | -2.535528 | 0.0132272 |
table_stats(Model_ArimaDummy, caption = "Estad\u00EDsticos del ARIMA, incluyendo dummies") %>%
kable_styling(latex_options = "hold_position")| sigma_2 | logLik | AIC | AICc | BIC |
|---|---|---|---|---|
| 8.38e-05 | 260.1834 | -506.3668 | -504.7893 | -489.7806 |
table_resids(Model_ArimaDummy,
caption = "Estad\u00EDsticos de residuales ARIMA incl. Dummies: Training set error measures") %>%
kable_styling(latex_options = "hold_position")| ME | RMSE | MAE | MPE | MAPE | MASE | ACF1 |
|---|---|---|---|---|---|---|
| -2.3e-06 | 0.0087972 | 0.0068175 | -479.8207 | 754.6297 | 0.5431329 | -0.0106213 |
LjungBox_ArimaDummy <- checkresiduals(Model_ArimaDummy)
Ljung-Box test
data: Residuals from Regression with ARIMA(2,0,1) errors
Q* = 11.897, df = 3, p-value = 0.007745
Model df: 6. Total lags used: 9autoplot(Model_ArimaDummy) + theme_minimal()
BoxLjung_ArimaDummy <- Box.test(Model_ArimaDummy$residuals, lag=5, fitdf=1, type="Lj")
BoxLjung_ArimaDummy
Box-Ljung test
data: Model_ArimaDummy$residuals
X-squared = 1.3672, df = 4, p-value = 0.8499evalResids_ArimaDummy <- cbind(data.frame(Model = "ARIMA incl. Dummies"), eval_resids)
jarque.bera.test(Model_ArimaDummy$residuals[!is.na(Model_ArimaDummy$residuals)])
Jarque Bera Test
data: Model_ArimaDummy$residuals[!is.na(Model_ArimaDummy$residuals)]
X-squared = 6.4372, df = 2, p-value = 0.04001Aunque el comportamiento de los residuos es menos aceptable tomando en cuenta los estadísticos de autocorrelación y normalidad, se ha mejorado respecto a estimación dentro de muestra al eliminar los valores atípicos más grandes.
b) Economía Estadounidense y Dummies
Adicionalmente, se incluye como variable explicativa el PIB de Estados Unidos de América (GDP) en un segundo modelo:
##### Evaluando modelo RegArima #####
dGDP_USA = ts_diff(log(GDP_USA))
dGDP_USA <- ts_xts(dGDP_USA)
dPIB_HN <- ts_xts(dPIB_HN)
Model_ArimaX = Arima(dPIB_HN, order = c(2, 0, 1), include.constant= TRUE,
xreg = ts_c(dGDP_USA, D09a, D09b))
table_coef(Model_ArimaX, caption = "Coeficientes del ARIMA, incluyendo dummies") %>%
kable_styling(latex_options = "hold_position")| Betas | Error_Estandar | t_Estadistic | p_value | |
|---|---|---|---|---|
| ar1 | -1.1338792 | 0.1030951 | -10.998376 | 0.0000000 |
| ar2 | -0.4246880 | 0.1007029 | -4.217237 | 0.0000661 |
| ma1 | 0.9999990 | 0.0518750 | 19.277079 | 0.0000000 |
| intercept | 0.0071342 | 0.0012061 | 5.914981 | 0.0000001 |
| dGDP_USA | 0.5733016 | 0.1748166 | 3.279446 | 0.0015557 |
| D09a | -0.0195547 | 0.0079605 | -2.456469 | 0.0162529 |
| D09b | -0.0154535 | 0.0062161 | -2.486047 | 0.0150550 |
table_stats(Model_ArimaX, caption = "Estad\u00EDsticos del ARIMA, incluyendo dummies") %>%
kable_styling(latex_options = "hold_position")| sigma_2 | logLik | AIC | AICc | BIC |
|---|---|---|---|---|
| 7.49e-05 | 265.1436 | -514.2872 | -512.23 | -495.3316 |
table_resids(Model_ArimaX,
caption = "Estad\u00EDsticos de residuales ARIMA incl. Regresores: Training set error measures") %>%
kable_styling(latex_options = "hold_position")| ME | RMSE | MAE | MPE | MAPE | MASE | ACF1 |
|---|---|---|---|---|---|---|
| -1.2e-05 | 0.0082602 | 0.0066943 | -713.3399 | 1016.809 | 0.5333182 | -0.0338943 |
LjungBox_ArimaX <- checkresiduals(Model_ArimaX)
Ljung-Box test
data: Residuals from Regression with ARIMA(2,0,1) errors
Q* = 17.748, df = 3, p-value = 0.0004959
Model df: 7. Total lags used: 10autoplot(Model_ArimaX) + theme_minimal()
BoxLjung_ArimaX <- Box.test(Model_ArimaX$residuals, lag=5, fitdf=1, type="Lj")
BoxLjung_ArimaX
Box-Ljung test
data: Model_ArimaX$residuals
X-squared = 2.4988, df = 4, p-value = 0.6449jarque.bera.test(Model_ArimaX$residuals[!is.na(Model_ArimaX$residuals)])
Jarque Bera Test
data: Model_ArimaX$residuals[!is.na(Model_ArimaX$residuals)]
X-squared = 1.0762, df = 2, p-value = 0.5838evalResids_ArimaX <- cbind(data.frame(Model = "ARIMA incl. Regresores"), eval_resids)6) Comparación de modelos
Prueba de Ljung-Box
Esta prueba se utiliza para verificar correlación serial en los residuales6.
\(Q^{\star}=T(T+2)\sum_{k=1}^h(T-k)^{-1}\hat\rho_k^2\)
\(T\) = número de observaciones;
\(h\) = máximo rezago a ser considerado;
\(\hat\rho_k\) es la autocorrelación para el rezago \(k\);
Valores altos de \(Q^{\star}\) sugieren que la autocorrelación no proviene de una serie que se comporta como ruido blanco. \(Q^{\star}\) tiene una distribución \(\chi^2\) con \((h-K)\) grados de libertad, donde \(K\) es el número de parámetros en el modelo. Si \(Q^{\star}\) es no significativo (p-values relativamente grandes), podemos concluir que los residuos no pueden distinguirse de un proceso ruido blanco.
LjungBox <- data.frame(Modelos = c("ARIMA(1,0,0)", "ARIMA Autom\u00E1tico",
"ARIMA con Dummies","Arima con Regresores"),
Q = c(LjungBox_ar1$statistic,LjungBox_Auto$statistic,
LjungBox_ArimaDummy$statistic,LjungBox_ArimaX$statistic),
df = c(LjungBox_ar1$parameter,LjungBox_Auto$parameter,
LjungBox_ArimaDummy$parameter,LjungBox_ArimaX$parameter),
p_value = c(LjungBox_ar1$p.value,LjungBox_Auto$p.value,
LjungBox_ArimaDummy$p.value,LjungBox_ArimaX$p.value))
kable(LjungBox,
caption = "Ljung-Box test") %>%
kable_styling(c("striped", "bordered")) %>%
kable_styling(latex_options = "hold_position")| Modelos | Q | df | p_value |
|---|---|---|---|
| ARIMA(1,0,0) | 6.869088 | 6 | 0.3331238 |
| ARIMA Automático | 3.954418 | 4 | 0.4122098 |
| ARIMA con Dummies | 11.896826 | 3 | 0.0077451 |
| Arima con Regresores | 17.747550 | 3 | 0.0004959 |
De acuerdo con estos resultados, el modelo cuyos residuales se comportan más como un proceso ruido blanco es el ARIMA seleccionado de manera automática.
Prueba de Box-Ljung
Este tests examina la hipótesis nula de independencia en una serie de tiempo dada; se trata de comprobar si un grupo de autocorrelaciones son simultáneamente cero. Esta prueba es asintóticamente equivalente a la prueba anterior; si el p-value es mayor a 0.05, los residuales son independientes y los errores pueden considerarse una aproximación a un proceso ruido blanco. Los parámetros para su cálculo son los mismos a los usados en la prueba anterior:
\(Q^{\star}=T\sum_{k=1}^K\hat\rho_k^2\)
El estadístico Q está distribuido (aproximadamente) como una \(\chi^2\) con \(K\) grados de libertad. De esta manera, si el valor calculado de \(Q^{\star}\) es mayor que el nivel crítíco (digamos, 5%), podemos asegurar con un 95% de confianza que los coeficientes de autocorrelación verdaderos \(\rho_1,...,\rho_k\) no todos son cero.
Conviene resaltar que para estas pruebas, no existe una guía clara para la elección de \(K\). Si \(K\) es pequeño, existe el peligro de no tomar en cuenta autocorrelaciones de mayor orden, pero si \(K\) es demasiado grande relativo al tamaño de la muestra, la distribución en muestras finitas está propensa a deteriorarse, divergiendo de la distribución \(\chi^2\).7
BoxLjung <- data.frame(Modelos = c("ARIMA(1,0,0)","ARIMA Autom\u00E1tico",
"ARIMA con Dummies","Arima con Regresores"),
ChiSqrt = c(BoxLjung_ar1$statistic,BoxLjung_Auto$statistic,
BoxLjung_ArimaDummy$statistic,BoxLjung_ArimaX$statistic),
df = c(BoxLjung_ar1$parameter,BoxLjung_Auto$parameter,
BoxLjung_ArimaDummy$parameter,BoxLjung_ArimaX$parameter),
p_value = c(BoxLjung_ar1$p.value,BoxLjung_Auto$p.value,
BoxLjung_ArimaDummy$p.value,BoxLjung_ArimaX$p.value))
kable(BoxLjung,
caption = "Prueba de Box-Ljung") %>%
kable_styling(c("striped", "bordered")) %>%
kable_styling(latex_options = "hold_position")| Modelos | ChiSqrt | df | p_value |
|---|---|---|---|
| ARIMA(1,0,0) | 3.0870372 | 4 | 0.5433671 |
| ARIMA Automático | 0.9580848 | 4 | 0.9160837 |
| ARIMA con Dummies | 1.3671895 | 4 | 0.8498771 |
| Arima con Regresores | 2.4987951 | 4 | 0.6448516 |
Al igual que con la prueba anterior, el modelo cuyos residuales se aproximan más a un proceso ruido blanco es el ARIMA automático.
Evaluación del pronóstico, residuales dentro de muestra
Los siguientes estadísticos fueron tomados en cuenta para evaluar si los modelos observan un buen pronóstico dentro de muestra, considerando un modelo ARIMA con un término de error igual a \(\varepsilon_t=Y_t-\hat Y_t\):
- Error promedio (ME)
\(ME=\sum_{t=1}^T\varepsilon_t\)
- Raíz del Error Cuadrático Medio (RMSE)
\(RMSE=\sqrt{\frac{\sum_{t=1}^T\varepsilon_t^2}{T}}\)
- Error Absoluto Medio (MAE)
\(MAE=\frac{\sum_{t=1}^T|\varepsilon_t|}{T}\)
- Error Porcentual Medio (MPE)
\(MPE=\frac{100\%}{T}\sum_{t=1}^T\frac{Y_t-\hat Y_t}{Y_t}\)
- Error Absoluto Porcentual Medio (MAPE)
\(MAPE=\frac{100\%}{T}\sum_{t=1}^T\bigg |\frac{Y_t-\hat Y_t}{Y_t}\bigg |\)
- Error Absoluto Escalado Medio (MASE)
\(MASE=\frac{\frac{1}{J}\sum_j|\varepsilon_j|}{\frac{1}{T-1}\sum_{t=2}^T\bigg|Y_t-Y_{t-1}\bigg|}\)
evalResids <- rbind(evalResids_ar1,evalResids_Auto,
evalResids_ArimaDummy,evalResids_ArimaX)
evalResids[,2:8] <- round(evalResids[,2:8],6)
kable(evalResids,
caption = "Evaluaci\u00F3n de Residuales: Pronóstico dentro de Muestra") %>%
kable_styling(c("striped", "bordered")) %>%
kable_styling(latex_options = "hold_position")| Model | ME | RMSE | MAE | MPE | MAPE | MASE | ACF1 |
|---|---|---|---|---|---|---|---|
| ar1 | 1.0e-06 | 0.010837 | 0.008611 | -881.9949 | 1303.3830 | 0.686048 | -0.002867 |
| automático | -1.4e-05 | 0.010185 | 0.008102 | -1115.7595 | 1450.7883 | 0.645433 | 0.016298 |
| ARIMA incl. Dummies | -2.0e-06 | 0.008797 | 0.006818 | -479.8207 | 754.6297 | 0.543133 | -0.010621 |
| ARIMA incl. Regresores | -1.2e-05 | 0.008260 | 0.006694 | -713.3399 | 1016.8091 | 0.533318 | -0.033894 |
En la tabla puede observarse que el modelo que tiene los menores valores de los estadísticos (mejor modelo) es el que incluye como regresores las dummies y GDP.
7) Resultado de pronóstico dentro de muestra
Se etima el modelo con menor RMSE usando ventanas expandibles (rolling windows) con datos pasados (primer trimestre 2004 a último trimestre 2019) y produciendo un pronóstico (“pseudo-out-of Sample”). Esta es una manera de mostrar la eficiencia del pronóstico con datos observados8.
##### Evaluando pronostico dentro de muestra, 4 trimestres
horizon = 4
AllFcst = HistYoY
startSeries = ts_start(dPIB_HN)
for(smpEnd in seq(as.Date("2004-01-01"), as.Date("2019-04-01"), by = "quarter")){
x = ts_span(dPIB_HN, NULL, as.Date(smpEnd))
xregres = ts_span(ts_c(dGDP_USA), NULL, as.Date(smpEnd))
refit = Arima(x, xreg = xregres, order = c(2, 0, 1), include.constant = T)
start.date = as.Date(smpEnd + 31 * 3)
end.date = as.Date(smpEnd + 31 * 3 * horizon)
xregres = dGDP_USA[paste(start.date,end.date,sep="::")]
fcst = forecast(refit, xreg = xregres, h = horizon)
temp1 = rbind(as.matrix(fcst$x), as.matrix(fcst$mean))
temp1b = xts(temp1, order.by = seq(startSeries, length.out= length(temp1), by = "quarter"))
temp1b[1] = 0
temp2 = exp(ts_cumsum(temp1b))
temp3 = ts_pcy(temp2)
fcst = ts_span(temp3, as.Date(smpEnd), NULL)
AllFcst = ts_c(AllFcst, fcst)
}
AllFcst = ts_span(AllFcst, "2000-01-01", "2019-10-01")
df = data.frame(as.Date(index(ts_span(ts_xts(HistYoY), "2000-01-01"))), AllFcst)
colnames(df)[1] = "Date"
meltdf <- reshape2::melt(df,id="Date")
p = ggplot(meltdf,aes(x=Date,
y=value,
colour=variable,
group=variable)) +
geom_line() +
theme_minimal() +
theme(legend.position="none") +
geom_hline(yintercept=4, linetype="solid", color = "black") +
geom_hline(yintercept=3, linetype="dashed", color = "black") +
geom_hline(yintercept=5, linetype="dashed", color = "black") +
ggtitle("Pron\u00F3stico de PIB dentro de muestra: 4 trimestres") +
theme(
axis.title.x = element_blank(),
axis.title.y = element_blank())
if(getOutputFormat() == 'html_document') {
ggplotly(p) %>%
layout(title = "Pron\u00F3stico de PIB dentro de muestra: 4 Trimestres",
xaxis = list(title = ""),
yaxis = list (title = "%"))
} else {
p
}8) Fan chart de PIB
Este ejercicio implica la construcción de un modelo de pronóstico y evaluación mediante fanchart. Adicionalmente, se calcula la probabilidad de eventos particulares utilizando simulación de Montecarlo.
8.1) Estimación de un modelo ARMA y fan chart
Se estima el modelo en diferencias logarítmicas. Alternativamente, se puede estimar en tasas de crecimiento interanual (ts_pc). Note que la función autoplot automáticamente produce un intervalo de pronóstico. Podemos especificar manualmente las probabilidades que queremos cubrir (la opción de niveles).
Con el fin de simular un impacto por el efecto del COVID-19 durante 2020, se realiza un pronóstico asumiendo comportamiento de GDP tomando en cuenta el dato del World Economic Outlook del Fondo Monetario Internacional (FMI), adicional a un impacto esperado de la economía interna similar a los dos primeros trimestres de 2009 (D09a = D09b = 1) del segundo al tercer trimestre de 2019 (asumiendo que en estos se espera el mayor impacto económico); adicionalmente, para 2021 se asume un crecimiento trimestral de GDP para 2021 de 1.175% (4.7% anual) y que la economía interna vuelva a sus niveles de producción normales (D09a = D09b = 0).
dPIB_HN <- 100*ts_diff(log(PIB_HN))
ModelTest = Arima(dPIB_HN, order = c(2, 0, 1), include.constant= TRUE,
xreg = ts_c(dGDP_USA, D09a, D09b))
inicio <- "2020-01-01"
fin <- "2021-10-01"
dates_for <- seq(as.Date(inicio), as.Date(fin), by = "quarter")
dGDP_USA <- data.frame(dates_for, dGDP_USA = c(-0.0424,-0.1,0.0575,0.0325,rep(0.01175, 4)))
dGDP_USA = xts::xts(dGDP_USA[, 2], order.by = dates_for)
D09a = data.frame(dates_for, D09a = c(0,rep(1, 2),0,rep(0, 4)))
D09b = data.frame(dates_for, D09b = c(0,rep(1, 2),0,rep(0, 4)))
xregres = ts_span(ts_c(dGDP_USA, D09a, D09b), NULL, as.Date(fin))
fcstPIB_HN = forecast(ModelTest, xreg = xregres, h = 8, level = c(50, 70, 80, 95))
IntervalRange = fcstPIB_HN$upper-fcstPIB_HN$lower
colnames(IntervalRange) = c("50%", "70%", "80%", "95%")
if(getOutputFormat() == 'html_document') {
dfHist <- data.frame(date=as.Date(time(dPIB_HN)), HistYoY=as.matrix(dPIB_HN))
dfMid <- data.frame(date=as.Date(time(fcstPIB_HN$mean)), AutoFcstYoY=as.matrix(fcstPIB_HN$mean))
dfLow <- data.frame(date=as.Date(time(fcstPIB_HN$lower)), AutoFcstYoY=as.matrix(fcstPIB_HN$lower))
dfUp <- data.frame(date=as.Date(time(fcstPIB_HN$upper)), AutoFcstYoY=as.matrix(fcstPIB_HN$upper))
df <- as.data.frame(merge(dfHist, dfMid, all.x = TRUE, all.y = TRUE, by="date"))
df <- as.data.frame(merge(df, dfLow, all.x = TRUE, all.y = TRUE, by="date"))
df <- as.data.frame(merge(df, dfUp, all.x = TRUE, all.y = TRUE, by="date"))
names(df) <- c("date","HistYoY","Mean","low50","low70","low80","low95","up50","up70","up80","up95")
plot_ly(df, x = ~date, y = ~HistYoY, name = 'Observado', type = 'scatter', mode = 'lines') %>%
add_trace(x = ~date, y = ~Mean, name = 'Mean', type = 'scatter', mode = 'lines') %>%
add_trace(x = ~date, y = ~low50, name = 'low50', type = 'scatter', mode = 'lines') %>%
add_trace(x = ~date, y = ~low70, name = 'low70', type = 'scatter', mode = 'lines') %>%
add_trace(x = ~date, y = ~up50, name = 'up50', type = 'scatter', mode = 'lines') %>%
add_trace(x = ~date, y = ~up70, name = 'up70', type = 'scatter', mode = 'lines') %>%
layout(title = list(text = paste0('PIB Honduras: Pron\u00F3stico ARIMA,
fanchart modelo con Regresores',
'<br>',
'<sup>',
'Variaci\u00F3n interanual')),
xaxis = list(title = ""),
yaxis = list (title = "%"))
} else {
autoplot(fcstPIB_HN) +
theme_minimal() +
theme(text=element_text(size=8))
}kable(fcstPIB_HN)| Point Forecast | Lo 50 | Hi 50 | Lo 70 | Hi 70 | Lo 80 | Hi 80 | Lo 95 | Hi 95 | |
|---|---|---|---|---|---|---|---|---|---|
| 2020 Q1 | -1.6779882 | -2.2651336 | -1.090843 | -2.5802066 | -0.7757698 | -2.7935828 | -0.5623936 | -3.3841428 | 0.0281665 |
| 2020 Q2 | -8.6172060 | -9.2105287 | -8.023883 | -9.5289166 | -7.7054954 | -9.7445376 | -7.4898744 | -10.3413109 | -6.8931011 |
| 2020 Q3 | 0.6024155 | -0.0100639 | 1.214895 | -0.3387316 | 1.5435626 | -0.5613144 | 1.7661454 | -1.1773559 | 2.3821869 |
| 2020 Q4 | 2.5121833 | 1.8650673 | 3.159299 | 1.5178128 | 3.5065538 | 1.2826426 | 3.7417241 | 0.6317631 | 4.3926036 |
| 2021 Q1 | 1.4205776 | 0.7509245 | 2.090231 | 0.3915763 | 2.4495790 | 0.1482158 | 2.6929395 | -0.5253319 | 3.3664871 |
| 2021 Q2 | 1.3763950 | 0.6983074 | 2.054483 | 0.3344331 | 2.4183569 | 0.0880074 | 2.6647826 | -0.5940238 | 3.3468138 |
| 2021 Q3 | 1.3848727 | 0.7051075 | 2.064638 | 0.3403329 | 2.4294124 | 0.0932976 | 2.6764478 | -0.5904210 | 3.3601663 |
| 2021 Q4 | 1.3940239 | 0.7142015 | 2.073846 | 0.3493962 | 2.4386516 | 0.1023401 | 2.6857077 | -0.5814360 | 3.3694838 |
Interpretación: El gráfico muestra la media de pronóstico (línea sólida) y el intervalo de pronóstico. Los números asociados con el intervalo de pronóstico representan la cobertura. Los intervalos representan una medida de la confianza que tenemos en el pronóstico: si son muy grandes, el pronóstico puntual no es muy exacto y vice versa.
Aunque es difícil de observar en muestras cortas, el intervalo de pronóstico se agranda para horizontes más largos (debido a que la varianza del error de pronóstico también se incrementa).
8.2) Calcular la varianza del error de pronóstico VEP y simular la densidad de pronóstico
a) Extraer la varianza del error de pronóstico y comparar con la varianza incondicional
Esta es una función útil para calcular la varianza del error de pronóstico asumiendo que dichos errores están normalmente distribuidos.
sigh2 = getForecastVariance(fcstPIB_HN)
autoplot(sigh2) +
theme_minimal() +
geom_hline(yintercept=var(dPIB_HN, na.rm = T), linetype="dashed") +
ggtitle("Varianza del error de pron\u00F3stico y varianza incondicional")
b) Simular la densidad de pronóstico y comparar con el intervalo de pronóstico analítico obtenido
Se realiza una simulación de la densidad de pronóstico asumiendo que el error de pronóstico está distribuido como una normal. Esto implica que \(y_{t+h} \thicksim N(y_{t+h|t}, sigh^2)\), esto es, el valor futuro del crecimiento del PIB está distribuido como una normal con media igual al pronóstico puntual y varianza igual a la varianza del error de pronóstico.
NSim = 100
fcsth = fcstPIB_HN$mean
SimFcst = matrix(NA, nrow = 8, ncol = NSim)
for (h in 1:8){
SimFcst[h, ] = rnorm(NSim, fcsth[h], sqrt(sigh2[h]))
}
SimFcst = xts(SimFcst, order.by = as.Date(index(ts_xts(fcsth))))Podemos graficar las simulaciones. Cada línea es un posible resultado. El promedio será igual a la media de pronóstico y podemos calcular percentiles que sean iguales a los intervalos de pronóstico.
plot(SimFcst)
ci95 =t(apply(SimFcst, 1, quantile, probs = c(0.025, 0.975), na.rm = TRUE) )
ci95 = xts(ci95, order.by = as.Date(index(ts_xts(fcsth))))allci = ts_c(ci95,fcstPIB_HN$lower[, "95%"], fcstPIB_HN$upper[, "95%"])
plot(allci)
c) Calcular probabilidad de crecimiento dentro de un rango del PIB en Q2 2020
La ventaja de sumular la incertidumbre del pronóstico es que podemos estimar la probabilidad de eventos particulares.
# Probabilidad de crecimiento menor o igual a -2% en el primer trimestre de 2020
(paste0("Probabilidad de crecimiento menor a -2% en 2020q1: ",PNeg2020 = mean(SimFcst["2020-01-01",] <= -2)))[1] "Probabilidad de crecimiento menor a -2% en 2020q1: 0.32"# Probabilidad de crecimiento del PIB menor a -4% in Q1 2020
(paste0("Probabilidad de crecimiento menor a -4% en 2020q1: ",PLarge2020 = mean(SimFcst["2020-01-01",] <= -4)))[1] "Probabilidad de crecimiento menor a -4% en 2020q1: 0"# Probailidad de crecimiento del PIB entre -2% y -4% in Q1 2020
(paste0("Probabilidad de crecimiento entre -4% y -2% en 2020q1: ",
PBetween2020 = mean(SimFcst["2020-01-01",] <= -2 & SimFcst["2020-01-01",] >= -4)))[1] "Probabilidad de crecimiento entre -4% y -2% en 2020q1: 0.32"# Graficar barras
PNeg = ts_rowMeans(SimFcst <= -2)
PLarge = ts_rowMeans(SimFcst <= -4)
PBetween = ts_rowMeans(SimFcst <= -2 & SimFcst >= -4)
prob.data = data.frame(index(PNeg), PNeg, PLarge, PBetween)
colnames(prob.data) = c("Date", "<= -2%", "<= -4%", "entre -2% y -4%")
prob.data <- melt(prob.data,id.vars = "Date")
names(prob.data) <- c("Date","Rangos","Valor")
p <- ggplot(prob.data, aes(x = Date, y = Valor, fill=Rangos)) +
geom_bar(stat='identity', position=position_dodge()) +
theme_minimal() +
theme(axis.text.x = element_blank()) +
ggtitle("Probabilidades de Crecimientos Trimestrales por Rangos")
if(getOutputFormat() == 'html_document') {
ggplotly(p)
} else {
p
}d) Probabilidad que el crecimiento anual del PIB en 2020 sea menor a -4%
Paso 1: Calcular el nivel para cada simulación y agregarlo a una frecuencia anual;
# Paso 1: Calcular el nivel para cada simulacion y agregarlo a una frecuencia anual
for (s in 1:NSim){
# Agregar datos historicos al pronostico
temp = ts_bind(dPIB_HN, SimFcst[,s])
temp[1] = 0
# Cumulative sum of log-differences divided by 100
temp = ts_cumsum(temp)#/100)
# Calculate exponential of series
temp = exp(temp/100)
# Aggregate to annual frequency
temp = ts_frequency(temp, to ="year", aggregate="sum")
# Calculate growth rates of annual series
temp = ts_pc(temp)
# Save result in collection of time series
if (s == 1){
SimFcstGrt = temp
} else {
SimFcstGrt = ts_c(SimFcstGrt, temp)
}
}Paso 2: Calcular la probabilidad;
# Step 2: Calculate the probability
p2020 = mean(SimFcstGrt["2020-01-01",] <= -4)
p2021 = mean(SimFcstGrt["2021-01-01",] <= -4)
# Plot some interval forecasts for annual CPI growth
ci95 = as.xts(t(apply(SimFcstGrt, 1, quantile, probs = c(0.025, .5, 0.975), na.rm = TRUE) ))
df <- data.frame(date=as.Date(time(ci95)), HistYoY=as.matrix(ci95))
rownames(df) <- c()
names(df) <- c("date","Min","Mediana","Max")
if(getOutputFormat() == 'html_document') {
plot_ly(df, x = ~date, y = ~Min, name = 'Min', type = 'scatter', mode = 'lines') %>%
add_trace(x = ~date, y = ~Mediana, name = 'Mediana', type = 'scatter', mode = 'lines', fill = 'tonexty') %>%
add_trace(x = ~date, y = ~Max, name = 'Max', type = 'scatter', mode = 'lines', fill = 'tonexty') %>%
layout(title = list(text = paste0('PIB Honduras: Pron\u00F3stico ARIMA, Fanchart modelo RegArima',
'<br>',
'<sup>',
'Variaci\u00F3n anual')),
xaxis = list(title = ""),
yaxis = list (title = "%"))
} else {
ts_plot(
`Superior` = ci95[,"97.5%"],
`Mediana` = ci95[,"50%"],
`Inferior` = ci95[,"2.5%"],
`PIB`= SimFcstGrt["1981-01-01/2021-12-01",1],
title = "PIB Honduras",
subtitle = "Tasa de crecimiento anual"
)
}Los resultados indican que la probabilidad que la tasa de crecimiento anual del PIB sea menor a -4% a finales de 2020 es 94% y para 2021 de 0%.
Anexos
Los Anexos 1 y 2 se basan en Pyndick & Rubinfeld (2001), capítulo 16, pg.515 en adelante.
Anexo 1: Caminata aleatoria
El ejemplo más simple de una serie de tiempo estocástica es el proceso de caminata aleatoria; en este proceso, cada cambio sucesivo en \(y_t\) es extraído en forma independiente de una distribución de probabilidad con media cero. Por tanto, \(y_t\) está determinado por
\[\begin{equation} \tag{A1.1} \label{eqn:A1.1} Y_t=Y_{t-1}+\varepsilon_t \end{equation}\]
con \(E(\varepsilon_t)=0\) y \(E(\varepsilon_t,\varepsilon_s)=0\) para \(t\neq s\)
Consideremos el caso en el que se requiere un pronóstico para un proceso de caminata aleatoria. El pronóstico está dado por
\[\begin{equation} \tag{A1.2} \label{eqn:A1.2} \hat Y_{T+1}=E(Y_{T+1}|Y_T,\dots,Y_1)=Y_T+E(\varepsilon_{T+1})=Y_T \end{equation}\]
El pronóstico dos períodos adelante es:
\[\begin{equation} \tag{A1.3} \label{eqn:A1.3} \hat Y_{T+2}=E(Y_{T+2}|Y_T,\dots,Y_1)=E(Y_{T+1}+\varepsilon_{T+2})=E(Y_{T}+\varepsilon_{T+1}+\varepsilon_{T+2})=Y_t \end{equation}\]
Del mismo modo, el pronóstico \(l\) periodos adelante también es \(Y_T\).
Aunque el pronóstico \(\hat Y_{T+1}\) será el mismo sin importar cuán grande sea \(l\), la varianza del error de pronóstico crecerá conforme \(l\) se haga mayor. Para un periodo el error de pronóstico está dado por:
\[\begin{equation} \tag{A1.4} \label{eqn:A1.4} e_1=Y_{T+1}-\hat Y_{T+1}=Y_T+\varepsilon_{T+1}-Y_t=\varepsilon_{T+1} \end{equation}\]
y su varianza es \(E(\varepsilon^2_{T+1})=\sigma^2_{\varepsilon}\). Para el pronóstico de dos periodos
\[\begin{equation} \tag{A1.5} \label{eqn:A1.5} e_2=Y_{T+2}-\hat Y_{T+2}=Y_T+\varepsilon_{T+1}+\varepsilon_{T+2}-Y_t=\varepsilon_{T+1}+\varepsilon_{T+2} \end{equation}\]
y su varianza es:
\[\begin{equation} \tag{A1.6} \label{eqn:A1.6} E[(\varepsilon_{T+1}+\varepsilon_{T+2})^2]=E(\varepsilon_{T+1}^2)+E(\varepsilon_{T+2}^2)+2E(\varepsilon_{T+1}\varepsilon_{T+2}) \end{equation}\]
Dado que \(\varepsilon_{T+1}\) y \(\varepsilon_{T+2}\) son independientes, el tercer término de la ecuación \(\eqref{eqn:A1.6}\) es cero y la varianza del error es \(2\sigma_{\varepsilon}^2\). De igual forma, para el pronóstico del periodo \(l\), la varianza del error es \(l\sigma^2_{\varepsilon}\). Así, el error estándar del pronóstico se incrementa con la raíz cuadrada de \(l\). Por tanto se pueden obtener intervalos de confianza para nuestros pronósticos y estos intervalos se volverán más amplios conforme se incremente el horizonte del pronóstico.
Anexo 2: Criterio ADF
Suponga que creemos que una variable \(Y_t\), la cual ha estado creciendo con el tiempo, puede describirse por la siguiente ecuación:
\[\begin{equation} \tag{A2.1} \label{eqn:A2.1} Y_t=\alpha+\beta\cdot t+\rho\cdot Y_{t-1}+\varepsilon_t \end{equation}\]
Una posibilidad es que \(Y_t\) ha estado creciendo debido a que tiene una tendencia positiva (\(\beta_1>0\)) pero sería estacionaria después de eliminar la tendencia (es decir, \(\rho < 1\)).
Otra posibilidad es que \(Y_t\) ha estado creciendo debido a que sigue una caminata aleatoria con un rumbo positivo (es decir, \(\alpha>0\), \(\beta=0\) y \(\rho=1\)); en este caso trabajamos con \(\Delta Y_t\).
La eliminación de la tendencia no haría estacionaria la serie y la inclusión de \(Y_t\) en una regresión (aún si se elimina la tendencia) podría conducir a resultados imprecisos.
Podriamos pensar que la ecuación \(\eqref{eqn:A2.1}\) pudiera ser estimada con OLS y que la estadística \(t\) en \(\hat\rho\) se puede usar para probar si \(\hat\rho\) es significativamente diferente de 1. Sin embargo, si el valor verdadero de \(\rho\) en efecto es 1, el estimador OLS está sesgado hacia cero. Por tanto, al usar OLS de esta manera puede conducirnos a rechazar incorrectamente la hipótesis de la caminata aleatoria.
Dickey y Fuller derivaron la distribución para el estimador \(\hat\rho\) que se cumple cuando \(\rho=1\) y generaron estadísticas para una prueba F simple de la hipótesis de caminata aleatoria; es decir, la hipótesis de que \(\beta=0\) y \(\rho=1\). Suponga que \(Y_t\) puede describirse con la ecuación \(\eqref{eqn:A2.1}\)
Usando OLS, ejecutamos primero la regresión sin restricción
\[\begin{equation} \tag{A2.2} \label{eqn:A2.2} Y_t-Y_{t-1}=\alpha+\beta\cdot t+(\rho-1)\cdot Y_{t-1} +\varepsilon_t \end{equation}\]
y luego la regresión restringida
\[\begin{equation} \tag{A2.3} \label{eqn:A2.3} Y_t-Y_{t-1}=\alpha \end{equation}\]
Luego calculamos la razón F estándar9, donde \(ESS_R\) y \(ESS_{UR}\) son las sumas de cuadrados de los residuales en las regresiones restringida y sin restricciones, respectivamente, \(N\) es el número de observaciones, \(k\) es el número de parámetros estimados en la regresión sin restriccion y \(q\) es el número de restricciones de parámetro) para probar si las restricciones (\(\beta=0\), \(\rho=1\)) se cumplen; esta razón no está distribuida como una distribución F estándar bajo la hipótesis nula, de acuerdo con la distribución calculada por Dickey y Fuller, los valores críticos son mucho mayores que los de la tabla F estándar.
Un problema con la ecuación \(\eqref{eqn:A2.1}\) es que hace la suposición implícita de que no hay correlación serial de ninguna clase en el término de error \(\varepsilon_t\). Nos gustaría permitir, con frecuencia, una correlación serial en \(\varepsilon_t\) y todavía probar por una raíz unitaria. Esto puede hacerse con la prueba de Dickey-Fuller aumentada. Esta prueba se lleva a cabo expandiendo la ecuación \(\eqref{eqn:A2.1}\) para incluir cambios rezagados en \(Y_t\) en el lado derecho de la ecuación:
\[\begin{equation} \tag{A2.4} \label{eqn:A2.4} Y_t=\alpha+\beta\cdot t+\rho\cdot Y_{t-1}+\sum_{j=1}^p\lambda_j\Delta Y_{t-j}+\varepsilon_t \end{equation}\]
donde \(\Delta Y_t=Y_t-Y_{t-1}\).
La prueba de raíz unitaria se realiza así:
Usando OLS, primero se ejecuta la regresión sin restricción
\[\begin{equation} \tag{A2.5} \label{eqn:A2.5} Y_t-Y_{t-1}=\alpha+\beta\cdot t+(\rho-1)\cdot Y_{t-1}+\sum_{j=1}^p\lambda_j\Delta Y_{t-j} \end{equation}\]
y luego la regresión restringida
\[\begin{equation} \tag{A2.6} \label{eqn:A2.6} Y_t-Y_{t-1}=\alpha+\sum_{j=1}^p\lambda_j\Delta Y_{t-j} \end{equation}\]
Entonces, se calcula una razón F estándar para probar si se cumplen las restricciones (\(\beta=0\), \(\rho=1\)). Una vez más, debemos usar la distribución tabulada por Dickey y Fuller.
Debe tomarse en cuenta que la prueba ADF solo nos permite rechazar (o dejar de rechazar) la hipótesis de que una variable no es una caminata aleatoria. Una falla en rechazar (en especial en un nivel de significancia alto) solo proporcina una evidencia débil a favor de la hipótesis de la caminata aleatoria.
Anexo 3: Referencias
https://www.rapidtables.com/code/text/unicode-characters.html
https://www.andrewheiss.com/blog/2018/03/08/amelia-broom-huxtable/
http://haozhu233.github.io/kableExtra/awesome_table_in_html.html
https://plotly-r.com/plotly_book.pdf
https://bookdown.org/yihui/rmarkdown/
Hayashi, F., “Econometrics”, First Ed., Princeton University Press, 2000.
Pindyck, RS, Rubinfeld, D., “Econometría Modelos y Pronósticos”, 4ta. Ed., McGraw-Hill, 2001.