SESIÓN 6. Análisis de Serie Temporal Multivariante en R

Autor/a

Gerson Rivera

Fecha de publicación

29 agosto 2024

Etapas del Análisis TS Multivariante

  1. Análisis exploratorios de los datos.
  2. Dividir la serie en conjuntos de entrenamiento y prueba.
  3. Prueba de estacionariedad.
  4. Transformar la serie de entrenamiento (si es necesario).
  5. Construir un modelo VAR sobre las series transformadas.
  6. Causalidad de Granger.
  7. Diagnóstico del modelo.
  8. Realizar pronósticos utilizando el modelo elegido.
  9. Transformación inversa del pronóstico a la escala original.
  10. Realizar una evaluación del pronóstico.

Detalles Introductorios

Generar un dataframe aleatorio

x = rnorm(100, 1)
y = rnorm(100, 30)
z = rnorm(100, 500)
xyz = data.frame(x, y, z)
class(xyz)
[1] "data.frame"

Convertirlo a mts

mymts = ts(xyz,
           frequency = 12,
           start = c(1940, 4))
class(mymts)
[1] "mts"    "ts"     "matrix" "array"

Herramientas exploratorias básicas

plot(mymts)
head(mymts)
               x        y        z
[1,] -0.09952983 29.37945 500.4254
[2,]  0.28179870 27.77379 500.6494
[3,]  1.84926994 31.13288 501.8384
[4,]  0.42153761 30.40711 501.4630
[5,]  3.35560524 30.56498 499.8612
[6,]  1.57449048 28.89361 499.9377
#class(mymts)
library(TSstudio)
ts_plot(mymts)

Caso de Estudio: Datos macroeconómicos

Cargando los datos reales

  • La información de la serie corresponde a datos macroeconómicos de Estados Unidos, contiene dos variables de series de tiempo sobre la producción (PIB - en inglés GDP) y el desempleo.
library(readxl)
dat <- read_excel("blanchQua.xlsx")
View(dat)
#library(readr)
#dat = read_csv("blanchQua.csv")
#View(dat)
class(dat)
[1] "tbl_df"     "tbl"        "data.frame"
head(dat)
# A tibble: 6 × 3
  dates       GDP         U
  <dbl>     <dbl>     <dbl>
1 1948.  85443952 -14374997
2 1948.  12100904 -67461753
3 1948.  -8505847 -24506939
4 1949. -20765859  78511467
5 1949. -14918681  19614029
6 1949. -14913989  27710244

Paquetes importantes

  • Se realiza instalación o activación de los principales paquetes en R que tienen el modelo VAR
  • Observación. En ambos la función se llama igual pero los parámetros son diferentes
#install.packages("vars")
#install.packages("MTS")
library(vars)
Cargando paquete requerido: MASS
Cargando paquete requerido: strucchange
Cargando paquete requerido: zoo

Adjuntando el paquete: 'zoo'
The following objects are masked from 'package:base':

    as.Date, as.Date.numeric
Cargando paquete requerido: sandwich
Cargando paquete requerido: urca
Cargando paquete requerido: lmtest
library(MTS)  

Adjuntando el paquete: 'MTS'
The following object is masked from 'package:vars':

    VAR

Librería para el test ADF de estacionariedad

library(tseries)
Registered S3 method overwritten by 'quantmod':
  method            from
  as.zoo.data.frame zoo 
head(dat)
# A tibble: 6 × 3
  dates       GDP         U
  <dbl>     <dbl>     <dbl>
1 1948.  85443952 -14374997
2 1948.  12100904 -67461753
3 1948.  -8505847 -24506939
4 1949. -20765859  78511467
5 1949. -14918681  19614029
6 1949. -14913989  27710244

Análisis exploratorio

Convertir a objeto ts

gdp <- ts(dat$GDP, start = c(1948, 2), freq = 4)
une <- ts(dat$U, start = c(1948, 2), freq = 4)

Gráfico con plot

dat.bv <- cbind(gdp, une)
head(dat.bv)
           gdp       une
