Quiz 2 Econometría II - María José Guerrero - María José Otero

Descripción de los datos

Para el presente trabajo se utilizarán 3 series de tiempo extraídas de un archivo de GoogleSheets específico. En la serie de tiempo se podrán encontrar las siguientes variables, todas con una incidencia mensual desde el año 1990 hasta el 2023.

• UNRATE (Unemployment Rate): Representa el número de desempleados como porcentaje de la fuerza laboral de Estados Unidos. Los datos de la fuerza laboral están restringidos a personas mayores de 16 años, que actualmente residen en 1 de los 50 estados

• T10Y2YM (10-Year Treasury Constant Maturity Minus 2-Year Treasury Constant Maturity): También conocido como el CMT Index, indica el valor teórico de un Tesoro de Estados Unidos que se basa en los valores recientes de los bonos del Tesos de EEUU subastados

• FEDFUND (Federal Funds Effective Rate): Es la tasa de interés a la que las instituciones de depósito (bancos, cajas de ahorro y cooperativas de credito) intercambian fondos federales (saldos mantenidos en los bancos de la Reserva Federal) entre sí durante la noche. Cuando una institución de depósito tiene saldos excedentes en su cuenta de reserva, presta a otros bancos que necesitan saldos mayores. En términos más simples, un banco con exceso de efectivo, lo que a menudo se conoce como liquidez, prestará a otro banco que necesite aumentar rápidamente la liquidez.

Gráfica inicial de la serie para analizar estacionariedad

subdata <- subset(brent, year>="1990/01/01" & year<="2023/01/01")
View(subdata)

subdata2 <- subset(brent, year>="1990/01/01" & year<="2022/09/01")
View(subdata2)


plot.ts(subdata[, c("FEDFUND", "T10Y2YM", "UNRATE")], main = "Análisis gráfico")

En general, se puede observar que las tres variables aparentan ser relativamente variables, siendo la Tasa de Interés de los Bonos la que resalta con mayor variabilidad a lo largo del tiempo. Por su parte, esta última variable parece ser una serie no estacionaria en media y varianza.

Por otro lado, únicamente desde el lado del análisis gráfico, se podría decir que la Tasa de Desempleo es una serie estacionaria en media pero no en varianza y que la Tasa Efectiva de la Reserva podría tratarse de otra serie no estacionaria en media y varianza

Análisis Correlogramas

ggarrange((ggAcf(subdata$UNRATE, main = "AC de UNRATE", col="purple")),
          (ggPacf(subdata$UNRATE, main = "PAC de UNRATE", col="purple")),
          (ggAcf(subdata$FEDFUND, main = "AC de FEDFUND", col="purple")),
          (ggPacf(subdata$FEDFUND, main = "PAC de FEDFUND", col="purple")),
          (ggAcf(subdata$T10Y2YM, main = "AC de T10Y2YM", col="purple")),
          (ggPacf(subdata$T10Y2YM, main = "PAC de T10Y2YM", col="purple"))
          + 
            rremove("x.text"), ncol = 2, nrow = 3)

Podemos observar los siguientes resultados:

• La serie de la Tasa de Desempleo corresponde a un proceso autorregresivo de primer orden (AR(1))

• La serie de la Tasa de Interés de los Bonos corresponde a un proceso de tercer orden (AR(3)) con picos significativos en los rezagos 9 y 16

• La serie de la Tasa Efectiva de la Reserva corresponde a un proceso autorregresivo de tercer orden (AR(3))

Test de Dickey-Fuller

DF.FEDFUND <- ur.df(subdata$FEDFUND, type = "trend", selectlags = "AIC") 
summary(DF.FEDFUND)

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression trend 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.95216 -0.03708 -0.00621  0.05071  0.52884 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.079e-02  3.108e-02   0.669   0.5040    
## z.lag.1     -7.929e-03  4.689e-03  -1.691   0.0917 .  
## tt          -1.385e-05  9.790e-05  -0.141   0.8875    
## z.diff.lag   6.305e-01  3.980e-02  15.841   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1441 on 391 degrees of freedom
## Multiple R-squared:  0.4084, Adjusted R-squared:  0.4039 
## F-statistic: 89.97 on 3 and 391 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -1.6908 2.005 2.9058 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau3 -3.98 -3.42 -3.13
## phi2  6.15  4.71  4.05
## phi3  8.34  6.30  5.36

