Modelos Univariados, con Desempleo Total de los Estados Unidos de América
Modelos Univariados, con Desempleo Latino

Modelos Univariados, con Desempleo Total de los Estados Unidos de América

Este trabajo tiene como objetivo desarrollar un ejercicio utilizando el ingreso por remesas familiares 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 mas 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. Los resultados no representan ninguna postura institucional, solamente se publican con fines de mostrar una aplicación de pronóstico con modelos ARIMA.

Se elige utilizar los datos sin ajuste estacional debido a que la serie de remesas solamente se tiene en datos históricos. Para tomar en cuenta esta circunstancia dentro de las estimaciones, en los proceso de ARIMA automático y con regresores se incorpora la estimación de coeficientes estacionales.

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(sarima)
library(tidyverse)
library(tsbox)
library(x12)
library(xts)

source("R/functions.R")
set.seed(123)

1) Importar los datos sin ajuste estacional

Se utilizan las remesas familiares y la tasa de desempleo de Estados Unidos de América.

##### Definir fechas #####
inicio <- "2000-01-01"
fin <- "2020-03-01"
dates <- seq(as.Date(inicio), as.Date(fin), by = "month")

##### Remesas Honduras, serie de tiempo #####
Data_HN <- readRDS("RData/Remesas_Mensual.rds") %>%
  dplyr::select(Fecha, Valor)
Remesas_HN = xts::xts(Data_HN[, 2], order.by = dates)
Remesas_HN = ts_ts(ts_span(Remesas_HN, inicio, fin))
dRemesas_HN = ts_diff(log(Remesas_HN))

##### Unemployment USA #####
# https://fred.stlouisfed.org/series/UNRATENSA
fredr_set_key("8828b5fd4b7cf9795f597aca7954c957")
Unemp_USA = fredr::fredr(series_id = "UNRATENSA",
                observation_start = as.Date(inicio)
                )
Unemp_USA = xts::xts(Unemp_USA[, 3], order.by = dates)
Unemp_USA = ts_ts(ts_span(Unemp_USA, inicio, fin))

##### Industrial Production Index USA #####
# https://fred.stlouisfed.org/series/INDPRO
fredr_set_key("8828b5fd4b7cf9795f597aca7954c957")
INDPRO_USA = fredr::fredr(series_id = "INDPRO",
                observation_start = as.Date(inicio)
                )
INDPRO_USA = xts::xts(INDPRO_USA[, 3], order.by = dates)
INDPRO_USA = ts_ts(ts_span(INDPRO_USA, inicio, fin))

2) Inspección de los datos: Remesas_HN

# if(getOutputFormat() == 'html_document') {
  plot_ly(Data_HN, 
          x = ~Fecha, 
          y = ~log(Valor),
          type = 'scatter', mode = 'lines') %>%
    layout(title = "Remesas Honduras (en logs)",
           xaxis = list(title = ""),
           yaxis = list (title = ""))

# } else {
#   ts_plot(
#     `Remesas_HN`= log(Remesas_HN),
#     title = "Remesas Honduras (en logs)",
#     subtitle = "Serie Mensual"
#     )
# }

# 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 = "Remesas Honduras: Tasa de crecimiento mensual",
           xaxis = list(title = ""),
           yaxis = list (title = "%"))

#  } else {
#   ts_plot(
#     `Remesas_HN`= tsbox::ts_pc(Remesas_HN),
#     title = "Remesas Honduras",
#     subtitle = "Tasa de crecimiento respecto al mes anterior"
#     )
# }

# if(getOutputFormat() == 'html_document') {
  plot_ly(Data_HN, x = ~Fecha, y = ~c(rep(NA,12),diff(log(Valor),12))*100,
          type = 'scatter', mode = 'lines') %>%
    layout(title = "Remesas Honduras: Tasa de crecimiento interanual",
           xaxis = list(title = ""),
           yaxis = list (title = "%"))

#  } else {
#   ts_plot(
#     `Remesas_HN`= tsbox::ts_pcy(Remesas_HN),
#     title = "Remesas Honduras",
#     subtitle = "Tasa de crecimiento interanual"
#     )
# }

Conclusiones preliminares: podemos observar que la serie es no estacionaria. Debido a esto, estimaremos el modelo en primeras diferencias. La tasa de crecimiento mensual (que se aproxima a través de la primera diferencia logarítmica) aparenta ser estacionaria. En la serie de variación interanual se observa un aparente cambio estructural a partir de 2010.

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

Se realizó una prueba de Dickey Fuller Aumentada (ADF, por sus siglas en inglés) 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ístico. 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 Remesas en niveles es claramente no estacionario, mientras que su tasa de crecimiento (dlogs) es estacionaria. La misma situación se presenta para el Remesas de Estados Unidos de América (Unemp).

