library(quantmod)
## Loading required package: xts
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## Loading required package: TTR
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(tseries)
library(fImport)
## Warning: package 'fImport' was built under R version 4.3.3
## Loading required package: timeDate
## Loading required package: timeSeries
## Warning: package 'timeSeries' was built under R version 4.3.3
## 
## Attaching package: 'timeSeries'
## The following object is masked from 'package:zoo':
## 
##     time<-
## The following objects are masked from 'package:graphics':
## 
##     lines, points
library(httr)
library(RCurl)
library(rjson)
library(parsedate)
## Warning: package 'parsedate' was built under R version 4.3.3
library(urca)
#Datos
#---------------
getSymbols("BTC-USD", quote="Close", from="2016-01-01", periodicty = "daily")
## [1] "BTC-USD"
plot(`BTC-USD`[,4], type="l")

Sin tendencia y sin constante

precio_BTC <- na.omit(`BTC-USD`[,4])
# Test ADF sin tendencia ni constante
#---------------------------------------
y_nada=ur.df(precio_BTC, type = "none", lags =0)
summary(y_nada)
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression none 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7594.7  -147.0     1.6   161.2  7265.1 
## 
## Coefficients:
##          Estimate Std. Error t value Pr(>|t|)
## z.lag.1 0.0007178  0.0006534   1.099    0.272
## 
## Residual standard error: 889.6 on 3006 degrees of freedom
## Multiple R-squared:  0.0004013,  Adjusted R-squared:  6.878e-05 
## F-statistic: 1.207 on 1 and 3006 DF,  p-value: 0.2721
## 
## 
## Value of test-statistic is: 1.0986 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62

1.1027 > -1.95 – Acepto H0 –> Raíz unitaria (NO ESTACIONARIA)

Con constante “drift”

# Test ADF sin tendencia ni constante
#-------------------------------------------
y_drift=ur.df(precio_BTC, type="drift", lags=0)
summary(y_drift)
## 
## ############################################### 
## # 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 
## -7580.5  -156.6   -19.2   150.0  7268.0 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.210e+01  2.340e+01   0.945    0.345
## z.lag.1     7.633e-05  9.424e-04   0.081    0.935
## 
## Residual standard error: 889.6 on 3005 degrees of freedom
## Multiple R-squared:  2.183e-06,  Adjusted R-squared:  -0.0003306 
## F-statistic: 0.006561 on 1 and 3005 DF,  p-value: 0.9354
## 
## 
## Value of test-statistic is: 0.081 1.0495 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.43 -2.86 -2.57
## phi1  6.43  4.59  3.78

0.0855 > -2.86 –> Acepto H0 –> Raíz unitaria (NO ESTACIONARIA)

Con tendencia y constante

# Test ADF con tendencia ni constante
#-------------------------------------------
y_trend=ur.df(precio_BTC, type="trend", lags=0)
summary(y_trend)
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression trend 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7525.6  -171.2     6.4   142.3  7291.5 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -21.410210  33.164244  -0.646   0.5186  
## z.lag.1      -0.001958   0.001448  -1.353   0.1763  
## tt            0.053134   0.028713   1.851   0.0643 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 889.2 on 3004 degrees of freedom
## Multiple R-squared:  0.001141,   Adjusted R-squared:  0.0004758 
## F-statistic: 1.716 on 2 and 3004 DF,  p-value: 0.18
## 
## 
## Value of test-statistic is: -1.3526 1.8417 1.7155 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau3 -3.96 -3.41 -3.12
## phi2  6.09  4.68  4.03
## phi3  8.27  6.25  5.34

-1.3486 > -3.41 –> Acepto H0 –> Raíz unitaria (NO ESTACIONARIA)

Sin constante y tendencia: No estacionaria Con constante y sin tendencia: No estacionaria Con constante y tendencia: No estacionaria

En al menos una especificación resultó que la serie del bitcoin no es estacionaria.

Cointegración