Como se logra observar el test estadístico es de 2.9058 y el phi 3 toma valores de 8.34, 6.30 y 5.36 para cada uno de los niveles de significancia, esto demuestra que la hipótesis nula del modelo se rechaza y se presenta una serie no estacionaria.

DF.UNRATE <- ur.df(subdata$UNRATE, type = "trend", selectlags = "AIC") 
summary(DF.UNRATE)

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression trend 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.6837 -0.1375 -0.0267  0.0901 10.2017 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.3472409  0.1147545   3.026 0.002642 ** 
## z.lag.1     -0.0553526  0.0166868  -3.317 0.000995 ***
## tt          -0.0001461  0.0002500  -0.585 0.559206    
## z.diff.lag   0.0530132  0.0506017   1.048 0.295444    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.566 on 391 degrees of freedom
## Multiple R-squared:  0.02845,    Adjusted R-squared:  0.02099 
## F-statistic: 3.816 on 3 and 391 DF,  p-value: 0.01021
## 
## 
## Value of test-statistic is: -3.3171 3.7427 5.6005 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau3 -3.98 -3.42 -3.13
## phi2  6.15  4.71  4.05
## phi3  8.34  6.30  5.36

Se puede evidenciar que test estadístico es de 5.6005 y el phi 3 toma valores de 8.34, 6.30 y 5.36 para cada uno de los niveles de significancia, esto demuestra que la hipótesis nula del modelo se rechaza para todos los percentiles excepto del decimo, es decir, se presenta una serie no estacionaria para el percentil uno y el quinto pero para el decimo es estacionario.

DF.T10Y2YM <- ur.df(subdata$T10Y2YM, type = "trend", selectlags = "AIC") 
summary(DF.T10Y2YM)

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression trend 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.73643 -0.07320 -0.00426  0.06185  0.55638 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.021e-02  1.621e-02   1.863   0.0632 .  
## z.lag.1     -1.507e-02  7.822e-03  -1.926   0.0548 .  
## tt          -7.429e-05  5.990e-05  -1.240   0.2156    
## z.diff.lag   3.171e-01  4.815e-02   6.584 1.47e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1353 on 391 degrees of freedom
## Multiple R-squared:  0.1097, Adjusted R-squared:  0.1029 
## F-statistic: 16.06 on 3 and 391 DF,  p-value: 7.28e-10
## 
## 
## Value of test-statistic is: -1.9263 1.772 2.6384 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau3 -3.98 -3.42 -3.13
## phi2  6.15  4.71  4.05
## phi3  8.34  6.30  5.36

Es posible observar que el test estadístico es de 2.6384 y el phi 3 toma valores de 8.34, 6.30 y 5.36 para cada uno de los niveles de significancia, esto demuestra que la hipótesis nula del modelo se rechaza debido a que el valor dado por el test es menor que los presentados por phi3, por lo que se presenta una serie no estacionaria.

Luego de evidenciar que todas las variables son no estacionarias se pretende diferenciarlas las veces que sean necesarias hasta que se llegue a una serie estacionaria en todos sus niveles de significancia.

#FEDFUND
subdata <- subdata %>%
  mutate(dFEDFUND = c(NA,diff(subdata$FEDFUND)))

DF1 <- ur.df(subdata[!is.na(subdata$dFEDFUND),]$dFEDFUND, type="none", selectlags ="AIC")
summary(DF1)

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression none 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.94584 -0.03660  0.00169  0.04034  0.47424 
## 
## Coefficients:
##            Estimate Std. Error t value Pr(>|t|)    
## z.lag.1    -0.30275    0.04282  -7.070 7.14e-12 ***
## z.diff.lag -0.16937    0.04997  -3.389 0.000772 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1429 on 392 degrees of freedom
## Multiple R-squared:  0.2051, Adjusted R-squared:  0.2011 
## F-statistic: 50.59 on 2 and 392 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -7.0702 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62

Después de haber diferenciado la serie una vez, se procede a volver a realizar el test de Dickey Fuller para identificar si hay raíz unitaria. Para este caso, se observa que el test da un resultado de -7.0702 lo que es mayor a los niveles de significancia de los diferentes percentiles -2.58, -1.95, -1.62. Demostrando así que ya es una serie estacionaria

