Se presenta una visión de los comandos vec en R a través de un ejemplo extenso. Tenemos datos mensuales sobre los precios de venta promedio de viviendas en cuatro ciudades de Texas: Austin, Dallas, Houston y San Antonio. En el conjunto de datos, estos precios medios de viviendas están contenidos en las variables de austin, dallas, houston y sa. Las series comienzan en enero de 1990 y van hasta diciembre de 2003, para un total de 168 observaciones. Para el presente ejemplo trabajaremos con las series Dallas y Houston.

library(foreign)
txhprice=read.dta("c:/datos/txhprice.dta")
attach(txhprice)
txh<-data.frame(dallas,houston)
txhs<-ts(txh,frequency = 12,start = 1990)
plot.ts(txhs)

Grafico de las series

El siguiente grafico representa nuestros datos.

plot.ts(txhs)

A excepción de la caída a finales de 1991, los precios de la vivienda en Dallas y Houston parecen estrechamente relacionadas. Aunque los precios medios en las dos ciudades serán diferentes debido a las variaciones de recursos y otros factores, si los mercados de vivienda se convierten en demasiado diferentes, la gente y las empresas migran, con lo que los precios medios de la vivienda vuelven hacia atrás entre sí. Por lo tanto, se espera que la serie de precios medios de vivienda en Houston esté cointegrada con la serie de los precios medios de vivienda en Dallas.

Paso 1: Confirmación de raiz unitaria

El gráfico indican que las series tienen tendencia y en potencia son procesos I(1). En un mercado competitivo, los precios actuales y pasados contienen toda la información disponible, así que el precio de mañana será un paseo aleatorio de los precios de hoy. Se opta por usar la prueba Dickey - Fuller para investigar la presencia de una raíz unitaria en cada serie.

Prueba para Dallas
library(urca)
df.test=ur.df(dallas, type=c("drift"),lags=0)
summary(df.test)
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.107791 -0.023694  0.001706  0.024089  0.133754 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  0.24088    0.20509   1.175    0.242
## z.lag.1     -0.02006    0.01731  -1.159    0.248
## 
## Residual standard error: 0.04161 on 165 degrees of freedom
## Multiple R-squared:  0.008074,   Adjusted R-squared:  0.002062 
## F-statistic: 1.343 on 1 and 165 DF,  p-value: 0.2482
## 
## 
## Value of test-statistic is: -1.1589 1.1768 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1  6.52  4.63  3.81
df.test1=ur.df(diff(dallas), type=c("drift"),lags=0)
summary(df.test1)
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.11625 -0.02423  0.00173  0.02262  0.11509 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.004313   0.003085   1.398    0.164    
## z.lag.1     -1.324594   0.073867 -17.932   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03963 on 164 degrees of freedom
## Multiple R-squared:  0.6622, Adjusted R-squared:  0.6602 
## F-statistic: 321.6 on 1 and 164 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -17.9321 160.7813 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1  6.52  4.63  3.81
Prueba para Houston
library(urca)
df.test=ur.df(houston, type=c("drift"),lags=0)
summary(df.test)
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.288315 -0.024380  0.000381  0.032488  0.281205 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  0.36274    0.23723   1.529    0.128
## z.lag.1     -0.03062    0.02026  -1.511    0.133
## 
## Residual standard error: 0.05627 on 165 degrees of freedom
## Multiple R-squared:  0.01365,    Adjusted R-squared:  0.007674 
## F-statistic: 2.284 on 1 and 165 DF,  p-value: 0.1326
## 
## 
## Value of test-statistic is: -1.5112 1.6289 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1  6.52  4.63  3.81
df.test1=ur.df(diff(houston), type=c("drift"),lags=0)
summary(df.test1)
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.259018 -0.026154 -0.000476  0.028616  0.163980 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.006239   0.003800   1.642    0.103    
## z.lag.1     -1.510858   0.067161 -22.496   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.04883 on 164 degrees of freedom
## Multiple R-squared:  0.7553, Adjusted R-squared:  0.7538 
## F-statistic: 506.1 on 1 and 164 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -22.4961 253.0371 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1  6.52  4.63  3.81

