Introduccion

Pasos para el modelado de una serie temporal:

• Dibujar gráficamente la serie.
• Transformación de la serie (la más común son los logaritmos) para hacerla mas estacional, sobre todo en varianza.
• Aplicar la función VAR_select para saber el orden del modelo VAR que minimiza los errores.
• Aplicar la función VAR para hacer el modelo de predicción.
• Comprobación de la aleatoriedad de los residuos con la función serial.test.
• Usar forecast para estimar valores futuros según el modelo.

Carga y gráficos de la serie

Para esta aproximación al modelado de series temporales multivariables con VAR, usaremos la serie Canada, que contiene cuatro indicadores económicos trimestrales de Canada desde 1980 al año 2001:

• prod - indice de producción
• e - nivel de empleo
• U - nivel de desempleo en %
• rw - salario real

help(Canada)

Vamos a cargar la serie temporal y a hacer un grafico para ver si aspecto:

data(Canada)
ts <- Canada
plot.ts(ts)

Vemos que no se aprecia estacionalidad a simple vista. Si que se aprecia tendencia en los cuatro indicadores.

Transformación de la serie y diferencias

Vamos a transformar la serie con logaritmos, para estabilizar la varianza y poder aplicar mejor un modelo.

ts.log <- log(ts)
plot.ts(ts.log)

Test de correlación

Vamos a calcular la autocorrelación de las series, asi como la correlación cruzada, para comprobar si las series están correladas entre si, de tal manera que podemos establecer modelos predictivos que saquen beneficio de esas correlaciones para hacer predicciones.

cor(ts) #matriz de correlación de la serie original

##               e       prod          rw           U
## e     1.0000000  0.8724332  0.92730670 -0.41975053
## prod  0.8724332  1.0000000  0.74783386 -0.49722710
## rw    0.9273067  0.7478339  1.00000000 -0.05856635
## U    -0.4197505 -0.4972271 -0.05856635  1.00000000

cor(ts.log) #matriz de correlación de la serie en forma logaritmica

##               e       prod          rw           U
## e     1.0000000  0.8712838  0.92351887 -0.40558682
## prod  0.8712838  1.0000000  0.73897848 -0.49372870
## rw    0.9235189  0.7389785  1.00000000 -0.03145907
## U    -0.4055868 -0.4937287 -0.03145907  1.00000000

#en ambas series se mantiene la correlación cruzada con varlores similares

Aplicación del modelo VAR

Vamos ahora a aplicar el modelado con VAR. Para ello, primero hemos de estimar el orden del modelo, es decir, cuantas muestras anteriores vamos a integrar en nuestro modelo. PAra hacer esto usamos la función VAR_select. Para la selección del orden usamos el SC (que es igual al BIC) más que el AIC (que se usa para series univariantes), ya que ofrece un orden menor, y por el principio de parsimonia elegimos los modelos más simples.

#type = c(conts, trend, both, none), para incluir constante y tendencia o ambos

lag_estimator <- VARselect(ts.log,type="both") #lag=1 según el SC (que es igual al BIC)
lag_estimator$selection

## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      3      2      1      3

lag <- lag_estimator$selection[3] #Seleccionamos según el BIC
lag

## SC(n) 
##     1

Una vez tenemos el orden (lag) del modelo, podemos usar la función VAR para crear el modelo de predicción.

A1 <- VAR(ts.log, p = lag, type = "both")
A1