#UNRATE
subdata <- subdata %>%
  mutate(dUNRATE = c(NA,diff(subdata$UNRATE)))

DF2 <- ur.df(subdata[!is.na(subdata$dUNRATE),]$dUNRATE, type="none", selectlags ="AIC")
summary(DF2)

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression none 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.0133 -0.1107 -0.0121  0.0893 10.2753 
## 
## Coefficients:
##            Estimate Std. Error t value Pr(>|t|)    
## z.lag.1    -1.07947    0.07013  -15.39   <2e-16 ***
## z.diff.lag  0.10695    0.05022    2.13   0.0338 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.57 on 392 degrees of freedom
## Multiple R-squared:  0.4935, Adjusted R-squared:  0.4909 
## F-statistic: 190.9 on 2 and 392 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -15.392 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62

Se ha diferenciado la serie una vez por lo que se procede a volver a realizar el test de Dickey Fuller para identificar si hay raíz unitaria. Para este caso, se observa que el test da un resultado de -15.392 lo que es mayor a los niveles de significancia de los diferentes percentiles -2.58, -1.95, -1.62. Demostrando así que ya es una serie estacionaria

#T10Y2YM
subdata <- subdata %>%
  mutate(dT10Y2YM = c(NA,diff(subdata$T10Y2YM)))

DF3 <- ur.df(subdata[!is.na(subdata$dT10Y2YM),]$dT10Y2YM, type="none", selectlags ="AIC")
summary(DF3)

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression none 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.70626 -0.07086 -0.00426  0.05688  0.54500 
## 
## Coefficients:
##            Estimate Std. Error t value Pr(>|t|)    
## z.lag.1    -0.75139    0.05890 -12.758   <2e-16 ***
## z.diff.lag  0.09271    0.05024   1.845   0.0657 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1352 on 392 degrees of freedom
## Multiple R-squared:  0.3501, Adjusted R-squared:  0.3468 
## F-statistic: 105.6 on 2 and 392 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -12.758 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62

## Se diferencio una vez para volver la serie estacionaria 

Por último se ha diferenciado la serie una vez para así proceder a volver a realizar el test de Dickey Fuller con el cual se logrará identificar si hay raíz unitaria en el modelo. Para este caso, se observa que el test da un resultado de -12.758 lo que es mayor a los niveles de significancia de los diferentes percentiles -2.58, -1.95, -1.62. Demostrando así que ya es una serie estacionaria.

En conclusión, hemos pasado cada una de las variables escogidas que presentaban series no estacionarias a series con estacionaridad.

Modelo de Rezago Distribuido

Un Modelo de Rezagos Distribuidos (ADL) son regresiones estándar de mínimos cuadrados que incluyen retrasos tanto de la variable dependiente como de las variables explicativas como regresores. Es un método para examinar las relaciones de cointegración entre variables.

Un modelo econométrico de este tipo es utilizado en series temporales en las que una o varias variables explicativas pueden tener efectos sobre la variable dependiente tras uno o más periodos.

Es por esto que se va iniciar haciendo el analisis del modelo con los rezagos de las variables esogidas.

# Dos rezagos en T10Y2YM
rezagod <- dlm(formula = dFEDFUND ~ dT10Y2YM, 
               data = data.frame(subdata), q = 2)
summary(rezagod)

## 
## Call:
## lm(formula = as.formula(model.formula), data = design)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.80893 -0.06632  0.00458  0.08397  0.59548 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.011844   0.007905  -1.498 0.134882    
## dT10Y2YM.t  -0.508590   0.058592  -8.680  < 2e-16 ***
## dT10Y2YM.1  -0.211334   0.061603  -3.431 0.000667 ***
## dT10Y2YM.2  -0.312495   0.058532  -5.339 1.59e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1569 on 390 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.3004, Adjusted R-squared:  0.2951 
## F-statistic: 55.83 on 3 and 390 DF,  p-value: < 2.2e-16
## 
## AIC and BIC values for the model:
##         AIC       BIC
## 1 -335.4749 -315.5931