[1,]  85443952 -14374997
[2,]  12100904 -67461753
[3,]  -8505847 -24506939
[4,] -20765859  78511467
[5,] -14918681  19614029
[6,] -14913989  27710244
ts_plot(dat.bv)
library(lmtest)
grangertest(gdp,une)
Granger causality test

Model 1: une ~ Lags(une, 1:1) + Lags(gdp, 1:1)
Model 2: une ~ Lags(une, 1:1)
  Res.Df Df      F  Pr(>F)  
1    155                    
2    156 -1 2.8316 0.09444 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
grangertest(une,gdp)
Granger causality test

Model 1: gdp ~ Lags(gdp, 1:1) + Lags(une, 1:1)
Model 2: gdp ~ Lags(gdp, 1:1)
  Res.Df Df      F Pr(>F)  
1    155                   
2    156 -1 5.5805 0.0194 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#library(urca)
#jo=ca.jo(dat.bv,type='eigen')
#summary(jo)

Gráfico con autoplot

library(ggplot2)
#install.packages("ggfortify")
library(ggfortify)
autoplot(dat.bv )

plot.ts(dat.bv)
plot(dat.bv)

Dividir la serie (entrenamiento y prueba)

library(dplyr)

Adjuntando el paquete: 'dplyr'
The following object is masked from 'package:MASS':

    select
The following objects are masked from 'package:stats':

    filter, lag
The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union
dim(dat.bv)
[1] 159   2
n_obs=10
end=dim(dat.bv)[1]

X_train = dat.bv[1:(end-n_obs),]
X_test = dat.bv[(end-n_obs+1):end,]
dim(X_test)
[1] 10  2
dim(X_train)
[1] 149   2

Prueba de estacionariedad

apply(X_train, 2, adf.test) # Estamos en columnas (2)
Warning in FUN(newX[, i], ...): p-value smaller than printed p-value
$gdp

    Augmented Dickey-Fuller Test

data:  newX[, i]
Dickey-Fuller = -4.8587, Lag order = 5, p-value = 0.01
alternative hypothesis: stationary


$une

    Augmented Dickey-Fuller Test

data:  newX[, i]
Dickey-Fuller = -4.0081, Lag order = 5, p-value = 0.01078
alternative hypothesis: stationary

\underline{Nota}. El valor 2 es para especificar que lo queremos aplicar por columnas

  • Todos los p-valores son < 0.05, no hay que diferenciar las series, puesto que cumplen con estacionariedad.

VAR modeling

Identificación del orden del modelo

library(vars)

VARselect(X_train, type = "none", lag.max = 10)
$selection
AIC(n)  HQ(n)  SC(n) FPE(n) 
     1      1      1      1 

$criteria
                  1            2            3            4            5
AIC(n) 6.950704e+01 6.954474e+01 6.956378e+01 6.956397e+01 6.959292e+01
HQ(n)  6.954136e+01 6.961337e+01 6.966673e+01 6.970124e+01 6.976450e+01
SC(n)  6.959149e+01 6.971363e+01 6.981712e+01 6.990175e+01 7.001514e+01
FPE(n) 1.536482e+30 1.595546e+30 1.626347e+30 1.626892e+30 1.675080e+30
                  6            7            8            9           10
AIC(n) 6.961073e+01 6.965910e+01 6.969151e+01 6.966769e+01 6.968798e+01
HQ(n)  6.981663e+01 6.989931e+01 6.996604e+01 6.997653e+01 7.003114e+01
SC(n)  7.011740e+01 7.025022e+01 7.036707e+01 7.042769e+01 7.053243e+01
FPE(n) 1.705811e+30 1.791258e+30 1.851523e+30 1.809515e+30 1.848632e+30
library(MTS)

#library(urca)
#co=ca.jo(X_train)
#summary(jo)

#MTS::VARorder(X_train,12)

Creando el modelo

var.a <- vars::VAR(X_train,                   
                   lag.max = 10,                   
                   ic = "AIC",                   
                   type = "const")

summary(var.a)

