1. Correlation bw Two Non-stationary TS

When comparing two non-statioanry time series, we can’t use correlation as a guide for deciding if they are linearly dependent or not.

set.seed(4234)

# Say, X1 is ARIMA(1,1,1)
X1  = arima.sim(n = 250, list(order = c(1,1,1), ar = c(0.8), ma = c(-.24) ))

# Say, X2 is another ARIMA(1,1,1), completely independent
X2  = arima.sim(n = 250, list(order = c(1,1,1), ar = c(0.7), ma = c(-.14) ))

# Also say, X3 is linearly dependant on X1 with noise
X3  = 5 + .5*X1 + rnorm(251)


# They have nothing to do with each other, but correlation doesn't says so --
plot(X1, type="l", main=paste("cor=", round(cor(X1,X2), 4) ), ylim=c(-100,50) )
lines(X2, col="red")

cor(X1,X2)
## [1] 0.8559513
cor(X1,X3)
## [1] 0.9979338


1b. Simulate 1000 times to see overall behavior

COR1 = 0
COR2 = 0
for (i in 1:1000) {

    X1  = arima.sim(n = 250, list(order = c(1,1,1), ar = c(0.8), ma = c(-.24) ))
    X2  = arima.sim(n = 250, list(order = c(1,1,1), ar = c(0.8), ma = c(-.24) ))
    X3  = 5 + .5*X1 + rnorm(251)
    COR1[i] = cor(X1,X2)
    COR2[i] = cor(X1,X3)

}
hist(COR1, xlim=c(-1,1),
 main="correlation between two independent non-stationary TS")

sort(COR1)[c(50,950)]
## [1] -0.8080696  0.7940176
hist(COR2, xlim=c(-1,1),
  main="correlation between two linearly depenent non-stationary TS")

sort(COR2)[c(50,950)]
## [1] 0.9821016 0.9987706


2 Engle-Granger Method

E-G Method says to look at the resudual after regressing one non-stationary TS to another.

library(tseries)

# Say, Y1 and Y2 are uncorrelated Random Walks
e1 = rnorm(100, 0, 1)
e2 = rnorm(100, 0, 1)
Y1 = cumsum(e1)
Y2 = cumsum(e2)     # Y1 and Y2 are totally unrelated.


# Say, X1 and X2 are correlated Random Walks
e1 = rnorm(100, 0, 1)
e2 = rnorm(100, 0, 1)
X1 = cumsum(e1)
X2 = .7*X1 + e2 + 3   # X2 is linearly dependent on X1 + (stationary noise)


# Do the regression on Ys
RegY = lm(Y2 ~ Y1)
adf.test(RegY$residuals)  # H0: Non-stationarity can't be rejected -> residuals are non-stationary
## 
##  Augmented Dickey-Fuller Test
## 
## data:  RegY$residuals
## Dickey-Fuller = -3.0132, Lag order = 4, p-value = 0.1569
## alternative hypothesis: stationary
# Do the regression on Xs
RegX = lm(X2 ~ X1)
adf.test(RegX$residuals)  # H0: Non-stationarity is rejected. -> residuals are stationary
## Warning in adf.test(RegX$residuals): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  RegX$residuals
## Dickey-Fuller = -6.1648, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary

2b. Simulate 1000 times for overall behavior

Engle-Granger Method uses the above characteristic to decide over two unrelated non-stationary TS (as Y1 and Y2) vs two linearly dependent non-stationary TS (as X1 and X2).

options(warn = -1) # Turns off all warning messages

# Overall, of course, they are uncorrelated
EG = matrix(rep(0,1000*2), 1000, 2)
for (i in 1:1000) {

    e1 = rnorm(100, 0, 1);   Y1 = cumsum(e1)
    e2 = rnorm(100, 0, 1);   Y2 = cumsum(e2)

    e1 = rnorm(100, 0, 1);   X1 = cumsum(e1)
    e2 = rnorm(100, 0, 1);   X2 = .7*X1 + e2 + 3

    RegY  = lm(Y2 ~ Y1)
    RegX  = lm(X2 ~ X1)
    EG[i,] = c( adf.test(RegY$residuals)$statistic,
                 adf.test(RegX$residuals)$statistic)

}
hist(EG[,1], col = cm.colors(30,   alpha = 0.5), xlim=c(-6,2), ylim=c(0,300) )
par(new=T)

