Let’s say we have some fluctuating ‘scariness’ time series that represents how frightening a horror movie is every second. This is let’s say 80 observations sampled from a random normal distribution (mean=5, sd=2)
set.seed(123)
TS_HowScary <- rnorm(80, 15, 3)
head(TS_HowScary)
## [1] 13.31857 14.30947 19.67612 15.21153 15.38786 20.14519
Now let’s Lag the ‘scariness’ time series by some fixed interval (3 units). We will call this the ‘Heart Rate’ time series. In this simulation. Heart Rate is a perfect lagged function of Scariness. You know that any peak in ‘scariness’ will evoke the same peak in ‘heart rate’ three units later becuase that’s how we modeled the data. When you lag a TS, you will create missing values for the first N of the lag function. You need to extrapolate to fill in values so that both time series are of equal length
# creates a lagged version of teh first time series (shifted by 3 turns)
TS_HeartRate <- lag(TS_HowScary, 3)
# count missing obsercations sum(is.na(TS_HeartRate))
# rule=2 allows for extarpolation across leading or trailing missing
# observations
TS_HeartRate <- zoo::na.approx(TS_HeartRate, rule = 2)
# count missing obsercations sum(is.na(TS_HeartRate))
Plot them. In this case, HR is a perfect lagged function of Scariness level with the first three observations imputed/extrapolated.
TS_HowScary <- as.ts(TS_HowScary)
TS_HeartRate <- as.ts(TS_HeartRate)
# plot the first TS then add the second using lines
plot(TS_HowScary, col = "gray", ylim = c(0, 30))
lines(TS_HeartRate, col = "darkred", lty = 2, lwd = 2)
The syntax of the Granger test is
lmtest::grangertest(X,Y, order=1) where X is the first time
series and y is the second time series. Order is the number of lags from
the first time series to test at, Default is lag=1.
oes Heart rate cause scariness? Most certainly not! So the test of heart rate grannger causing scariness should be ns. Let’s run it and find out.
library(lmtest)
# Does the scary time series 'granger cause' the HR time series?
res1 <- lmtest::grangertest(TS_HeartRate, TS_HowScary, order = 3)
print(res1)
## Granger causality test
##
## Model 1: TS_HowScary ~ Lags(TS_HowScary, 1:3) + Lags(TS_HeartRate, 1:3)
## Model 2: TS_HowScary ~ Lags(TS_HowScary, 1:3)
## Res.Df Df F Pr(>F)
## 1 70
## 2 73 -3 0.7251 0.5404
In this case, we know that it should because we simulated the data so that HR is a perfect lagged function of Scariness
res2 <- lmtest::grangertest(TS_HowScary, TS_HeartRate, order = 3)
print(res2)
## Granger causality test
##
## Model 1: TS_HeartRate ~ Lags(TS_HeartRate, 1:3) + Lags(TS_HowScary, 1:3)
## Model 2: TS_HeartRate ~ Lags(TS_HeartRate, 1:3)
## Res.Df Df F Pr(>F)
## 1 70
## 2 73 -3 9.5459e+33 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Voila! There you have it… Scariness Granger Causes HR in this simuklated example at p<.00000000000001