Los rezagos dados por el modelo demuestran que los valores p < 2e-16 para los coeficientes de dT10Y2YM.t, dT10Y2YM.1 y dT10Y2YM.2 indican que todas estas variables predictoras son altamente significativas y tienen un efecto negativo en la variable de respuesta. Es decir, un aumento en cualquiera de estas variables predictoras se asocia con una disminución en la variable de respuesta.

# Modelo de rezago distribuido: método de Koyck
koyck <- koyckDlm(subdata$dT10Y2YM , subdata$dFEDFUND)
summary(koyck, diagnostic=T)

## 
## Call:
## "Y ~ (Intercept) + Y.1 + X.t"
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.818363 -0.064612  0.004479  0.062796  0.420666 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.005184   0.006809  -0.761   0.4469    
## Y.1          0.543633   0.052221  10.410   <2e-16 ***
## X.t         -0.491048   0.205231  -2.393   0.0172 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1342 on 392 degrees of freedom
## Multiple R-Squared: 0.4858,  Adjusted R-squared: 0.4832 
## Wald test: 155.2 on 2 and 392 DF,  p-value: < 2.2e-16 
## 
## Diagnostic tests:
##                  df1 df2  statistic      p-value
## Weak instruments   1 392 23.4357538 1.861385e-06
## Wu-Hausman         1 391  0.1898209 6.633062e-01
## 
##                                alpha       beta       phi
## Geometric coefficients:  -0.01136002 -0.4910481 0.5436332

Realizando el modelo Koyck se puede inferir que la ecuacion del modelo tiene forma : Y =-0.005184 + 0.543633Y.1 - 0.491048X.t, estos coeficientes dados por el modelo indican que un aumento de una unidad en Y.1 se asocia con un aumento de 0.543633 unidades en la tasa de interés institucional, mientras que un aumento de una unidad en X.t se asocia con una disminución de 0.491048 unidades en la tasa de interés institucional o FEDFUND.

koyck2 <- lm( dFEDFUND~ lag(dFEDFUND) + dT10Y2YM,
             data=subdata)
summary(koyck2)

## 
## Call:
## lm(formula = dFEDFUND ~ lag(dFEDFUND) + dT10Y2YM, data = subdata)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.84082 -0.06013  0.00451  0.05982  0.43717 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -0.004848   0.006737  -0.720    0.472    
## lag(dFEDFUND)  0.559476   0.037223  15.030  < 2e-16 ***
## dT10Y2YM      -0.404452   0.048548  -8.331 1.36e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1336 on 392 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.4899, Adjusted R-squared:  0.4873 
## F-statistic: 188.3 on 2 and 392 DF,  p-value: < 2.2e-16

La variable FEFUND tiene un coeficiente de 0.5595. Este coeficiente es estadísticamente significativo y positivo, lo que indica que un aumento en la variable independiente se asocia con un aumento en la variable dependiente en el período siguiente.Por otra parte, la variable “dT10Y2YM” tiene un coeficiente de -0.4044. Este coeficiente es estadísticamente significativo y negativo, lo que indica que un aumento en la variable independiente se asocia con una disminución en la variable dependiente.

h <- durbinH(koyck2, "lag(dFEDFUND)")
abs(h)

## [1] 3.961956

2*pnorm(-abs(h)) # Calculando el pvalor

## [1] 7.433819e-05

Calculando el h de Durbin en un modelo autorregresivo, se puede afirmar que la hipótesis nula del modelo es H0: No autocorrelación, rechazo si |h|> 1.96 o si pvalor<5%. En este caso h es de 3.96 lo que quiere decir que la hipótesis nula se rechaza y el modelo puede estar presentando autocorrelacion.

medianar <- -(log(2)/log(0.5178846))
medianar

## [1] 1.053411

Inferimos con los resultados dados que el 50% del cambio total de las tasas de interés interbancarias se logra en más de un periodo.

mediar <- 0.5178846/(1-0.5178846)
mediar

## [1] 1.074192

Se requiere un promedio de mas de un periodo para que el efecto de los cambios de los bonos de EEUU afecte la tasa interbancaria.

Modelo ADL

Comparación modelos ADL

Para hacer el modelo ADL, se pondrá en comparación varias estimaciones para así lograr escoger el modelo ADL más optimo para realizar predicciones.