hist(EG[,2], col = heat.colors(30, alpha = 0.5), xlim=c(-6,2), ylim=c(0,300) )

par(new=F)

sum(EG[,1]< -3)
## [1] 177
sum(EG[,2]< -3)
## [1] 993


CODE ONLY

##--------------------------------------
## 1. Correlation bw Two Non-stationary TS

# When comparing two non-statioanry time series, we can't use correlation as
# a guide for deciding if they are linearly dependent or not.

set.seed(4234)

# Say, X1 is ARIMA(1,1,1)
X1  = arima.sim(n = 250, list(order = c(1,1,1), ar = c(0.8), ma = c(-.24) ))

# Say, X2 is another ARIMA(1,1,1), completely independent
X2  = arima.sim(n = 250, list(order = c(1,1,1), ar = c(0.7), ma = c(-.14) ))

# Also say, X3 is linearly dependant on X1 with noise
X3  = 5 + .5*X1 + rnorm(251)


# They have nothing to do with each other, but correlation doesn't says so --
plot(X1, type="l", main=paste("cor=", round(cor(X1,X2), 4) ), ylim=c(-100,50) )
lines(X2, col="red")
cor(X1,X2)
cor(X1,X3)


###----------
### 1b. Simulate 1000 times to see overall behavior

COR1 = 0
COR2 = 0
for (i in 1:1000) {

    X1  = arima.sim(n = 250, list(order = c(1,1,1), ar = c(0.8), ma = c(-.24) ))
    X2  = arima.sim(n = 250, list(order = c(1,1,1), ar = c(0.8), ma = c(-.24) ))
    X3  = 5 + .5*X1 + rnorm(251)
    COR1[i] = cor(X1,X2)
    COR2[i] = cor(X1,X3)

}
hist(COR1, xlim=c(-1,1),
 main="correlation between two independent non-stationary TS")
sort(COR1)[c(50,950)]
hist(COR2, xlim=c(-1,1),
  main="correlation between two linearly depenent non-stationary TS")
sort(COR2)[c(50,950)]


##--------------------------------------
## 2 Engle-Granger Method

#E-G Method says to look at the resudual after regressing one non-stationary TS to another.

library(tseries)

# Say, Y1 and Y2 are uncorrelated Random Walks
e1 = rnorm(100, 0, 1)
e2 = rnorm(100, 0, 1)
Y1 = cumsum(e1)
Y2 = cumsum(e2)     # Y1 and Y2 are totally unrelated.


# Say, X1 and X2 are correlated Random Walks
e1 = rnorm(100, 0, 1)
e2 = rnorm(100, 0, 1)
X1 = cumsum(e1)
X2 = .7*X1 + e2 + 3   # X2 is linearly dependent on X1 + (stationary noise)


# Do the regression on Ys
RegY = lm(Y2 ~ Y1)
adf.test(RegY$residuals)  # H0: Non-stationarity can't be rejected -> residuals are non-stationary


# Do the regression on Xs
RegX = lm(X2 ~ X1)
adf.test(RegX$residuals)  # H0: Non-stationarity is rejected. -> residuals are stationary


###----------
### 2b. Simulate 1000 times for overall behavior

# Engle-Granger Method uses the above characteristic to
# decide over two unrelated non-stationary TS (as Y1 and Y2) vs
# two linearly dependent non-stationary TS (as X1 and X2).


options(warn = -1) # Turns off all warning messages

# Overall, of course, they are uncorrelated
EG = matrix(rep(0,1000*2), 1000, 2)
for (i in 1:1000) {

    e1 = rnorm(100, 0, 1);   Y1 = cumsum(e1)
    e2 = rnorm(100, 0, 1);   Y2 = cumsum(e2)

    e1 = rnorm(100, 0, 1);   X1 = cumsum(e1)
    e2 = rnorm(100, 0, 1);   X2 = .7*X1 + e2 + 3

    RegY  = lm(Y2 ~ Y1)
    RegX  = lm(X2 ~ X1)
    EG[i,] = c( adf.test(RegY$residuals)$statistic,
                 adf.test(RegX$residuals)$statistic)

}
hist(EG[,1], col = cm.colors(30,   alpha = 0.5), xlim=c(-6,2), ylim=c(0,300) )
par(new=T)

hist(EG[,2], col = heat.colors(30, alpha = 0.5), xlim=c(-6,2), ylim=c(0,300) )
par(new=F)

sum(EG[,1]< -3)
sum(EG[,2]< -3)