VAR Estimation Results:
========================= 
Endogenous variables: gdp, une 
Deterministic variables: const 
Sample size: 148 
Log Likelihood: -5557.256 
Roots of the characteristic polynomial:
0.5536 0.1061
Call:
vars::VAR(y = X_train, type = "const", lag.max = 10, ic = "AIC")


Estimation results for equation gdp: 
==================================== 
gdp = gdp.l1 + une.l1 + const 

        Estimate Std. Error t value Pr(>|t|)   
gdp.l1 1.591e-01  7.829e-02   2.032  0.04396 * 
une.l1 1.947e-01  6.725e-02   2.895  0.00437 **
const  2.869e+06  2.814e+06   1.020  0.30963   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Residual standard error: 33870000 on 145 degrees of freedom
Multiple R-Squared: 0.08384,    Adjusted R-squared: 0.07121 
F-statistic: 6.635 on 2 and 145 DF,  p-value: 0.001749 


Estimation results for equation une: 
==================================== 
une = gdp.l1 + une.l1 + const 

         Estimate Std. Error t value Pr(>|t|)    
gdp.l1  1.074e-01  8.308e-02   1.293    0.198    
une.l1  5.006e-01  7.137e-02   7.014 8.17e-11 ***
const  -2.707e+06  2.986e+06  -0.906    0.366    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Residual standard error: 35940000 on 145 degrees of freedom
Multiple R-Squared: 0.2645, Adjusted R-squared: 0.2544 
F-statistic: 26.08 on 2 and 145 DF,  p-value: 2.12e-10 



Covariance matrix of residuals:
           gdp        une
gdp  1.147e+15 -1.035e+14
une -1.035e+14  1.292e+15

Correlation matrix of residuals:
         gdp      une
gdp  1.00000 -0.08507
une -0.08507  1.00000
#var.b <- MTS::VAR(X_train,p=2)
#summary(var.b)
#var.c <- MTS::VAR(X_train,p=1)
#summary(var.c)

Causalidad de Granger

library(vars)
causality(var.a, cause = c("gdp"))
$Granger

    Granger causality H0: gdp do not Granger-cause une

data:  VAR object var.a
F-Test = 1.6722, df1 = 1, df2 = 290, p-value = 0.197


$Instant

    H0: No instantaneous causality between: gdp and une

data:  VAR object var.a
Chi-squared = 1.0635, df = 1, p-value = 0.3024
causality(var.a, cause = c("une"))
$Granger

    Granger causality H0: une do not Granger-cause gdp

data:  VAR object var.a
F-Test = 8.3825, df1 = 1, df2 = 290, p-value = 0.004076


$Instant

    H0: No instantaneous causality between: une and gdp

data:  VAR object var.a
Chi-squared = 1.0635, df = 1, p-value = 0.3024

Diagnosis del modelo

  • Se utiliza Portmanteau test para objetos var
bv.serial= serial.test(var.a)
bv.serial

    Portmanteau Test (asymptotic)

data:  Residuals of VAR object var.a
Chi-squared = 65.502, df = 60, p-value = 0.2918
  • Deseamos obtener un p\_valor > 0.05

Posibles soluciones si p\_valor< 0.05

  • Cambiar el orden del modelo.
  • Cambiar el tipo de modelo.
  • Añadir otro paso de diferenciación.
  • Transformar con logaritmos.
plot(bv.serial, names = "gdp")

plot(bv.serial, names = "une")

Forecasting usando el modelo VAR

predictions <- predict(var.a, n.ahead = 10, ci = 0.95)
plot(predictions, names = "gdp")

predictions <- predict(var.a, n.ahead = 10, ci = 0.95)
plot(predictions, names = "une")

Alternativa gráfica

fanchart(predictions, names = "gdp")

fanchart(predictions, names = "une")

predictions$fcst$gdp[,1]
 [1] -4209701.1  -678481.2   705577.2  1411239.8  1795470.2  2007497.7
 [7]  2124803.4  2189735.8  2225681.4  2245580.6

Evaluando el modelo

pred=predictions$fcst
rmse=sqrt(mean((X_test[,1]-pred$gdp)^2))
cat('RMSE gdp: ', rmse)
RMSE gdp:  72015608
rmse=sqrt(mean((X_test[,2]-pred$une)^2))
cat('RMSE une: ', rmse)
RMSE une:  86200366

