La idea de la siguiente entrada es mostrar una estimación es timar un VAR donde se incluya un factor extraído mediante componentes principales como variable dependiente. Puntualmente, se usa la data EuStockMarkets, para estimar la variable DAX en función del primer componente identificado a partir del resto de cotizaciones disponibles en la base (SMI, CAC, FTSE).
library(dplyr)
library(vars)
library(forecast)
library(tidyverse)
EuStockMarkets %>% head()
Time Series:
Start = c(1991, 130)
End = c(1991, 135)
Frequency = 260
DAX SMI CAC
FTSE
1991.496 1628.75 1678.1 1772.8
2443.6
1991.500 1613.63 1688.5 1750.5
2460.2
1991.504 1606.51 1678.6 1718.0
2448.2
1991.508 1621.04 1684.1 1708.1
2470.4
1991.512 1618.16 1686.6 1723.1
2484.7
1991.515 1610.61 1671.6 1714.3
2466.8
Paso 1. Obtener series estacionarias.
En primer lugar, diferenciamos las series para tratar el tema de raíz unitaria. Este es un ejemplo no se trata este tema con mayor rigor, sin embargo, en ejercicios reales es necesario hacer el test.
as.tibble() %>%
mutate(across(
where(is.numeric), ~(log(.)-dplyr::lag(log(.)))*100)
) %>%
na.omit()
# A tibble: 1,859 x 4
DAX SMI CAC FTSE
<dbl> <dbl> <dbl> <dbl>
1 -0.933 0.618 -1.27 0.677
2 -0.442 -0.588 -1.87 -0.489
3 0.900 0.327 -0.578 0.903
4 -0.178 0.148 0.874 0.577
5 -0.468 -0.893 -0.512 -0.723
6 1.24 0.674 1.17 0.852
7 0.576 1.22 1.31 0.821
8 -0.287 -0.359 -0.194 0.0837
9 0.635 1.10 0.0171 -0.523
10 0.118 0.436 0.313 1.40
# ... with 1,849 more rows
Paso 2. Obtener el componente principal del conjunto de variables de interés.
Ahora obtengo el primer componente de los retornos estandarizados. Recordando que solo corresponde a las tres últimas series. Obtengo los ponderadores correspondientes y recupero el primer componente principal, luego de revisar el porcentaje de la varianza explicada por cada uno de los componentes. De la función prcomp podemos extraer directamente el primer componente asociado. El mismo, explica el 74.4% de la variabilidad del vector de variable de interés.
summary(componteP)
Importance of components:
PC1 PC2 PC3
Standard deviation 1.4944 0.6482 0.5887
Proportion of Variance 0.7444 0.1401 0.1155
Cumulative Proportion 0.7444 0.8845 1.0000
# extract principal component
CP1 <- componteP$x[,1]
Paso 3. Estimo el FAVAR
Aquí estimamos el FAVAR aunque omitido las ultimas 5 observaciones, esto para poder comparar las proyecciones con los valores observados posteriormente. En un ejercicio mensual de pronósticos no es necesario hacer este paso, solo se hace cuando se pretende evaluar la efectividad histórica de un modelo determinado.
data_var <- head(data_var, -6)
CP1
[1,] -0.9326550 0.001671168
[2,] -0.4422175 1.716175165
[3,] 0.9003794 -0.367691926
[4,] -0.1778217 -0.930842533
[5,] -0.4676712 1.224903643
[6,] 1.2427042 -1.561823450
akaikelag <- VARselect(data_var, lag.max = 5, type = "const")$selection[1]
var_model <- VAR(data_var, lag=akaikelag)
Paso alternativo. Proyecciones alternativas con representación gráfica.
Finalmente, se puede usar la estimación para obtener proyecciones. Aquí proyectamos como ejemplo 6 periodos y adherimos distintos modelos alternativos para comparar esas proyecciones, correspondientes a las ultimas 6 observaciones (recuerde que anteriormente omitimos estas últimas 6 observaciones los datos usados para estimar el modelo).
# modelos alternativos
vecm_model <- VECM(data_var,
lag=akaikelag) #modelo VECM
tvar_model <- TVAR(data_var,
lag=akaikelag, nthresh=1, thDelay=1) #modelo TVAR
ar_model <-
arima(data_var[,1], order = c(3,0,1))
ets_model <-
ets(data_var[,1])
# pronósticos
h <- 6
var.pred <-
predict(var_model, n.ahead = h)
vecm.pred <-
predict(vecm_model, n.ahead = h)
tvar.pred <-
predict(tvar_model, n.ahead = h)
p_Ar <- predict(ar_model,
n.ahead = h)$pred
p_ets <-
forecast::forecast(ets_model, h=h)$mean
Un
gráfico para comparar los pronósticos.
## compare forecasts on a
plot
n <- 30
plot(1:n,
tail(data_returns$DAX, n), type="l")
lines((n-h+1):n,
var.pred$fcst$X[,1], lty=2, col=2)
lines((n-h+1):n,
vecm.pred[,1], lty=2, col=3)
lines((n-h+1):n,
tvar.pred[,1], lty=2, col=4)
lines((n-h+1):n, p_Ar, lty=2,
col=5)
lines((n-h+1):n, p_ets,
lty=2, col=6)
legend("bottomleft",
lty=c(1,rep(2,5)), col=1:6,
legend=c("Obs.",
"var", "vecm", "tvar", "Arima",
"Ets"))
Paso alternativo. Recuperar los niveles de la proyección en tasas de variación.
Algo importante a notar es que hemos hecho el pronóstico de una tasa de variación, sin embargo, es más fácil comunicar niveles de la variable. Para recuperar los niveles usamos el vector de valores proyectados y el último valor observado:
last_obs <-
tail(EuStockMarkets[,"DAX"],6)[1] #el primero de las ultimas 6 obs
proy <-
var.pred$fcst$X[,1]
v_hat <-
last_obs*cumprod(1+proy/100)
[1] 5601.966 5605.684
5609.410 5613.140 5616.872 5620.606
plot(1:n,
tail(EuStockMarkets[,"DAX"], n), type="l")
lines((n-h+1):n, v_hat, lty=2,
col=2)
