gsm.csv | |

File Size: | 16 kb |

File Type: | csv |

# Cointegration and VECM

The following packages contain a lot of useful functions for Time Series analysis:

```
library(vars)
library(astsa)
library(pracma)
library(urca)
library(tseries)
library(quantmod)
library(tseries)
library(egcm)
library(tsDyn)
```

## Cointegration

Throughout the following, somewhat simplified discussion, we will consider pairs of time series, though more general multivariate time series can be considered. We will also focus on using R to do the modelling and forecasting rather than on the mathematical details.

Suppose we have two series \(Y\) and \(X\). If the two series are stationary, then any relationship between the two can be investiagted using vector autoregression (see previous notes). So, let’s suppose that we have two non-stationary series *but whose first differences are stationary*. Such series are called

**integrated of order \(1\)**.

We can check for stationarity using the “Augmented Dicky-Fuller” test as in the following example. The series \(X\) below is certainly not stationary:

```
set.seed(123)
X <- ts(cumsum(rnorm(500,0,1)))
plot(X, col = "green")
```

`adf.test(X)`

```
##
## Augmented Dickey-Fuller Test
##
## data: X
## Dickey-Fuller = -1.9203, Lag order = 7, p-value = 0.612
## alternative hypothesis: stationary
```

The \(p\) value in the test is not “small enough” to reject the null hypothesis of non-stationarity.

However, the same test on the series of first differences gives:

```
X.diff <- diff(X,1)
adf.test(X.diff)
```

```
##
## Augmented Dickey-Fuller Test
##
## data: X.diff
## Dickey-Fuller = -8.3887, Lag order = 7, p-value = 0.01
## alternative hypothesis: stationary
```

In the test above, the \(p\) value is sufficiently small to reject the null hypothesis of non-stationarity.

So, now we know how to check for stationarity, let’s return to the case in which we have two series, \(X\) and \(Y\), which are both integrated of order 1. Then we can check for “cointegration”.

Before we give a more formal definition, let’s try to understand what this means intuitively. Recall that a stationary series oscillates around its mean; It’s like the series is tied by a rubber band to its mean value and as such can’t arbitrarily wander away from it.

Two series are said to be co-integrated if the difference (*) between them is a stationary series. In this case, it’s like the two series are linked by a rubber band which prevents them from arbitrarily wandering off from each other. Another way of putting this is to say that there is a long run relationship between \(X\) and \(Y\).

Yet another way of understaning cointegrated series is to say that both are driven by the same underlying stochastic process. We can easily simulate this:

```
US <- ts(cumsum(rnorm(500,0,1)))
X <- US + rnorm(500,0,1)
Y <- US + rnorm(500,0,5)
min <- min(min(X), min(Y))
max <- max(max(X), max(Y))
plot(X, col = "green", ylim = c(min,max), ylab="")
lines(Y, col = "red")
```

Before checking for cointegration, we should check that the two series are appropriate for the analysis: they must be integrated of order 1. Thus, we check the first differences of each series for stationarity:

```
X.diff <- diff(X,1)
Y.diff <- diff(Y,1)
adf.test(X.diff)
```

```
##
## Augmented Dickey-Fuller Test
##
## data: X.diff
## Dickey-Fuller = -8.1547, Lag order = 7, p-value = 0.01
## alternative hypothesis: stationary
```

`adf.test(Y.diff)`

```
##
## Augmented Dickey-Fuller Test
##
## data: Y.diff
## Dickey-Fuller = -12.723, Lag order = 7, p-value = 0.01
## alternative hypothesis: stationary
```

The \(p\) value in both cases indicates stationarity in the first differences and hence that both \(X\) and \(Y\) are integrated of order \(1\). We can now check for cointegration: the **difference** between the series should be stationary:

`adf.test(Y-X)`

```
##
## Augmented Dickey-Fuller Test
##
## data: Y - X
## Dickey-Fuller = -9.3888, Lag order = 7, p-value = 0.01
## alternative hypothesis: stationary
```