# ADL (1,1)
ADL11 <- ardlDlm(formula = dFEDFUND ~ dT10Y2YM, 
                 data = data.frame(subdata), p = 1, q = 1)
summary(ADL11)

## 
## Time series regression with "ts" data:
## Start = 3, End = 397
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.84345 -0.05982  0.00348  0.05878  0.43712 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.004959   0.006749  -0.735    0.463    
## dT10Y2YM.t  -0.399274   0.050030  -7.981 1.63e-14 ***
## dT10Y2YM.1  -0.023616   0.054205  -0.436    0.663    
## dFEDFUND.1   0.552393   0.040653  13.588  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1338 on 391 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.4902, Adjusted R-squared:  0.4863 
## F-statistic: 125.3 on 3 and 391 DF,  p-value: < 2.2e-16

#ADL (1,2)
ADL12 <- ardlDlm(formula = dFEDFUND ~ dT10Y2YM, 
                 data = data.frame(subdata), p = 1, q = 2)
summary(ADL12)

## 
## Time series regression with "ts" data:
## Start = 4, End = 397
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.83860 -0.06036  0.00122  0.05663  0.46521 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.004094   0.006627  -0.618    0.537    
## dT10Y2YM.t  -0.409656   0.049162  -8.333 1.37e-15 ***
## dT10Y2YM.1  -0.033590   0.053192  -0.631    0.528    
## dFEDFUND.1   0.425101   0.049951   8.510 3.81e-16 ***
## dFEDFUND.2   0.194416   0.046001   4.226 2.96e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1311 on 389 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.5125, Adjusted R-squared:  0.5075 
## F-statistic: 102.2 on 4 and 389 DF,  p-value: < 2.2e-16

#ADL (2,1)
ADL21 <- ardlDlm(formula = dFEDFUND ~ dT10Y2YM, 
                 data = data.frame(subdata), p = 2, q = 1)
summary(ADL21)

## 
## Time series regression with "ts" data:
## Start = 4, End = 397
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.83682 -0.06058  0.00187  0.05987  0.43826 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.005463   0.006715  -0.814  0.41637    
## dT10Y2YM.t  -0.416939   0.050170  -8.311 1.61e-15 ***
## dT10Y2YM.1   0.007520   0.055064   0.137  0.89144    
## dT10Y2YM.2  -0.138887   0.051503  -2.697  0.00731 ** 
## dFEDFUND.1   0.521715   0.041952  12.436  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1329 on 389 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.4995, Adjusted R-squared:  0.4943 
## F-statistic: 97.04 on 4 and 389 DF,  p-value: < 2.2e-16

#ADL (2,2)
ADL22 <- ardlDlm(formula = dFEDFUND ~ dT10Y2YM, 
                 data = data.frame(subdata), p = 2, q = 2)
summary(ADL22)

## 
## Time series regression with "ts" data:
## Start = 4, End = 397
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.83539 -0.05749  0.00000  0.05754  0.46298 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.004501   0.006622  -0.680 0.497113    
## dT10Y2YM.t  -0.418386   0.049439  -8.463 5.41e-16 ***
## dT10Y2YM.1  -0.014700   0.054619  -0.269 0.787968    
## dT10Y2YM.2  -0.079210   0.053458  -1.482 0.139225    
## dFEDFUND.1   0.422358   0.049909   8.463 5.42e-16 ***
## dFEDFUND.2   0.171892   0.048380   3.553 0.000428 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1309 on 388 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.5152, Adjusted R-squared:  0.509 
## F-statistic: 82.47 on 5 and 388 DF,  p-value: < 2.2e-16

Luego de tener los diferentes modelos ADL se pretende entrar a hacer comparaciones con cada uno de los modelos y de alli escoger el mejor.

AIC(ADL11); AIC(ADL12); AIC(ADL21); AIC(ADL22)

## [1] -462.2412

## [1] -475.752

## [1] -465.3606

## [1] -475.9752

BIC (ADL11); BIC(ADL12);BIC(ADL21); BIC(ADL22)

## [1] -442.3468

## [1] -451.8939

## [1] -441.5025

## [1] -448.1407