# Criterio en una sola tabla
df <- list(log_Remesas_HN  = as.matrix(log(Remesas_HN)),
           dlog_Remesas_HN = as.matrix(ts_diff(log(Remesas_HN))),
           log_Unemp_USA  = as.matrix(log(Unemp_USA)),
           dlog_Unemp_USA = as.matrix(ts_diff(log(Unemp_USA))),
           log_INDPRO_USA  = as.matrix(log(INDPRO_USA)),
           dlog_INDPRO_USA = as.matrix(ts_diff(log(INDPRO_USA))))
tabla_adf(df, maxLags = 10, lagSelection = "BIC") %>%
  kable_styling(latex_options = "hold_position")

Resultados de Prueba CADF
transf	Nombre	Stat	Drift	None	Trend
log	INDPRO	est.t	-2.6026	0.0140	-3.1058
log	INDPRO	p.value	0.0939	0.6862	0.1074
log	INDPRO	sel.lag	4.0000	3.0000	4.0000
log	Remesas	est.t	-3.0540	2.8099	-2.1174
log	Remesas	p.value	0.0316	0.9989	0.5330
log	Remesas	sel.lag	2.0000	2.0000	2.0000
log	Unemp	est.t	-1.9406	-0.0462	-2.0210
log	Unemp	p.value	0.3133	0.6664	0.5864
log	Unemp	sel.lag	9.0000	9.0000	9.0000
dlog	INDPRO	est.t	-4.2620	-4.2957	-4.1933
dlog	INDPRO	p.value	0.0007	0.0000	0.0054
dlog	INDPRO	sel.lag	2.0000	2.0000	2.0000
dlog	Remesas	est.t	-7.1226	-16.3213	-8.6287
dlog	Remesas	p.value	0.0000	0.0000	0.0000
dlog	Remesas	sel.lag	10.0000	1.0000	10.0000
dlog	Unemp	est.t	-3.8211	-3.8393	-3.9823
dlog	Unemp	p.value	0.0031	0.0001	0.0105
dlog	Unemp	sel.lag	8.0000	8.0000	8.0000

4) Modelos de pronóstico ARIMA

Debido a que la serie 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)¹.