The low \(p\) value confirms cointegration since the difference between the two series is itself a stationary series.

There is a function in the urca package: ca.jo() which implements a test for cointegration known as Johansen’s Test:

```
YX <- cbind(Y,X)
j.test <- ca.jo(YX)
summary(j.test)
```

```
##
## ######################
## # Johansen-Procedure #
## ######################
##
## Test type: maximal eigenvalue statistic (lambda max) , with linear trend
##
## Eigenvalues (lambda):
## [1] 0.395467597 0.006783299
##
## Values of teststatistic and critical values of test:
##
## test 10pct 5pct 1pct
## r <= 1 | 3.39 6.50 8.18 11.65
## r = 0 | 250.64 12.91 14.90 19.19
##
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
##
## Y.l2 X.l2
## Y.l2 1.000000 1.00000
## X.l2 -1.010875 20.87859
##
## Weights W:
## (This is the loading matrix)
##
## Y.l2 X.l2
## Y.d -1.15290540 -0.0007514497
## X.d 0.05164347 -0.0007792290
```

The table of results looks complicated. The part to focus on for now is:

`j.test@teststat[c(1,2)]`

`## [1] 3.389592 250.643404`

`j.test@cval`

```
## 10pct 5pct 1pct
## r <= 1 | 6.50 8.18 11.65
## r = 0 | 12.91 14.90 19.19
```

The second “test statistic”, \(250.64\) is to be compared with the critical values in the second row of the table above. Since the test statistic is greater than \(12.91\), \(14.90\) and \(19.19\), we can reject the null hyothesis that there is no cointegration at the \(90\%\), \(95\%\) and \(99\%\) confidence levels. Of course, for this simulated data, we already knew that cointegration should hold.

(*) To summarise the definition of cointegration for a pair of series \(Y\) and \(X\), we must have:

- Both \(X\) and \(Y\) are integrated of order \(1\).
- Their difference, \(Y-X\) must be stationary.

This definition is a little too restrictive. In fact a pair of series \(X\) and \(Y\) are cointegrated if:

- Both \(X\) and \(Y\) are integrated of order \(1\).
- There exists \(\beta\) such that \(Y+\beta X\) is stationary.

Here’s another manufactured example to illustrate:

```
set.seed(123)
US <- ts(cumsum(rnorm(500,0,1)))
X <- US + rnorm(500,0,1)
Y <- 5*US + rnorm(500,0,1)
min <- min(min(X), min(Y))
max <- max(max(X), max(Y))
plot(X, col = "green", ylim = c(min,max), ylab="")
lines(Y, col = "red")
```

```
YX <- cbind(Y,X)
j.test <- ca.jo(YX)
summary(j.test)
```

```
##
## ######################
## # Johansen-Procedure #
## ######################
##
## Test type: maximal eigenvalue statistic (lambda max) , with linear trend
##
## Eigenvalues (lambda):
## [1] 0.343255773 0.004691507
##
## Values of teststatistic and critical values of test:
##
## test 10pct 5pct 1pct
## r <= 1 | 2.34 6.50 8.18 11.65
## r = 0 | 209.39 12.91 14.90 19.19
##
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
##
## Y.l2 X.l2
## Y.l2 1.00000 1.0000000
## X.l2 -5.00816 0.2290517
##
## Weights W:
## (This is the loading matrix)
##
## Y.l2 X.l2
## Y.d -0.0002634315 -0.013806678
## X.d 0.2069663467 -0.002478258
```

Comparing the test statistic \(209.39\) with the critical values \(12.91\), \(14.90\) and \(19.19\) means that we have sufficient evidence to accept cointegration. Also notice the cointegration relation:

`j.test@V[,1]`

```
## Y.l2 X.l2
## 1.00000 -5.00816
```