H0: No hay relación de cointegración H1: Al menos existe una relación de cointegración

library(Quandl)
# Datos para cointegración
#----------------------------------------------
Poil  <- Quandl("FRED/MCOILWTICO", type ="ts", collapse = "monthly", start_date = "1995-01-01", end_date = "2017-12-01")
Pgas  <- Quandl("FRED/GASREGCOVM", type ="ts", collapse = "monthly", start_date = "1995-01-01", end_date = "2017-12-01")
#Calculamos los logaritmos
#----------------------------
LPoil <- log(Poil)
LPgas <- log(Pgas)
library(vars)
## Warning: package 'vars' was built under R version 4.3.3
## Loading required package: MASS
## Loading required package: strucchange
## Warning: package 'strucchange' was built under R version 4.3.3
## Loading required package: sandwich
## Loading required package: lmtest
## 
## Attaching package: 'lmtest'
## The following object is masked from 'package:RCurl':
## 
##     reset
#Fusión de las dos variables en un DF
#---------------------------------------
v1 <- cbind(LPoil, LPgas)
#Var (modelo de vectores autorregresivos)
#----------------------------------------------
var <- VARselect(v1, type = "trend")
var
## $selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##     10      4      3     10 
## 
## $criteria
##                    1             2             3             4             5
## AIC(n) -1.100552e+01 -1.131266e+01 -1.138368e+01 -1.140805e+01 -1.140308e+01
## HQ(n)  -1.097305e+01 -1.125854e+01 -1.130791e+01 -1.131063e+01 -1.128401e+01
## SC(n)  -1.092469e+01 -1.117794e+01 -1.119508e+01 -1.116555e+01 -1.110670e+01
## FPE(n)  1.660974e-05  1.221741e-05  1.137992e-05  1.110633e-05  1.116214e-05
##                    6             7             8             9            10
## AIC(n) -1.141632e+01 -1.141124e+01 -1.144593e+01 -1.145337e+01 -1.147193e+01
## HQ(n)  -1.127560e+01 -1.124887e+01 -1.126192e+01 -1.124771e+01 -1.124462e+01
## SC(n)  -1.106605e+01 -1.100708e+01 -1.098789e+01 -1.094144e+01 -1.090611e+01
## FPE(n)  1.101598e-05  1.107305e-05  1.069659e-05  1.061883e-05  1.042536e-05
nlags <- var$selection["SC(n)"]
nlags
## SC(n) 
##     3
#Causalidad de Granger
#--------------------------------
v1.CA <- ca.jo(v1, ecdet = "trend", type = "trace", K=nlags, spec = "transitory", season = 12)
summary(v1.CA)
## 
## ###################### 
## # Johansen-Procedure # 
## ###################### 
## 
## Test type: trace statistic , with linear trend in cointegration 
## 
## Eigenvalues (lambda):
## [1]  1.037680e-01  1.873353e-02 -1.699366e-19
## 
## Values of teststatistic and critical values of test:
## 
##           test 10pct  5pct  1pct
## r <= 1 |  5.16 10.49 12.25 16.26
## r = 0  | 35.07 22.76 25.32 30.45
## 
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
## 
##               LPoil.l1    LPgas.l1    trend.l1
## LPoil.l1  1.0000000000  1.00000000   1.0000000
## LPgas.l1 -1.6267117480  3.75026554 -10.7065612
## trend.l1  0.0009487176 -0.01919425   0.1692144
## 
## Weights W:
## (This is the loading matrix)
## 
##             LPoil.l1     LPgas.l1      trend.l1
## LPoil.d -0.006283109 -0.008367640 -2.545462e-19
## LPgas.d  0.138252677 -0.003301629  8.623684e-20

Values of teststatistic and critical values of test:

      test 10pct  5pct  1pct

r <= 1 | 5.16 10.49 12.25 16.26

r = 0 | 35.07 22.76 25.32 30.45