\[\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.

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.

ARIMA automático: Evaluación de los coeficientes

Model_Auto = auto.arima(dRemesas_HN, 
                       max.p = 12, max.q = 12, d = 0, max.P = 12, max.Q = 12,
                       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")

Coeficientes del ARIMA automático
	Betas	Error_Estandar	t_Estadistic	p_value
ar1	-0.6462468	0.0686829	-9.409138	0.0000000
ar2	-0.3532919	0.0784738	-4.502037	0.0000105
ar3	0.1339324	0.0684185	1.957548	0.0514373
sar1	0.5359416	0.0659370	8.128083	0.0000000
sar2	0.2343666	0.0704676	3.325876	0.0010190

coefs_Auto <- coefs
table_stats(Model_Auto, caption = "Estad\u00EDsticos del ARIMA autom\u00E1tico") %>%
  kable_styling(latex_options = "hold_position")

Estadísticos del ARIMA automático
sigma_2	logLik	AIC	AICc	BIC
0.0048473	298.8777	-585.7554	-585.3979	-564.8217

stats_Auto <- cbind(data.frame(Model = "autom\u00E1tico"), stats)

Todos los coeficientes 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")

Estadísticos de residuales ARIMA automático: Training set error measures
ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
0.0040933	0.0688997	0.0516962	0.3856084	165.7589	0.7007538	0.0320732

evalResids_Auto <- cbind(data.frame(Model = "autom\u00E1tico"), eval_resids)

LjungBox_Auto <- checkresiduals(Model_Auto)


    Ljung-Box test

data:  Residuals from ARIMA(3,0,0)(2,0,0)[12] with zero mean
Q* = 31.288, df = 19, p-value = 0.03753

Model df: 5.   Total lags used: 24

autoplot(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 = 1.6243, df = 4, p-value = 0.8044

jarque.bera.test(Model_Auto$residuals[!is.na(Model_Auto$residuals)])


    Jarque Bera Test

data:  Model_Auto$residuals[!is.na(Model_Auto$residuals)]
X-squared = 44.493, df = 2, p-value = 2.181e-10

## OJO MPE y MAPE son diferentes al summary

los residuos no presentan autocorrelación y no se ajustan a una distribución normal,. Nuevamente, se transforman los resultados a tasas de crecimiento interanual.

ARIMA automático: Pronóstico fuera de muestra

For_Auto = forecast(Model_Auto, h = 12) 
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("Remesas Honduras: Pron\u00F3stico ARIMA autom\u00E1tico") +
  labs(caption = "rojo = pron\u00F3stico, azul = estimado, negro = observado",
       x = "", y = "Variaci\u00F3n mensual (%)") +
  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
# }

Variaciones interanuales:

# Paso 1: Unir datos historicos y de pronostico
dfHist <- tsbox::ts_pcy(Remesas_HN)
dfHist <- ts_df(dfHist)
names(dfHist) <- c("date","dfHist")
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_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 = ~dfHist, name = 'Observado', type = 'scatter', mode = 'lines') %>%
    add_trace(x = ~date, y = ~FcstYoY_Auto, name = 'ARIMA Autom\u00E1tico', type = 'scatter', mode = 'lines') %>%
    layout(title = list(text = paste0('Remesas Honduras: Pron\u00F3sticos ARIMA')),
           xaxis = list(title = ""),
           yaxis = list (title = "Variaci\u00F3n interanual (%)"))

#  } else {
#   ts_plot(
#     `Hist.`= ts_span(dfHist, "2000-01-01"),
#     `Autom.` = FcstYoY_Auto,
#     title = "Honduras: Remesas y Pron\u00F3stico ARIMA",
#     subtitle = "Variaci\u00F3n interanual (%)"
#   )
# }

# Tabla de resultados
names(dfFor_Auto) <- c("Fecha","Autom\u00E1tico")
kable(dfFor_Auto,
      caption = "Resultados de pron\u00F3stico: Variaci0nes Interanuales") %>%
  kable_styling(c("striped", "bordered")) %>%
  kable_styling(latex_options = "hold_position")

Resultados de pronóstico: Variaci0nes Interanuales
Fecha	Automático
2020-04-01	-3.069109
2020-05-01	-3.602523
2020-06-01	-9.604693
2020-07-01	-8.655567
2020-08-01	-7.348518
2020-09-01	-12.296534
2020-10-01	-7.926696
2020-11-01	-7.586392
2020-12-01	-10.263048
2021-01-01	-3.783451
2021-02-01	-6.246696
2021-03-01	4.618832

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: Economía Estadounidense

Adicionalmente, se incluye como variable explicativa la tasa de desempleo de Estados Unidos de América (Unemp_USA) en un segundo modelo:

##### Evaluando modelo RegArima #####
dUnemp_USA = ts_diff(log(Unemp_USA))
dUnemp_USA <- ts_xts(dUnemp_USA)
dINDPRO_USA = ts_diff(log(INDPRO_USA))
dINDPRO_USA <- ts_xts(dINDPRO_USA)
dRemesas_HN <- ts_xts(dRemesas_HN)
Model_ArimaX = Arima(dRemesas_HN, order = c(3,0,0), seasonal = list(order = c(2,0,0)), include.constant= TRUE,
                     xreg = ts_c(dUnemp_USA))
table_coef(Model_ArimaX, caption = "Coeficientes del ARIMA, incluyendo Desempleo USA Global") %>%
  kable_styling(latex_options = "hold_position")

Coeficientes del ARIMA, incluyendo Desempleo USA Global
	Betas	Error_Estandar	t_Estadistic	p_value
ar1	-0.6606540	0.1218141	-5.4234607	0.0000001
ar2	0.0249549	0.1597180	0.1562433	0.8759721
ar3	0.3531588	0.0997515	3.5403874	0.0004794
sar1	0.0502057	0.1265481	0.3967324	0.6919159
sar2	-0.2704840	0.1243163	-2.1757721	0.0305444
intercept	0.0114596	0.0033257	3.4457216	0.0006717
value	-0.5423556	0.0637132	-8.5124574	0.0000000

table_stats(Model_ArimaX, caption = "Estad\u00EDsticos del ARIMA, incluyendo  Desempleo USA Global") %>%
  kable_styling(latex_options = "hold_position")

Estadísticos del ARIMA, incluyendo Desempleo USA Global
sigma_2	logLik	AIC	AICc	BIC
0.0067451	264.6755	-513.351	-512.733	-485.4395

table_resids(Model_ArimaX, 
             caption = "Estad\u00EDsticos de residuales ARIMA incl. Regresores: Training set error measures") %>%
  kable_styling(latex_options = "hold_position")

Estadísticos de residuales ARIMA incl. Regresores: Training set error measures
ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
0.0006584	0.0809321	0.0659862	-1.013889	173.0477	0.8944582	0.0494775

LjungBox_ArimaX <- checkresiduals(Model_ArimaX)


    Ljung-Box test

data:  Residuals from Regression with ARIMA(3,0,0) errors
Q* = 13.512, df = 3, p-value = 0.00365

Model df: 7.   Total lags used: 10

autoplot(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.3717, df = 4, p-value = 0.6677

jarque.bera.test(Model_ArimaX$residuals[!is.na(Model_ArimaX$residuals)])


    Jarque Bera Test

data:  Model_ArimaX$residuals[!is.na(Model_ArimaX$residuals)]
X-squared = 3.9001, df = 2, p-value = 0.1423

evalResids_ArimaX <- cbind(data.frame(Model = "ARIMA incl. Desempleo USA Global"), eval_resids)

6) Comparación de modelos

Prueba de Ljung-Box

Esta prueba se utiliza para verificar correlación serial en los residuales².

\(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 Autom\u00E1tico","Arima con Desempleo USA Global"),
                       Q = c(LjungBox_Auto$statistic,LjungBox_ArimaX$statistic),
                       df = c(LjungBox_Auto$parameter,LjungBox_ArimaX$parameter),
                       p_value = c(LjungBox_Auto$p.value,LjungBox_ArimaX$p.value))
kable(LjungBox,
      caption = "Ljung-Box test") %>%
  kable_styling(c("striped", "bordered")) %>%
  kable_styling(latex_options = "hold_position")

Ljung-Box test
Modelos	Q	df	p_value
ARIMA Automático	31.28753	19	0.0375335
Arima con Desempleo USA Global	13.51215	3	0.0036503

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\).³

BoxLjung <- data.frame(Modelos = c("ARIMA Autom\u00E1tico","Arima con Desempleo USA Global"),
                       ChiSqrt = c(BoxLjung_Auto$statistic,BoxLjung_ArimaX$statistic),
                       df = c(BoxLjung_Auto$parameter,BoxLjung_ArimaX$parameter),
                       p_value = c(BoxLjung_Auto$p.value,BoxLjung_ArimaX$p.value))
kable(BoxLjung,
      caption = "Prueba de Box-Ljung") %>%
  kable_styling(c("striped", "bordered")) %>%
  kable_styling(latex_options = "hold_position")

Prueba de Box-Ljung
Modelos	ChiSqrt	df	p_value
ARIMA Automático	1.624347	4	0.8044098
Arima con Desempleo USA Global	2.371739	4	0.6677404

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_Auto,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")

Evaluación de Residuales: Pronóstico dentro de Muestra
Model	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
automático	0.004093	0.068900	0.051696	0.385608	165.7589	0.700754	0.032073
ARIMA incl. Desempleo USA Global	0.000658	0.080932	0.065986	-1.013889	173.0477	0.894458	0.049478

En la tabla puede observarse que el modelo que tiene los menores valores de los estadísticos (mejor modelo) es el ARIMA automático; sin embargo, se usará el modelo con regresores para simular impactos de la tasa de desempleo de los Estados Unidos en el pronóstico.

7) Fan chart de Remesas

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.