La prueba realizada muestra que en las series no son estacionarias a nivel peros si lo son en primeras diferencias.

La prueba de cointegración que utilizaremos incluye el caso en el que todas las variables son estacionarias, por lo que podemos aplazar la prueba formal de raíz unitaria hasta que probemos la cointegración. Las tendencias temporales en los datos parecen ser aproximadamente lineal, por lo que especificaremos esto al modelar estas series, el cual es el valor por defecto en una VEC.

Paso 2: Identificacion del Numero de Retardos a Incluir en el VEC

Para la prueba de cointegración o el ajuste de los VECMs, hay que especificar el número de rezagos a incluir. Basándose en el trabajo de Tsay (1984) y Paulsen (1984), Nielsen (2001) ha demostrado que los métodos implementados en varsoc pueden ser utilizados para determinar el orden de retardos para un modelo VAR con variables I(1). El orden del VEC correspondiente es siempre uno menos que el VAR. El comnado vec hace este ajuste automáticamente, por lo que siempre se referirá al orden del VAR subyacente. La salida siguiente usa VARorder para determinar el orden de rezago del VAR de los precios medios de las viviendas en Dallas y Houston.

library("MTS")
VARorder(txhs,maxp=10)
## selected order: aic =  9 
## selected order: bic =  2 
## selected order: hq =  2 
## Summary table:  
##        p      AIC      BIC       HQ     M(p) p-value
##  [1,]  0  -9.4086  -9.4086  -9.4086   0.0000  0.0000
##  [2,]  1 -12.7386 -12.6642 -12.7084 521.8347  0.0000
##  [3,]  2 -12.8645 -12.7158 -12.8042  26.4670  0.0000
##  [4,]  3 -12.8650 -12.6419 -12.7744   7.2375  0.1239
##  [5,]  4 -12.8399 -12.5424 -12.7192   3.3510  0.5009
##  [6,]  5 -12.8238 -12.4519 -12.6729   4.6126  0.3294
##  [7,]  6 -12.8575 -12.4113 -12.6764  11.7542  0.0193
##  [8,]  7 -12.8792 -12.3585 -12.6679   9.8694  0.0427
##  [9,]  8 -12.9063 -12.3113 -12.6648  10.5009  0.0328
## [10,]  9 -12.9596 -12.2902 -12.6879  13.9763  0.0074
## [11,] 10 -12.9330 -12.1892 -12.6311   2.8648  0.5807

Para cada longitud de rezago, la tabla reporta, el Criterio de informacion de Akaike (AIC), el criterio informacion bayesiana de Schwarz (BIC), Criterio de informacion de Hannan y Quinn (HQ) y el estadístico Chi Cuadrado (M(p)). Los criterios se basan en la teoria de la informacion y se supone que indican la informacion relativa perdida cuando los datos se ajustan usando diferentes especificaciones. La longitud de rezago que produce el valor minimo del estadistico de la informacion es la especificacion elegida. En este ejemplo, los criterios BIC y HQ eligen 2 rezagos. Puede ser el caso que este conjunto de pruebas estadisticas no dejen una direccion clara. Segun Lutkepohl (2005), el BIC y HQ proporcionan estimaciones consistentes del verdadero orden del retardo, mientras que AIC sobreestima el orden del retardo con una alta probabilidad.

1.3. Paso 3: Identificación del Número de Relaciones Cointegrantes

Las pruebas de cointegración implementadas en el paquete Vars se basan en el método de Johansen. Si el log likelihood del modelo no restringido que incluye las ecuaciones de cointegración es significativamente diferente del log likelihood del modelo restringido que no incluye las ecuaciones de cointegración, se rechaza la hipótesis nula de no cointegración. Aquí usamos ca.jo para determinar el número de ecuaciones de cointegración. ca.jo lleva a cabo el procedimiento de Johansen en un conjunto determinado de datos. Los estadísticos “trace” o “eigen” son reportados y la matriz de vectores propios, así como la matriz de carga.