test = 35.07 vc(5%) = 25.32 –> Rechazo H0

Al rechazarse H0, existe al menos una relación de cointegracón.

test = 5.16 < vc(5%) = 12.25 –> Acepto H0 r<= 1 relaciones de cointegración

r>0 & r <= –> r=1 existe una relación de cointegración.

v1.CA <- ca.jo(v1, ecdet = "trend", type = "eigen", K=nlags, spec = "transitory", season = 12)
summary(v1.CA)
## 
## ###################### 
## # Johansen-Procedure # 
## ###################### 
## 
## Test type: maximal eigenvalue statistic (lambda max) , with linear trend in cointegration 
## 
## Eigenvalues (lambda):
## [1]  1.037680e-01  1.873353e-02 -1.699366e-19
## 
## Values of teststatistic and critical values of test:
## 
##           test 10pct  5pct  1pct
## r <= 1 |  5.16 10.49 12.25 16.26
## r = 0  | 29.91 16.85 18.96 23.65
## 
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
## 
##               LPoil.l1    LPgas.l1    trend.l1
## LPoil.l1  1.0000000000  1.00000000   1.0000000
## LPgas.l1 -1.6267117480  3.75026554 -10.7065612
## trend.l1  0.0009487176 -0.01919425   0.1692144
## 
## Weights W:
## (This is the loading matrix)
## 
##             LPoil.l1     LPgas.l1      trend.l1
## LPoil.d -0.006283109 -0.008367640 -2.545462e-19
## LPgas.d  0.138252677 -0.003301629  8.623684e-20
     test 10pct  5pct  1pct

r <= 1 | 5.16 10.49 12.25 16.26 r = 0 | 29.91 16.85 18.96 23.65

test = 29.91 > 18.96 –> rechazo H0 –> existe al menos una relación de cointegración

test 5.16 < 12.25 –> Acepto H0 –> Existe una relación única de cointegración.

Existe una relación de cointegración

Eigenvectors, normalised to first column: (These are the cointegration relations)

          LPoil.l1    LPgas.l1    trend.l1

LPoil.l1 1.0000000000 1.00000000 1.0000000 LPgas.l1 -1.6267117480 3.75026554 -10.7065612 trend.l1 0.0009487176 -0.01919425 0.1692144

Weights W: (This is the loading matrix)

        LPoil.l1     LPgas.l1                      trend.l1

LPoil.d -0.006283109 -0.008367640 -0.00000000000000000025454624 LPgas.d 0.138252677 -0.003301629 0.00000000000000000008623684

Causalidad en sentido de Granger

tsData <- EuStockMarkets[,1:2]
head(tsData)
##          DAX    SMI
## [1,] 1628.75 1678.1
## [2,] 1613.63 1688.5
## [3,] 1606.51 1678.6
## [4,] 1621.04 1684.1
## [5,] 1618.16 1686.6
## [6,] 1610.61 1671.6
library(lmtest)
#Causalidad en sentido Granger
#------------------------------------
resultado_granger1 <- grangertest(DAX ~SMI, order = 3, data = tsData)
summary(resultado_granger1)
##      Res.Df           Df           F             Pr(>F)        
##  Min.   :1850   Min.   :-3   Min.   :8.497   Min.   :1.32e-05  
##  1st Qu.:1851   1st Qu.:-3   1st Qu.:8.497   1st Qu.:1.32e-05  
##  Median :1852   Median :-3   Median :8.497   Median :1.32e-05  
##  Mean   :1852   Mean   :-3   Mean   :8.497   Mean   :1.32e-05  
##  3rd Qu.:1852   3rd Qu.:-3   3rd Qu.:8.497   3rd Qu.:1.32e-05  
##  Max.   :1853   Max.   :-3   Max.   :8.497   Max.   :1.32e-05  
##                 NA's   :1    NA's   :1       NA's   :1

F = 8.497 y si p valor es 0.0000132 < 0.05. Rechazo H0. Existe causalidad en sentido de Granger.