7.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 dUnemp_USA igual a 5% desde abril de 2020 hasta diciembre 2020.

dRemesas_HN <- 100*ts_diff(log(Remesas_HN))
ModelTest = Arima(dRemesas_HN, order = c(3,0,0), seasonal = list(order = c(2,0,0)), include.constant= TRUE,
                     xreg = ts_c(dUnemp_USA))
inicio <- "2020-04-01"
fin <- "2020-12-01"
dates_for <- seq(as.Date(inicio), as.Date(fin), by = "month")
dUnemp_USA <- data.frame(dates_for, dUnemp_USA = c(rep(0.05, 9)))
dUnemp_USA = xts::xts(dUnemp_USA[, 2], order.by = dates_for)
xregres = ts_span(ts_c(dUnemp_USA), NULL, as.Date(fin))

fcstRemesas_HN = forecast(ModelTest, xreg = xregres, h = 9, level = c(50, 70, 80, 95))
IntervalRange = fcstRemesas_HN$upper-fcstRemesas_HN$lower
colnames(IntervalRange) = c("50%", "70%", "80%", "95%")

# if(getOutputFormat() == 'html_document') {
  dfHist <- data.frame(date=as.Date(time(dRemesas_HN)), HistYoY=as.matrix(dRemesas_HN))
  dfMid <- data.frame(date=as.Date(time(fcstRemesas_HN$mean)), AutoFcstYoY=as.matrix(fcstRemesas_HN$mean))
  dfLow <- data.frame(date=as.Date(time(fcstRemesas_HN$lower)), AutoFcstYoY=as.matrix(fcstRemesas_HN$lower))
  dfUp <- data.frame(date=as.Date(time(fcstRemesas_HN$upper)), AutoFcstYoY=as.matrix(fcstRemesas_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('Pron\u00F3stico ARIMA, 
                                      Intervalos de Confianza modelo con Regresores',
                                      '<br>',
                                      '<sup>',
                                      'Variaci\u00F3n interanual')),
           xaxis = list(title = ""),
           yaxis = list (title = "%"))

 # } else {
 #  autoplot(fcstRemesas_HN) +
 #    theme_minimal() +
 #    theme(text=element_text(size=8))
 # }