De los cuatro modelos, ADL11 tiene el valor AIC más bajo (-462.2412), lo que indica que es el mejor modelo de los cuatro. Los valores AIC de los otros modelos son todos más altos, lo que sugiere que son peores modelos que ADL11. Lo mismo pasa para BIC, el modelo con BIC más bajo es el ADL(11), por lo que en ambos casos el modelo ADL(11).

Modelo ADL(11)

ADL11 <- ardlDlm(formula = dFEDFUND ~ dT10Y2YM, 
                 data = subdata, p = 1, q = 1)
summary(ADL11)

## 
## Time series regression with "ts" data:
## Start = 3, End = 397
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.84345 -0.05982  0.00348  0.05878  0.43712 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.004959   0.006749  -0.735    0.463    
## dT10Y2YM.t  -0.399274   0.050030  -7.981 1.63e-14 ***
## dT10Y2YM.1  -0.023616   0.054205  -0.436    0.663    
## dFEDFUND.1   0.552393   0.040653  13.588  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1338 on 391 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.4902, Adjusted R-squared:  0.4863 
## F-statistic: 125.3 on 3 and 391 DF,  p-value: < 2.2e-16

Analizando los resultados de la regresión se puede notar que los valores de p para todos los coeficientes fueron menor que 0.05, lo que indica que son estadísticamente significativos.

predict.AD11 <- dLagM::forecast(ADL11, x = c(-0.06, -0.21, -0.06, -0.01), h = 4, interval = T)
predict.AD11

## $forecasts
##       95% LB   Forecast    95% UB
## 1 -0.1365257 0.14628433 0.3709256
## 2 -0.1494162 0.16111235 0.4424667
## 3 -0.1971214 0.11295455 0.4104797
## 4 -0.2516804 0.06284631 0.3816475
## 
## $call
## forecast.ardlDlm(model = ADL11, x = c(-0.06, -0.21, -0.06, -0.01), 
##     h = 4, interval = T)
## 
## attr(,"class")
## [1] "forecast.ardlDlm" "dLagM"

Para cada uno de los cuatro períodos futuros, se presentan tres valores: el límite inferior del intervalo de confianza del 95%, el pronóstico y el límite superior del intervalo de confianza del 95%. Por ejemplo, para el primer período futuro, el límite inferior del intervalo de confianza del 95% es -0.1110775, el pronóstico es 0.14628433 y el límite superior del intervalo de confianza del 95% es 0.3956422.

error.predict <- predict.AD11[["forecasts"]][["Forecast"]][1] - subdata$FEDFUND[394]
error.predict

## [1] -2.933716

Por otra parte, el error que se presenta en el modelo es de -2.9337. Esto significa que la predicción realizada subestimó el valor real de la variable en ese período en 2.933716 unidades.

Es con esto en mente que se procede a realizar las predicciones para los cuatro periodos determinados.

data_FF <- data.frame(obs1=subdata$FEDFUND[1:393]) 
NAs <- data.frame(obs1=rep(NA,4))
data_FF <- rbind(data_FF, NAs)

FF394 <- 0.14628433 + subdata$FEDFUND[393]
FF395 <- 0.16111235 + FF394
FF396 <- 0.11295455 + FF395
FF397 <-  0.06284631 + FF396

FF1 <- data.frame(pred1=c(subdata$FEDFUND[393], FF394, FF395, FF396, FF397))
Otras_celdas1 <- data.frame(pred1=rep(NA, 392))
FF1 <- rbind(Otras_celdas1, FF1)
data_FF <- data.frame(data_FF, FF1)

#--------------------------------------------------------------------------------
uFF394 <- 0.4221359+subdata$FEDFUND[393]
uFF395 <- 0.4613492+uFF394
uFF396 <- 0.4241844+uFF395
uFF397 <-  0.4066458+uFF396

uFF1 <- data.frame(u1=c(subdata$FEDFUND[393], uFF394, uFF395, uFF396, uFF397))

Otras_celdasu <- data.frame(u1=rep(NA, 392))

uFF1 <- rbind(Otras_celdasu, uFF1)
data_FF <- data.frame(data_FF, uFF1)
#----------------------------------------------------------------------------------------
lFF394 <- -0.1239314 + subdata$FEDFUND[393]
lFF395 <- -0.1297740 + lFF394
lFF396 <- -0.2097703 + lFF395
lFF397 <- -0.2720789 + lFF396