library(vars)
## Loading required package: MASS
## Loading required package: strucchange
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## Loading required package: sandwich
## Loading required package: lmtest
## 
## Attaching package: 'vars'
## The following object is masked from 'package:MTS':
## 
##     VAR
johansen<-ca.jo(txhs, type="eigen",spec="transitory",ecdet="none",K=2)
summary(johansen)
## 
## ###################### 
## # Johansen-Procedure # 
## ###################### 
## 
## Test type: maximal eigenvalue statistic (lambda max) , with linear trend 
## 
## Eigenvalues (lambda):
## [1] 0.244975014 0.001074729
## 
## Values of teststatistic and critical values of test:
## 
##           test 10pct  5pct  1pct
## r <= 1 |  0.18  6.50  8.18 11.65
## r = 0  | 46.65 12.91 14.90 19.19
## 
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
## 
##             dallas.l1 houston.l1
## dallas.l1   1.0000000   1.000000
## houston.l1 -0.8675936   1.038759
## 
## Weights W:
## (This is the loading matrix)
## 
##            dallas.l1   houston.l1
## dallas.d  -0.3038799 -0.002749322
## houston.d  0.5027143 -0.002774421

En la salida anterior, rechazamos firmemente la hipótesis nula de no cointegración y no rechazamos la hipótesis nula de a lo sumo una ecuación de cointegración. Por lo tanto concluimos que hay una ecuación de cointegración en el modelo bivariado.

1.4. Paso 4: Estimacion del Modelo VEC

Habiendo determinado que existe una ecuación de cointegración entre las series Dallas y Houston, ahora queremos estimar los parámetros de una cointegración bivariada VEC para estas dos series utilizando cajorls.

vecm1<-cajorls(johansen,r=1)
vecm1
## $rlm
## 
## Call:
## lm(formula = substitute(form1), data = data.mat)
## 
## Coefficients:
##              dallas.d  houston.d
## ect1         -0.30388   0.50271 
## constant      0.51883  -0.84564 
## dallas.dl1   -0.16473  -0.06197 
## houston.dl1  -0.09984  -0.33284 
## 
## 
## $beta
##                  ect1
## dallas.l1   1.0000000
## houston.l1 -0.8675936
summary(vecm1$rlm)
## Response dallas.d :
## 
## Call:
## lm(formula = dallas.d ~ ect1 + constant + dallas.dl1 + houston.dl1 - 
##     1, data = data.mat)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.107171 -0.024624  0.002201  0.023697  0.119788 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## ect1        -0.30388    0.09085  -3.345 0.001023 ** 
## constant     0.51883    0.15388   3.372 0.000934 ***
## dallas.dl1  -0.16473    0.08794  -1.873 0.062827 .  
## houston.dl1 -0.09984    0.06508  -1.534 0.126987    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03855 on 162 degrees of freedom
## Multiple R-squared:  0.1692, Adjusted R-squared:  0.1487 
## F-statistic: 8.247 on 4 and 162 DF,  p-value: 4.438e-06
## 
## 
## Response houston.d :
## 
## Call:
## lm(formula = houston.d ~ ect1 + constant + dallas.dl1 + houston.dl1 - 
##     1, data = data.mat)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.251415 -0.024548  0.001749  0.026632  0.137640 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## ect1         0.50271    0.10688   4.703 5.45e-06 ***
## constant    -0.84564    0.18104  -4.671 6.26e-06 ***
## dallas.dl1  -0.06197    0.10345  -0.599     0.55    
## houston.dl1 -0.33284    0.07657  -4.347 2.43e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.04535 on 162 degrees of freedom
## Multiple R-squared:  0.3737, Adjusted R-squared:  0.3582 
## F-statistic: 24.17 on 4 and 162 DF,  p-value: 1.082e-15