# kable(fcstRemesas_HN)

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).

7.2) Intervalos de Confianza

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(fcstRemesas_HN)
autoplot(sigh2) +
  theme_minimal() + 
  geom_hline(yintercept=var(dRemesas_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 Remesas 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 = fcstRemesas_HN$mean
SimFcst = matrix(NA, nrow = 9, ncol = NSim)
for (h in 1:9){
  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,fcstRemesas_HN$lower[, "95%"], fcstRemesas_HN$upper[, "95%"])
plot(allci)

c) Calcular probabilidad de crecimiento dentro de un rango del Remesas 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 -10% en 2020q1: ",PNeg2020 = mean(SimFcst["2020-04-01",] <= -10)))

[1] "Probabilidad de crecimiento menor a -10% en 2020q1: 0"

# Probabilidad de crecimiento del Remesas menor a -14% in Q1 2020
(paste0("Probabilidad de crecimiento menor a -14% en 2020q1: ",PLarge2020 = mean(SimFcst["2020-04-01",] <= -14)))

[1] "Probabilidad de crecimiento menor a -14% en 2020q1: 0"

# Probailidad de crecimiento del Remesas entre -10% y -14% in Q1 2020
(paste0("Probabilidad de crecimiento entre -10% y -14% en junio 2020: ",
        PBetween2020 = mean(SimFcst["2020-04-01",] <= -10 & SimFcst["2020-06-01",] >= -14)))

[1] "Probabilidad de crecimiento entre -10% y -14% en junio 2020: NaN"

# Graficar barras
PNeg = ts_rowMeans(SimFcst <= -10)
PLarge = ts_rowMeans(SimFcst <= -14)
PBetween = ts_rowMeans(SimFcst <= -10 & SimFcst >= -14)

prob.data = data.frame(index(PNeg), PNeg, PLarge, PBetween)
colnames(prob.data) = c("Date", "<= -10%", "<= -14%", "entre -10% y -14%")
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 Remesas en 2020 sea menor a -14%

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(dRemesas_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_10 = mean(SimFcstGrt["2020-01-01",] <= -10)
p2020_14 = mean(SimFcstGrt["2020-01-01",] <= -14)

# 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('Remesas 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%"],
#     `Remesas`= SimFcstGrt["1981-01-01/2021-12-01",1],
#     title = "Remesas Honduras",
#     subtitle = "Tasa de crecimiento anual"
#   )
# }

Los resultados indican que la probabilidad que la tasa de crecimiento anual del Remesas sea menor a -10% a finales de 2020 es 71% y menor a -14% a finales de 2020 es 54% . La mediana de los datos se ubica en -14%, valor similar al observado durante 2009; adicionalmente, la distribución del pronóstico es asimétrica, con valores de la distribución concentrados cerca de la mediana.

Modelos Univariados, con Desempleo Latino

1) Economía Estadounidense, Desempleo Latino

Adicionalmente, se incluye como variable explicativa la tasa de desempleo latino en los Estados Unidos de América (Unemp_USA_Latin) en un segundo modelo. Puede verse que existe una alta correlación entre el desempleo global y el latino:

## Desempleo latino
# Fuente: https://data.bls.gov/pdq/SurveyOutputServlet
Unemp_USA_Latin <- read_csv("XLSData/Unemp_USA_Latin.txt")
Unemp_USA_Latin <- Unemp_USA_Latin %>%
  tidyr::gather("month", "Latino", 2:ncol(.)) %>%
  dplyr::mutate(months = case_when(
    month == 'Jan' ~ 1,
    month == 'Feb' ~ 2,
    month == 'Mar' ~ 3,
    month == 'Apr' ~ 4,
    month == 'May' ~ 5,
    month == 'Jun' ~ 6,
    month == 'Jul' ~ 7,
    month == 'Aug' ~ 8,
    month == 'Sep' ~ 9,
    month == 'Oct' ~ 10,
    month == 'Nov' ~ 11,
    month == 'Dec' ~ 12),
    dates = as.Date(paste0(Year,"/",months,"/01"))) %>%
  dplyr::arrange(dates) %>%
  tidyr::drop_na()
## Desempleo global
inicio <- "1990-01-01"
fredr_set_key("8828b5fd4b7cf9795f597aca7954c957")
Unemp_USA = fredr::fredr(series_id = "UNRATENSA",
                observation_start = as.Date(inicio)
                ) 
Unemp_USA <- cbind(Unemp_USA,Unemp_USA_Latin$Latino) %>%
  select(-series_id)
names(Unemp_USA) <- c("Fecha","Global","Latino")
corr_Unemp <- cor(Unemp_USA$Global,Unemp_USA$Latino)

## Plots
Unemp_USA <- Unemp_USA %>%
  tidyr::gather("Desempleo", "Valor", 2:ncol(.))
plot1 <- ggplot(Unemp_USA) +
  geom_line(aes(x = Fecha, y = Valor, colour = Desempleo)) +
  theme_minimal()
plot1

## Serie de tiempo
inicio <- "2000-01-01"
fin <- "2020-03-01"
Unemp_USA_Latin <- Unemp_USA_Latin %>%
  filter(dates >= inicio,
         dates <= fin)
dates <- seq(as.Date(inicio), as.Date(fin), by = "month")
Unemp_USA_Latin = xts::xts(Unemp_USA_Latin[, 3], order.by = dates)
Unemp_USA_Latin = ts_ts(ts_span(Unemp_USA_Latin, inicio, fin))
dUnemp_USA_Latin = ts_diff(log(Unemp_USA_Latin))
dUnemp_USA_Latin <- ts_xts(dUnemp_USA_Latin)

El coeficiente de correlación entre las series es bastante alto, 0.8986883, asimismo se observa que a lo largo del tiempo el desempleo latino es mayor.

##### Evaluando modelo RegArima #####
Data_HN2 <- readRDS("RData/Remesas_Mensual.rds") %>%
  dplyr::select(Fecha, Valor) %>%
  dplyr::filter(Fecha <= as.Date("2020/04/01"))
Remesas_HN = xts::xts(Data_HN2[, 2], order.by = dates)
Remesas_HN = ts_ts(ts_span(Remesas_HN, inicio, fin))
dRemesas_HN = ts_diff(log(Remesas_HN))

Model_ArimaX_Latin = Arima(dRemesas_HN, order = c(3,0,0), seasonal = list(order = c(2,0,0)), 
                           include.constant= TRUE,
                     xreg = ts_c(dUnemp_USA_Latin))
table_coef(Model_ArimaX_Latin, caption = "Coeficientes del ARIMA, incluyendo Desempleo Latino") %>%
  kable_styling(latex_options = "hold_position")

Coeficientes del ARIMA, incluyendo Desempleo Latino
	Betas	Error_Estandar	t_Estadistic	p_value
ar1	-0.7044818	0.0679796	-10.363142	0.0000000
ar2	-0.4499321	0.0787151	-5.715955	0.0000000
ar3	0.0845068	0.0684107	1.235286	0.2179271
sar1	0.4968730	0.0647280	7.676324	0.0000000
sar2	0.2790474	0.0695787	4.010529	0.0000809
intercept	0.0128208	0.0076314	1.680007	0.0942517
value	-0.2647885	0.0666385	-3.973506	0.0000936

table_stats(Model_ArimaX_Latin, caption = "Estad\u00EDsticos del ARIMA, incluyendo Desempleo Latino") %>%
  kable_styling(latex_options = "hold_position")

Estadísticos del ARIMA, incluyendo Desempleo Latino
sigma_2	logLik	AIC	AICc	BIC
0.004555	307.2246	-598.4492	-597.8312	-570.5377

table_resids(Model_ArimaX_Latin, 
             caption = "Estad\u00EDsticos de residuales ARIMA incl. Regresores (Desempleo Latino): Training set error measures") %>%
  kable_styling(latex_options = "hold_position")

Estadísticos de residuales ARIMA incl. Regresores (Desempleo Latino): Training set error measures
ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
-0.0021192	0.0665077	0.0511233	-16.58912	164.0224	0.6929883	0.0376032

LjungBox_ArimaX_Latin <- checkresiduals(Model_ArimaX_Latin)


    Ljung-Box test

data:  Residuals from Regression with ARIMA(3,0,0)(2,0,0)[12] errors
Q* = 28.597, df = 17, p-value = 0.03843

Model df: 7.   Total lags used: 24

autoplot(Model_ArimaX_Latin) + theme_minimal()

BoxLjung_ArimaX_Latin <- Box.test(Model_ArimaX_Latin$residuals, lag=5, fitdf=1, type="Lj")
BoxLjung_ArimaX_Latin


    Box-Ljung test

data:  Model_ArimaX_Latin$residuals
X-squared = 2.9713, df = 4, p-value = 0.5626

jarque.bera.test(Model_ArimaX_Latin$residuals[!is.na(Model_ArimaX_Latin$residuals)])


    Jarque Bera Test

data:  Model_ArimaX_Latin$residuals[!is.na(Model_ArimaX_Latin$residuals)]
X-squared = 18.972, df = 2, p-value = 7.592e-05

evalResids_ArimaX_Latin <- cbind(data.frame(Model = "ARIMA incl. Regresores, Desempleo Latino"), eval_resids)

2) Comparación de modelos