Entrenando el modelo con todos los datos

VAR modeling

Identificación del orden del modelo

VARselect(dat.bv, type = "none", lag.max = 10)
$selection
AIC(n)  HQ(n)  SC(n) FPE(n) 
     1      1      1      1 

$criteria
                  1            2            3            4            5
AIC(n) 6.954718e+01 6.958299e+01 6.959022e+01 6.958581e+01 6.960660e+01
HQ(n)  6.957994e+01 6.964852e+01 6.968852e+01 6.971687e+01 6.977042e+01
SC(n)  6.962782e+01 6.974428e+01 6.983215e+01 6.990838e+01 7.000982e+01
FPE(n) 1.599395e+30 1.657757e+30 1.669890e+30 1.662739e+30 1.698002e+30
                  6            7            8            9           10
AIC(n) 6.962425e+01 6.966872e+01 6.970496e+01 6.968073e+01 6.970367e+01
HQ(n)  6.982084e+01 6.989806e+01 6.996707e+01 6.997560e+01 7.003131e+01
SC(n)  7.010811e+01 7.023322e+01 7.035011e+01 7.040651e+01 7.051010e+01
FPE(n) 1.728753e+30 1.808098e+30 1.875871e+30 1.832253e+30 1.876449e+30

Creando el modelo

var.a <- vars::VAR(dat.bv,                 
                   lag.max = 10,                   
                   ic = "AIC",                   
                   type = "const")

summary(var.a)

VAR Estimation Results:
========================= 
Endogenous variables: gdp, une 
Deterministic variables: const 
Sample size: 158 
Log Likelihood: -5935.824 
Roots of the characteristic polynomial:
0.5574 0.1062
Call:
vars::VAR(y = dat.bv, type = "const", lag.max = 10, ic = "AIC")


Estimation results for equation gdp: 
==================================== 
gdp = gdp.l1 + une.l1 + const 

        Estimate Std. Error t value Pr(>|t|)  
gdp.l1 1.599e-01  7.825e-02   2.044   0.0427 *
une.l1 1.601e-01  6.778e-02   2.362   0.0194 *
const  3.796e+06  2.830e+06   1.341   0.1818  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Residual standard error: 3.5e+07 on 155 degrees of freedom
Multiple R-Squared: 0.06199,    Adjusted R-squared: 0.04988 
F-statistic: 5.121 on 2 and 155 DF,  p-value: 0.007018 


Estimation results for equation une: 
==================================== 
une = gdp.l1 + une.l1 + const 

         Estimate Std. Error t value Pr(>|t|)    
gdp.l1  1.333e-01  7.919e-02   1.683   0.0944 .  
une.l1  5.037e-01  6.859e-02   7.344  1.1e-11 ***
const  -3.595e+06  2.864e+06  -1.255   0.2113    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Residual standard error: 35410000 on 155 degrees of freedom
Multiple R-Squared: 0.2726, Adjusted R-squared: 0.2632 
F-statistic: 29.04 on 2 and 155 DF,  p-value: 1.943e-11 



Covariance matrix of residuals:
           gdp        une
gdp  1.225e+15 -1.047e+14
une -1.047e+14  1.254e+15

Correlation matrix of residuals:
         gdp      une
gdp  1.00000 -0.08445
une -0.08445  1.00000

Causalidad de Granger

causality(var.a, cause = c("gdp"))
$Granger

    Granger causality H0: gdp do not Granger-cause une

data:  VAR object var.a
F-Test = 2.8316, df1 = 1, df2 = 310, p-value = 0.09343


$Instant

    H0: No instantaneous causality between: gdp and une

data:  VAR object var.a
Chi-squared = 1.1189, df = 1, p-value = 0.2902
causality(var.a, cause = c("une"))
$Granger

    Granger causality H0: une do not Granger-cause gdp

data:  VAR object var.a
F-Test = 5.5805, df1 = 1, df2 = 310, p-value = 0.01878