lFF1 <- data.frame(l1=c(subdata$FEDFUND[393], lFF394, lFF395, lFF396, lFF397))

Otras_celdasl <- data.frame(l1=rep(NA, 392))

lFF1 <- rbind(Otras_celdasl, lFF1)

data_FF <- data.frame(data_FF, lFF1)

data_FF <- data.frame(Date=subdata$year,data_FF)

Predicción

Después de haber hecho las predicciones la gráfica nos muestra que la tasa de interés interbancaria sigue con la tendencia a la alza en los periodos siguientes.

ggplot(data_FF)+
  geom_line(aes(x = Date, y = obs1, group =1, color = "Observada"), linewidth = 0.8)+
  geom_line(aes(x= Date, y= pred1, group =1, color="Predicha"), linewidth = 0.8)+
  geom_ribbon(aes(x =Date, y = pred1, group =1, ymin = l1, ymax = u1, fill="IC"), alpha = 0.1)+
  theme(legend.text = element_text(size = 10), text = element_text(size=8), legend.spacing.y = unit(-0.3, "cm"), legend.background=element_blank())+guides(shape = guide_legend(order = 2),col = guide_legend(order = 1)) + 
  scale_color_manual(name = "", values = c("Observada" = "darkblue", "Predicha" = "red")) + 
  scale_fill_manual(values = c(IC = "steelblue"), labels = c(IC = "IC 95%")) + 
  labs(x = "Fecha", y ="", fill = "") + ggtitle("FEDFUND")

Modelo ARMAX

Después de haber tenido el modelo ADL, se pretende realizar el modelo ARMAX, esto con el fin de realizar un paralelo entrelas predicciones dadas por ambos modelos. Para esto, con los correlogramas mostrados previamente podemos hacer un imaginario de los valores tomados por el ARIMA, es por esto que tendremos que hacer varias estimaciones de modelos.

subdata <- subdata |> 
  mutate(lagdT10 = lag(dT10Y2YM))

subdata3 <- subdata |>
  select(year, dFEDFUND, dT10Y2YM, lagdT10)

## ARMAX 2,1,1

ARMAX211 <- arima(subdata3$dFEDFUND, order = c(2,1,1), xreg = subdata3[,3:4])
summary(ARMAX211)

## 
## Call:
## arima(x = subdata3$dFEDFUND, order = c(2, 1, 1), xreg = subdata3[, 3:4])
## 
## Coefficients:
##          ar1     ar2      ma1  dT10Y2YM  lagdT10
##       0.4794  0.1772  -0.9954   -0.3349  -0.1028
## s.e.  0.0505  0.0507   0.0121    0.0481   0.0487
## 
## sigma^2 estimated as 0.01803:  log likelihood = 230.51,  aic = -449.02
## 
## Training set error measures:
##                       ME      RMSE        MAE MPE MAPE      MASE        ACF1
## Training set 0.008045327 0.1341091 0.08905719 NaN  Inf 0.9113365 -0.02133427

coeftest(ARMAX211)

## 
## z test of coefficients:
## 
##           Estimate Std. Error  z value  Pr(>|z|)    
## ar1       0.479355   0.050476   9.4966 < 2.2e-16 ***
## ar2       0.177241   0.050701   3.4958 0.0004726 ***
## ma1      -0.995397   0.012083 -82.3818 < 2.2e-16 ***
## dT10Y2YM -0.334926   0.048116  -6.9608 3.384e-12 ***
## lagdT10  -0.102764   0.048651  -2.1123 0.0346623 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# ARMAX 1,1,1 

ARMAX111 <- arima(subdata3$dFEDFUND, order = c(1,1,1), xreg = subdata3[,3:4])
summary(ARMAX111)

## 
## Call:
## arima(x = subdata3$dFEDFUND, order = c(1, 1, 1), xreg = subdata3[, 3:4])
## 
## Coefficients:
##          ar1      ma1  dT10Y2YM  lagdT10
##       0.3141  -0.8090   -0.3151  -0.1020
## s.e.  0.0850   0.0587    0.0494   0.0493
## 
## sigma^2 estimated as 0.01839:  log likelihood = 227.82,  aic = -445.64
## 
## Training set error measures:
##                       ME      RMSE        MAE MPE MAPE      MASE        ACF1
## Training set 0.002884091 0.1354477 0.09228673 NaN  Inf 0.9443849 -0.01530316