Prueba de Ljung-Box

LjungBox <- data.frame(Modelos = c("ARIMA Autom\u00E1tico","Arima con Desempleo USA Global",
                                   "Arima con Desempleo USA Latino"),
                       Q = c(LjungBox_Auto$statistic,LjungBox_ArimaX$statistic,
                             LjungBox_ArimaX_Latin$statistic),
                       df = c(LjungBox_Auto$parameter,LjungBox_ArimaX$parameter,
                              LjungBox_ArimaX_Latin$parameter),
                       p_value = c(LjungBox_Auto$p.value,LjungBox_ArimaX$p.value,
                                   LjungBox_ArimaX_Latin$p.value))
kable(LjungBox,
      caption = "Ljung-Box test") %>%
  kable_styling(c("striped", "bordered")) %>%
  kable_styling(latex_options = "hold_position")

Ljung-Box test
Modelos	Q	df	p_value
ARIMA Automático	31.28753	19	0.0375335
Arima con Desempleo USA Global	13.51215	3	0.0036503
Arima con Desempleo USA Latino	28.59711	17	0.0384268

Prueba de Box-Ljung

BoxLjung <- data.frame(Modelos = c("ARIMA Autom\u00E1tico","Arima con Desempleo USA Global",
                                   "Arima con Desempleo USA Latino"),
                       ChiSqrt = c(BoxLjung_Auto$statistic,BoxLjung_ArimaX$statistic,
                                   BoxLjung_ArimaX_Latin$statistic),
                       df = c(BoxLjung_Auto$parameter,BoxLjung_ArimaX$parameter,
                              BoxLjung_ArimaX_Latin$statistic),
                       p_value = c(BoxLjung_Auto$p.value,BoxLjung_ArimaX$p.value,
                                   BoxLjung_ArimaX_Latin$p.value))