$Instant

    H0: No instantaneous causality between: une and gdp

data:  VAR object var.a
Chi-squared = 1.1189, df = 1, p-value = 0.2902

Diagnosis del modelo

  • Se utiliza Portmanteau test para objetos var
bv.serial= serial.test(var.a)
bv.serial

    Portmanteau Test (asymptotic)

data:  Residuals of VAR object var.a
Chi-squared = 63.601, df = 60, p-value = 0.3509
  • Deseamos obtener un p\_valor > 0.05

Posibles soluciones si p\_valor< 0.05

  • Cambiar el orden del modelo.
  • Cambiar el tipo de modelo.
  • Añadir otro paso de diferenciación.
  • Transformar con logaritmos.
plot(bv.serial, names = "gdp")

plot(bv.serial, names = "une")

Forecasting usando el modelo VAR

predictions <- predict(var.a, n.ahead = 15, ci = 0.95)
plot(predictions, names = "gdp")

predictions <- predict(var.a, n.ahead = 15, ci = 0.95)
plot(predictions, names = "une")

Alternativa gráfica

fanchart(predictions, names = "gdp")

fanchart(predictions, names = "une")

Más en el futuro “n.ahead = 50” (perdemos precisión)

predictions <- predict(var.a, n.ahead = 50, ci = 0.95)
plot(predictions, names = "gdp")

predictions <- predict(var.a, n.ahead = 50, ci = 0.95)
plot(predictions, names = "une")

Otro gráfico

fanchart(predictions, names = "gdp") 

fanchart(predictions, names = "une") 

predictions$fcst 
$gdp
          fcst     lower    upper       CI
 [1,] 16448644 -52143536 85040824 68592180
 [2,]  6516941 -63683817 76717700 70200758
 [3,]  4659055 -65973064 75291175 70632119
 [4,]  4014225 -66752940 74781391 70767165
 [5,]  3696308 -67112914 74505531 70809223
 [6,]  3523510 -67298785 74345804 70822294
 [7,]  3427659 -67398697 74254015 70826356
 [8,]  3374281 -67453337 74201899 70827618
 [9,]  3344533 -67483477 74172543 70828010
[10,]  3327952 -67500180 74156083 70828132
[11,]  3318709 -67509461 74146879 70828170
[12,]  3313557 -67514624 74141739 70828181
[13,]  3310686 -67517499 74138871 70828185
[14,]  3309085 -67519101 74137271 70828186
[15,]  3308193 -67519994 74136379 70828187
[16,]  3307695 -67520491 74135882 70828187
[17,]  3307418 -67520769 74135605 70828187
[18,]  3307264 -67520923 74135450 70828187
[19,]  3307177 -67521009 74135364 70828187
[20,]  3307129 -67521057 74135316 70828187
[21,]  3307103 -67521084 74135289 70828187
[22,]  3307088 -67521099 74135275 70828187
[23,]  3307079 -67521107 74135266 70828187
[24,]  3307075 -67521112 74135262 70828187
[25,]  3307072 -67521114 74135259 70828187
[26,]  3307071 -67521116 74135258 70828187
[27,]  3307070 -67521117 74135257 70828187
[28,]  3307070 -67521117 74135256 70828187
[29,]  3307069 -67521117 74135256 70828187
[30,]  3307069 -67521118 74135256 70828187
[31,]  3307069 -67521118 74135256 70828187
[32,]  3307069 -67521118 74135256 70828187
[33,]  3307069 -67521118 74135256 70828187
[34,]  3307069 -67521118 74135256 70828187
[35,]  3307069 -67521118 74135256 70828187
[36,]  3307069 -67521118 74135256 70828187
[37,]  3307069 -67521118 74135256 70828187
[38,]  3307069 -67521118 74135256 70828187
[39,]  3307069 -67521118 74135256 70828187
[40,]  3307069 -67521118 74135256 70828187
[41,]  3307069 -67521118 74135256 70828187
[42,]  3307069 -67521118 74135256 70828187
[43,]  3307069 -67521118 74135256 70828187
[44,]  3307069 -67521118 74135256 70828187
[45,]  3307069 -67521118 74135256 70828187
[46,]  3307069 -67521118 74135256 70828187
[47,]  3307069 -67521118 74135256 70828187
[48,]  3307069 -67521118 74135256 70828187
[49,]  3307069 -67521118 74135256 70828187
[50,]  3307069 -67521118 74135256 70828187