SMI causa o explica el DAX

#Causalidad en sentido Granger
#------------------------------------
resultado_granger2 <- grangertest(SMI ~ DAX, order = 3, data = tsData)
summary(resultado_granger2)
##      Res.Df           Df           F             Pr(>F)       
##  Min.   :1850   Min.   :-3   Min.   :2.658   Min.   :0.04689  
##  1st Qu.:1851   1st Qu.:-3   1st Qu.:2.658   1st Qu.:0.04689  
##  Median :1852   Median :-3   Median :2.658   Median :0.04689  
##  Mean   :1852   Mean   :-3   Mean   :2.658   Mean   :0.04689  
##  3rd Qu.:1852   3rd Qu.:-3   3rd Qu.:2.658   3rd Qu.:0.04689  
##  Max.   :1853   Max.   :-3   Max.   :2.658   Max.   :0.04689  
##                 NA's   :1    NA's   :1       NA's   :1

F = 2.658 p valor = 0.04689 < 0.05 – > rechazo H0. Existe relación de causalidad del DAX hacia SMI

Dax causa SMI
SMI causa DAX

Relación de causalidad bidireccional

Ejmplo simulado de causalidad en sentido de Granger

set.seed(12345)
datos <- data.frame(
  tiempo = 1:100,
  serie_x = rnorm(100),
  serie_y = 0.5 * lag(rnorm(100)) + rnorm(100))
head(datos)
##   tiempo    serie_x    serie_y
## 1      1  0.5855288 -1.3241830
## 2      2  0.7094660 -1.2073713
## 3      3 -0.1093033  0.4547310
## 4      4 -0.4534972  0.3959846
## 5      5  0.6058875  0.9018910
## 6      6 -1.8179560 -0.1628122
# Prueba de Granger
#-------------------------------------------------
resultado_granger1 <- grangertest(serie_x ~serie_y, order = 3, data = datos)
summary(resultado_granger1)
##      Res.Df            Df           F             Pr(>F)      
##  Min.   :90.00   Min.   :-3   Min.   :1.553   Min.   :0.2064  
##  1st Qu.:90.75   1st Qu.:-3   1st Qu.:1.553   1st Qu.:0.2064  
##  Median :91.50   Median :-3   Median :1.553   Median :0.2064  
##  Mean   :91.50   Mean   :-3   Mean   :1.553   Mean   :0.2064  
##  3rd Qu.:92.25   3rd Qu.:-3   3rd Qu.:1.553   3rd Qu.:0.2064  
##  Max.   :93.00   Max.   :-3   Max.   :1.553   Max.   :0.2064  
##                  NA's   :1    NA's   :1       NA's   :1

F = 1.553 P-Valor ó 0.2064 > 0.05. Aceptamos H0 –> Y no causa a X

# Prueba de Granger
#-------------------------------------------------
resultado_granger2 <- grangertest(serie_y ~serie_x, order = 3, data = datos)
summary(resultado_granger2)
##      Res.Df            Df           F             Pr(>F)      
##  Min.   :90.00   Min.   :-3   Min.   :1.487   Min.   :0.2235  
##  1st Qu.:90.75   1st Qu.:-3   1st Qu.:1.487   1st Qu.:0.2235  
##  Median :91.50   Median :-3   Median :1.487   Median :0.2235  
##  Mean   :91.50   Mean   :-3   Mean   :1.487   Mean   :0.2235  
##  3rd Qu.:92.25   3rd Qu.:-3   3rd Qu.:1.487   3rd Qu.:0.2235  
##  Max.   :93.00   Max.   :-3   Max.   :1.487   Max.   :0.2235  
##                  NA's   :1    NA's   :1       NA's   :1

F = 1.487 P-Valor ó 0.2235 > 0.05. Aceptamos H0 –> X no causa a Y

X no causa a Y
Y no causa a X

No hay relación en sentido de Granger entre X e Y.