kable(BoxLjung,
      caption = "Prueba de Box-Ljung") %>%
  kable_styling(c("striped", "bordered")) %>%
  kable_styling(latex_options = "hold_position")

Prueba de Box-Ljung
Modelos	ChiSqrt	df	p_value
ARIMA Automático	1.624347	4.000000	0.8044098
Arima con Desempleo USA Global	2.371739	4.000000	0.6677404
Arima con Desempleo USA Latino	2.971327	2.971327	0.5626352

Evaluación del pronóstico, residuales dentro de muestra

evalResids <- rbind(evalResids_Auto,evalResids_ArimaX,evalResids_ArimaX_Latin)
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")

Evaluación de Residuales: Pronóstico dentro de Muestra
Model	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
automático	0.004093	0.068900	0.051696	0.385608	165.7589	0.700754	0.032073
ARIMA incl. Desempleo USA Global	0.000658	0.080932	0.065986	-1.013889	173.0477	0.894458	0.049478
ARIMA incl. Regresores, Desempleo Latino	-0.002119	0.066508	0.051123	-16.589119	164.0224	0.692988	0.037603

3) Fan chart de Remesas, Desempleo Latino

3.1) Estimación de un modelo ARMA y fan chart

Con el fin de simular un impacto por el efecto del COVID-19 durante 2020, se realiza un pronóstico asumiendo comportamiento de dUnemp_USA_Latin igual a 6% desde marzo de 2020 hasta diciembre 2020.

dRemesas_HN <- 100*ts_diff(log(Remesas_HN))
ModelTest = Arima(dRemesas_HN, order = c(3,0,0), seasonal = list(order = c(2,0,0)), include.constant= TRUE,
                     xreg = ts_c(dUnemp_USA_Latin))
inicio <- "2020-04-01"
fin <- "2020-12-01"
dates_for <- seq(as.Date(inicio), as.Date(fin), by = "month")
dUnemp_USA_Latin <- data.frame(dates_for, dUnemp_USA_Latin = c(rep(0.06, 9)))
dUnemp_USA_Latin = xts::xts(dUnemp_USA_Latin[, 2], order.by = dates_for)
xregres = ts_span(ts_c(dUnemp_USA_Latin), NULL, as.Date(fin))

fcstRemesas_HN_Latin = forecast(ModelTest, xreg = xregres, h = 9, level = c(50, 70, 80, 95))
IntervalRange = fcstRemesas_HN_Latin$upper-fcstRemesas_HN_Latin$lower
colnames(IntervalRange) = c("50%", "70%", "80%", "95%")