This is known as the “cointegrating vector”. Compare the values here with the coefficients used to construct the simulated series. The result is telling us that a value of has been found so that \(Y + \beta X\) is stationary and the estimate is that \(\beta = -5.00816\). This provides information on the long term relationship between the series; for every unit change in \(X\), we expect a \(5\) unit change in \(Y\).

Let’s construct another example, but this time the effect is felt after 5 lags.

```
set.seed(123)
US <- ts(cumsum(rnorm(505,0,1)))
X <- ts(US[1:500] + rnorm(500,0,1))
Y <- ts(5*US[6:505] + rnorm(500,0,1))
min <- min(min(X), min(Y))
max <- max(max(X), max(Y))
plot(X, col = "green", ylim = c(min,max), ylab="", xlim = c(0,550))
lines(Y, col = "red")
```

```
CP <- cbind(Y,X)
j.test <- ca.jo(CP, K=5)
summary(j.test)
```

```
##
## ######################
## # Johansen-Procedure #
## ######################
##
## Test type: maximal eigenvalue statistic (lambda max) , with linear trend
##
## Eigenvalues (lambda):
## [1] 0.58585030 0.00285626
##
## Values of teststatistic and critical values of test:
##
## test 10pct 5pct 1pct
## r <= 1 | 1.42 6.50 8.18 11.65
## r = 0 | 436.36 12.91 14.90 19.19
##
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
##
## Y.l5 X.l5
## Y.l5 1.000000 1.00000
## X.l5 -5.031411 7.06022
##
## Weights W:
## (This is the loading matrix)
##
## Y.l5 X.l5
## Y.d -0.01137203 -4.739790e-03
## X.d 0.19540961 4.931075e-05
```

Again, we see strong evidence for cointegration which we should since both series share the same underlying stochastic process.

The example above was created so that the series were lagged by \(5\) time steps in their response to their common underlying process \(US\). But how would we know the lag length, \(K\) in \[\text{j.test <- ca.jo(CP, K = 5)}\] when the data is not simulated as such?

There’s a function for that. \(\text{VARselect}\) will estimate the required lag length. Let’s re-run the previous example and pretend we didn’t already know to set the lag length at \(K = 5\):

```
CP <- cbind(Y,X)
VS <- VARselect(CP)
k <- VS$selection[1]
print(k)
```

```
## AIC(n)
## 5
```

```
j.test <- ca.jo(CP, K = k)
summary(j.test)
```

```
##
## ######################
## # Johansen-Procedure #
## ######################
##
## Test type: maximal eigenvalue statistic (lambda max) , with linear trend
##
## Eigenvalues (lambda):
## [1] 0.58585030 0.00285626
##
## Values of teststatistic and critical values of test:
##
## test 10pct 5pct 1pct
## r <= 1 | 1.42 6.50 8.18 11.65
## r = 0 | 436.36 12.91 14.90 19.19
##
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
##
## Y.l5 X.l5
## Y.l5 1.000000 1.00000
## X.l5 -5.031411 7.06022
##
## Weights W:
## (This is the loading matrix)
##
## Y.l5 X.l5
## Y.d -0.01137203 -4.739790e-03
## X.d 0.19540961 4.931075e-05
```

We can now attempt to make a prediction using a VECM model in R. A VECM model is similar to a VAR model in that it attempts to explain the series values using the previous values of both variables (short run effect) *and* using an error correcting term which accounts for the fact that the two series, being cointegrated, have a tendency to pull back toward each other (long run effect).

First, we ask R to build a VECM model as follows:

`vecm <- VECM(CP, lag = (k-1))`

Note that the number of lags entered here is \(1\) less than the number of lags quoted in the VAR model selection (see above $ k = 5$). This is due to the function working on the first differences of the series and hence the number of lags reduces by \(1\).

Next, having built the VECM, we extract the tail of the multivariate series. The number of \(YX\) pairs to extract should equal the number of lags required by the VAR (\(5\) in the present case).

```
YXtail <- tail(CP,k)
prediction <- predict(vecm, newdata = YXtail, n.ahead=20)
print(prediction)
```