$une
            fcst     lower    upper       CI
 [1,]   564802.1 -68845302 69974906 69410104
 [2,] -1118941.2 -79027858 76789976 77908917
 [3,] -3290486.0 -83612245 77031273 80321759
 [4,] -4631919.1 -85686600 76422762 81054681
 [5,] -5393561.5 -86674497 75887374 81280936
 [6,] -5819585.4 -87170684 75531514 81351099
 [7,] -6057211.4 -87430098 75315675 81372886
 [8,] -6189682.7 -87569337 75189972 81379654
 [9,] -6263524.9 -87645282 75118232 81381757
[10,] -6304685.4 -87687096 75077725 81382410
[11,] -6327628.5 -87710242 75054985 81382613
[12,] -6340417.2 -87723094 75042259 81382676
[13,] -6347545.7 -87730242 75035150 81382696
[14,] -6351519.2 -87734221 75031183 81382702
[15,] -6353734.1 -87736438 75028970 81382704
[16,] -6354968.7 -87737673 75027736 81382705
[17,] -6355656.8 -87738362 75027048 81382705
[18,] -6356040.4 -87738745 75026664 81382705
[19,] -6356254.2 -87738959 75026451 81382705
[20,] -6356373.4 -87739078 75026331 81382705
[21,] -6356439.8 -87739145 75026265 81382705
[22,] -6356476.9 -87739182 75026228 81382705
[23,] -6356497.5 -87739202 75026207 81382705
[24,] -6356509.0 -87739214 75026196 81382705
[25,] -6356515.4 -87739220 75026189 81382705
[26,] -6356519.0 -87739224 75026186 81382705
[27,] -6356521.0 -87739226 75026184 81382705
[28,] -6356522.1 -87739227 75026183 81382705
[29,] -6356522.7 -87739228 75026182 81382705
[30,] -6356523.1 -87739228 75026182 81382705
[31,] -6356523.3 -87739228 75026182 81382705
[32,] -6356523.4 -87739228 75026181 81382705
[33,] -6356523.4 -87739228 75026181 81382705
[34,] -6356523.5 -87739228 75026181 81382705
[35,] -6356523.5 -87739228 75026181 81382705
[36,] -6356523.5 -87739228 75026181 81382705
[37,] -6356523.5 -87739228 75026181 81382705
[38,] -6356523.5 -87739228 75026181 81382705
[39,] -6356523.5 -87739228 75026181 81382705
[40,] -6356523.5 -87739228 75026181 81382705
[41,] -6356523.5 -87739228 75026181 81382705
[42,] -6356523.5 -87739228 75026181 81382705
[43,] -6356523.5 -87739228 75026181 81382705
[44,] -6356523.5 -87739228 75026181 81382705
[45,] -6356523.5 -87739228 75026181 81382705
[46,] -6356523.5 -87739228 75026181 81382705
[47,] -6356523.5 -87739228 75026181 81382705
[48,] -6356523.5 -87739228 75026181 81382705
[49,] -6356523.5 -87739228 75026181 81382705
[50,] -6356523.5 -87739228 75026181 81382705

Tamaño de los intervalos de confianza

diff_IC_gdp=predictions$fcst$gdp[,3]-predictions$fcst$gdp[,2]
plot(diff_IC_gdp, 
     main="Longitud de los IC vs cantidad de pronósticos a futuro - GDP", xlab='Cantidad de datos pronosticados en el futuro', ylab='Longitud del IC')

diff_IC_une=predictions$fcst$une[,3]-predictions$fcst$une[,2]
plot(diff_IC_une, main="Longitud de los IC vs cantidad de pronósticos a futuro - UNE", xlab='Cantidad de datos pronosticados en el futuro', ylab='Longitud del IC')