coeftest(ARMAX111)

## 
## z test of coefficients:
## 
##           Estimate Std. Error  z value  Pr(>|z|)    
## ar1       0.314105   0.085044   3.6934 0.0002213 ***
## ma1      -0.809013   0.058742 -13.7723 < 2.2e-16 ***
## dT10Y2YM -0.315121   0.049363  -6.3838 1.727e-10 ***
## lagdT10  -0.101955   0.049326  -2.0670 0.0387374 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Después de haber planteado los modelos ARMAX, tendremos que ponerlos a competir para verficar cual de los modelos planteados es el mejor para realizar la predicción.

# AIC Y BIC 

AIC(ARMAX211); AIC(ARMAX111)

## [1] -449.022

## [1] -445.6432

BIC(ARMAX211); BIC(ARMAX111)

## [1] -425.1639

## [1] -425.7614

El modelo ARMAX211 tiene una AIC de -451.0984, mientras que el modelo ARMAX111 tiene una AIC de -452.0493. El modelo con el AIC más bajo es generalmente considerado como el mejor modelo. En este caso, el modelo ARMAX211 tiene un AIC ligeramente más bajo, lo que indica que puede ser el mejor modelo.Entonces, en este caso, el modelo ARMAX211 tiene un BIC más bajo (-427.225) en comparación con ARMAX111 (-432.1549), lo que indica que ARMAX211 es un modelo que se ajusta mejor.

Procedemos a estimar las predicciones del modelo

newxreg <- data.frame(dT10Y2YM = c(-0.6,-0.21,-0.6,-0.1),
                      dT10 = c(0.01,-0.06, -0.21, -0.06))

predict.ARMAX <- predict(ARMAX111, n.ahead = 4, newxreg = newxreg)

upper95 <- predict.ARMAX[["pred"]] + 1.96*predict.ARMAX[["se"]]
lower95 <- predict.ARMAX[["pred"]] - 1.96*predict.ARMAX[["se"]]

df <- data.frame(observado = subdata$FEDFUND[1:393])
NAs0 <- data.frame(observado=rep(NA, 4))
df <- rbind(df, NAs0)

pred <- data.frame(predicho=c(subdata$FEDFUND[393], predict.ARMAX$pred))
NAs1 <- data.frame(predicho=rep(NA, 392))
pred <- rbind(NAs1, pred)
df <- data.frame(df, pred)

upper <- data.frame(upper=c(subdata$FEDFUND[393], upper95))
NAs2 <- data.frame(upper=rep(NA,392))
upper <- rbind(NAs2, upper) 
df <- data.frame(df, upper)


lower <- data.frame(lower=c(subdata$FEDFUND[393], lower95))
NAs3 <- data.frame(lower=rep(NA, 392)) 
lower <- rbind(NAs3, lower)
df <- data.frame(df, lower)

df <- data.frame(Date = subdata$year, df)

Predicción

A diferencia de la prediccion pasda, en esta se logra observar que hay una disminución importante en la tasa interbancaria en el rango elegido para hacer la predicción, es así que nos da una caída en la tendencia que luege parece estabilizarce un poco.

ggplot(df) + 
  geom_line(aes(x = Date, y = observado, group = 1,color = "Observada"), linewidth = 0.8) + 
  geom_line(aes(x = Date, y = predicho, group = 1, color = "Predicha"), linewidth = 0.8) + 
  geom_ribbon(aes(x = Date, y = predicho, group = 1, ymin = lower, ymax = upper, fill="IC"), alpha = 0.1)+
  theme(legend.text = element_text(size = 10), text = element_text(size=8), legend.spacing.y = unit(-0.3, "cm"), legend.background=element_blank()) + 
  guides(shape = guide_legend(order = 2),col = guide_legend(order = 1)) + 
  scale_color_manual(name = "", values = c("Observada" = "darkblue", "Predicha" = "red")) + 
  scale_fill_manual(values = c(IC = "steelblue"), labels = c(IC = "IC 95%")) + 
  labs(x = "Fecha", y ="", fill = "") + ggtitle("Predicción del ARMAX(1,1,1)")