## 
## VAR Estimation Results:
## ======================= 
## 
## Estimated coefficients for equation e: 
## ====================================== 
## Call:
## e = e.l1 + prod.l1 + rw.l1 + U.l1 + const + trend 
## 
##           e.l1        prod.l1          rw.l1           U.l1          const 
##  1.30773463474  0.09184686492 -0.03785043213  0.00754697858 -2.44360084556 
##          trend 
## -0.00005698038 
## 
## 
## Estimated coefficients for equation prod: 
## ========================================= 
## Call:
## prod = e.l1 + prod.l1 + rw.l1 + U.l1 + const + trend 
## 
##          e.l1       prod.l1         rw.l1          U.l1         const 
## -0.1625268434  0.9665460783 -0.0100501314  0.0009555908  1.3689922883 
##         trend 
##  0.0001117859 
## 
## 
## Estimated coefficients for equation rw: 
## ======================================= 
## Call:
## rw = e.l1 + prod.l1 + rw.l1 + U.l1 + const + trend 
## 
##          e.l1       prod.l1         rw.l1          U.l1         const 
## -0.1143637952 -0.2144032546  0.9476913702 -0.0081457822  2.4038290429 
##         trend 
##  0.0001642841 
## 
## 
## Estimated coefficients for equation U: 
## ====================================== 
## Call:
## U = e.l1 + prod.l1 + rw.l1 + U.l1 + const + trend 
## 
##          e.l1       prod.l1         rw.l1          U.l1         const 
## -33.887455143  -5.889314344   4.091409461   0.218974500 244.154885942 
##         trend 
##   0.004845017

Ahora vamos a ver si los residuos son ruido blanco y no están correlados entre si. Aunque el p-value sale algo por encima de 0.05, lo damos por bueno.

serial.test(A1, lags.pt=10, type="PT.asymptotic")

## 
##  Portmanteau Test (asymptotic)
## 
## data:  Residuals of VAR object A1
## Chi-squared = 167.1, df = 144, p-value = 0.09121

Ahora hacemos un grafico de la predicción de cada una de las series temporales.

forecast(A1) %>%
  autoplot() + xlab("Year")

Otro modelo

En el modelo anterior hemos visto que los coeficientes de tendencia son muy bajos. Vamos a hacer un segundo modelo sin tener el cuenta en coeficiente de tendencia, para ver si mejoramos el test de los residuos.

#type = c(const, trend, both, none), para incluir constante y tendencia o ambos

lag_estimator_2 <- VARselect(Canada,type="const") #lag=1 según el SC (que es igual al BIC)

lag_2 <- lag_estimator_2$selection[3] #Seleccionamos según el BIC
lag_2

## SC(n) 
##     1

A2 <- VAR(Canada, p = lag_2, type = "const")
A2

## 
## VAR Estimation Results:
## ======================= 
## 
## Estimated coefficients for equation e: 
## ====================================== 
## Call:
## e = e.l1 + prod.l1 + rw.l1 + U.l1 + const 
## 
##          e.l1       prod.l1         rw.l1          U.l1         const 
##    1.17353629    0.14479389   -0.07904568    0.52438144 -192.56360758 
## 
## 
## Estimated coefficients for equation prod: 
## ========================================= 
## Call:
## prod = e.l1 + prod.l1 + rw.l1 + U.l1 + const 
## 
##         e.l1      prod.l1        rw.l1         U.l1        const 
##   0.08709510   1.01970070  -0.02629309   0.32299246 -81.55109611 
## 
## 
## Estimated coefficients for equation rw: 
## ======================================= 
## Call:
## rw = e.l1 + prod.l1 + rw.l1 + U.l1 + const 
## 
##        e.l1     prod.l1       rw.l1        U.l1       const 
##  0.06381103 -0.13551199  0.96872851 -0.19538479 11.61375726 
## 
## 
## Estimated coefficients for equation U: 
## ====================================== 
## Call:
## U = e.l1 + prod.l1 + rw.l1 + U.l1 + const 
## 
##         e.l1      prod.l1        rw.l1         U.l1        const 
##  -0.19293575  -0.08086896   0.07538624   0.47530976 186.80892410

#Comprobamos que los residuos son ruido blanco

serial.test(A2, lags.pt=10, type="PT.asymptotic")

## 
##  Portmanteau Test (asymptotic)
## 
## data:  Residuals of VAR object A2
## Chi-squared = 157.92, df = 144, p-value = 0.2022

forecast(A2) %>%
  autoplot() + xlab("Year")

Vemos qu en el test de los residuos, obtenemos un valor muy alto, y por tanto, los residuos no son ruido blanco. Por lo tanto, con quedaremos con el primer modelo.

forecast(A1) %>%
  autoplot() + xlab("Year")