```
## Y X
## 501 81.61519 16.82452
## 502 82.20587 16.11749
## 503 81.73098 17.00995
## 504 82.33367 17.85215
## 505 82.43317 16.65940
## 506 82.66853 16.77244
## 507 82.99438 16.86311
## 508 82.98251 16.79673
## 509 83.13479 16.88238
## 510 83.30965 16.91227
## 511 83.44183 16.96221
## 512 83.60618 17.02164
## 513 83.76285 17.02582
## 514 83.91915 17.05784
## 515 84.07955 17.09291
## 516 84.23126 17.12140
## 517 84.38608 17.15461
## 518 84.54208 17.18716
## 519 84.69667 17.21966
## 520 84.85210 17.25268
```

We now add the prediction to the time series plot:

```
plot(X, col = "green", ylim = c(min,max), ylab="", xlim = c(400,550))
lines(Y, col = "red")
Xp <- 501:520
lines(Xp,prediction[,1], col = "darkred")
lines(Xp,prediction[,2], col = "darkgreen")
polygon(c(501,501,550,550),c(min,max,max,min), col = rgb(0.7,0,0,0.1), lty = 0)
```

In the following example, we examine the real gold and silver prices from 1st January 1915 to 1st March 2019 for cointegration, build a VECM and make a prediction for the next 20 months. You can try this example for yourself by downloading and importing the CSV included in this webpage.

```
GSM <- read.csv("GSM.csv")
#extract the Gold and Silver real prices
G <- ts(GSM[,1])
S <- ts(GSM[,2])
min <- min(min(S), min(G))
max <- max(max(S), max(G))
plot(G, col = "green", ylim = c(min,max), ylab="")
lines(S, col = "red")
```

```
#check the series are integrated of order 1
adf.test(diff(G,1))
```

```
##
## Augmented Dickey-Fuller Test
##
## data: diff(G, 1)
## Dickey-Fuller = -8.2573, Lag order = 10, p-value = 0.01
## alternative hypothesis: stationary
```

`adf.test(diff(S,1))`

```
##
## Augmented Dickey-Fuller Test
##
## data: diff(S, 1)
## Dickey-Fuller = -11.296, Lag order = 10, p-value = 0.01
## alternative hypothesis: stationary
```

```
GS <- cbind(G,S)
#select lag length
VS <- VARselect(GS)
k <- VS$selection[1]
#test for cointegration
j.test <- ca.jo(GS, K = k)
summary(j.test)
```

```
##
## ######################
## # Johansen-Procedure #
## ######################
##
## Test type: maximal eigenvalue statistic (lambda max) , with linear trend
##
## Eigenvalues (lambda):
## [1] 0.0171291105 0.0007057782
##
## Values of teststatistic and critical values of test:
##
## test 10pct 5pct 1pct
## r <= 1 | 0.88 6.50 8.18 11.65
## r = 0 | 21.44 12.91 14.90 19.19
##
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
##
## G.l10 S.l10
## G.l10 1.00000 1.00000
## S.l10 -96.58899 -20.37708
##
## Weights W:
## (This is the loading matrix)
##
## G.l10 S.l10
## G.d -9.572174e-05 -0.0035718979
## S.d 4.181583e-04 -0.0001328027
```

```
#build VECM
GSvecm <- VECM(GS, lag = (k-1))
#extract tail of data and predict
GStail <- tail(GS,k)
prediction <- predict(GSvecm, newdata = GStail, n.ahead=20)
#add prediction to plot
x <- length(G)+1
plot(G, col = "green", ylim = c(min,max), ylab="", xlim = c((x-200),(x+50)))
lines(S, col = "red")
Xp <- x:(x+19)
lines(Xp,prediction[,1], col = "darkgreen")
lines(Xp,prediction[,2], col = "darkred")
polygon(c((x+1),(x+1),(x+50),(x+50)),c(min,max,max,min), col = rgb(0.7,0,0,0.1), lty = 0)
```