# if(getOutputFormat() == 'html_document') {
  dfHist <- data.frame(date=as.Date(time(dRemesas_HN)), HistYoY=as.matrix(dRemesas_HN))
  dfMid <- data.frame(date=as.Date(time(fcstRemesas_HN_Latin$mean)), AutoFcstYoY=as.matrix(fcstRemesas_HN_Latin$mean))
  dfLow <- data.frame(date=as.Date(time(fcstRemesas_HN_Latin$lower)), AutoFcstYoY=as.matrix(fcstRemesas_HN_Latin$lower))
  dfUp <- data.frame(date=as.Date(time(fcstRemesas_HN_Latin$upper)), AutoFcstYoY=as.matrix(fcstRemesas_HN_Latin$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('Pron\u00F3stico ARIMA, 
                                      Intervalos de Confianza modelo con Desempleo Latino',
                                      '<br>',
                                      '<sup>',
                                      'Variaci\u00F3n interanual')),
           xaxis = list(title = ""),
           yaxis = list (title = "%"))

 # } else {
 #  autoplot(fcstRemesas_HN) +
 #    theme_minimal() +
 #    theme(text=element_text(size=8))
 # }

3.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

sigh2 = getForecastVariance(fcstRemesas_HN_Latin)
autoplot(sigh2) +
  theme_minimal() + 
  geom_hline(yintercept=var(dRemesas_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

NSim = 100
fcsth = fcstRemesas_HN_Latin$mean
SimFcst = matrix(NA, nrow = 9, ncol = NSim)
for (h in 1:9){
  SimFcst[h, ] = rnorm(NSim, fcsth[h], sqrt(sigh2[h]))
}
SimFcst = xts(SimFcst, order.by = as.Date(index(ts_xts(fcsth))))

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,fcstRemesas_HN_Latin$lower[, "95%"], fcstRemesas_HN_Latin$upper[, "95%"])
plot(allci)

d) Probabilidad que el crecimiento anual del Remesas en 2020 sea menor a -14%

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(dRemesas_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_10 = mean(SimFcstGrt["2020-01-01",] <= -10)
p2020_14 = mean(SimFcstGrt["2020-01-01",] <= -14)

# 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('Remesas Honduras: Pron\u00F3stico considerando Desempleo Latino',
                                      '<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%"],
#     `Remesas`= SimFcstGrt["1981-01-01/2021-12-01",1],
#     title = "Remesas Honduras",
#     subtitle = "Tasa de crecimiento anual"
#   )
# }

Los resultados indican que la probabilidad que la tasa de crecimiento anual del Remesas sea menor a -10% a finales de 2020 es 60% y menor a -14% a finales de 2020 es 52%; adicionalmente, la distribución del pronóstico es asimétrica, con valores de la distribución concentrados cerca de la mediana.

Hayashi (2000)↩
https://otexts.com/fpp2/residuals.html ↩
Hayashi (2000), pgs.142-144↩

Pronosticar Remesas Mensuales en Honduras

Modelos Univariados, con Desempleo Total de los Estados Unidos de América

1) Importar los datos sin ajuste estacional

2) Inspección de los datos: Remesas_HN

3) Pruebas de raíz unitaria

4) Modelos de pronóstico ARIMA

ARIMA automático: Evaluación de los coeficientes

ARIMA automático: Análisis de los residuales

ARIMA automático: Pronóstico fuera de muestra

5) Extensión del modelo: Economía Estadounidense

6) Comparación de modelos

Prueba de Ljung-Box

Prueba de Box-Ljung

Evaluación del pronóstico, residuales dentro de muestra

7) Fan chart de Remesas

7.1) Estimación de un modelo ARMA y fan chart

7.2) Intervalos de Confianza

a) Extraer la varianza del error de pronóstico y comparar con la varianza incondicional

b) Simular la densidad de pronóstico y comparar con el intervalo de pronóstico analítico obtenido

c) Calcular probabilidad de crecimiento dentro de un rango del Remesas en Q2 2020

d) Probabilidad que el crecimiento anual del Remesas en 2020 sea menor a -14%

Modelos Univariados, con Desempleo Latino

1) Economía Estadounidense, Desempleo Latino

2) Comparación de modelos

Prueba de Ljung-Box

Prueba de Box-Ljung

Evaluación del pronóstico, residuales dentro de muestra

3) Fan chart de Remesas, Desempleo Latino

3.1) Estimación de un modelo ARMA y fan chart

3.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

b) Simular la densidad de pronóstico y comparar con el intervalo de pronóstico analítico obtenido

d) Probabilidad que el crecimiento anual del Remesas en 2020 sea menor a -14%