Granger Causality Test


##載入套件

library(lmtest)
library(dplyr)
library(ggplot2)
library(openxlsx)
library(vars)
library(tseries)
library(jtools)
library(forecast)
library(urca)


##建立時間序列

season <- read.xlsx("E:/summer program/new data/byseason1.xlsx", sheet = '內插')

cover.timeseries <- ts(season$aboveground距平, start = c(1997,1),frequency = 4)
leaf.timeseries <- ts(season$leaf距平, start = c(1997,1),frequency = 4)
biomass.timeseries <- ts(season$aboveground距平, start = c(1997,1),frequency = 4)
relative.timeseries <- ts(season$relative距平, start = c(1997,1),frequency = 4)
oni.timeseries <- ts(season$ONI0, start = c(1997,1),frequency = 4)
pmm.timeseries <- ts(season$PMM0, start = c(1997,1),frequency = 4)
nino4.timeseries <- ts(season$NINO40, start = c(1997,1),frequency = 4)
cp.timeseries <- ts(season$CPENSO0, start = c(1997,1),frequency = 4)
temp.timeseries <- ts(season$TX距平0, start = c(1997,1),frequency = 4)
rain.timeseries <- ts(season$PP距平0, start = c(1997,1),frequency = 4)
salt.timeseries <- ts(season$Salt距平0, start = c(1997,1),frequency = 4)
tempde.timeseries <- ts(season$TXde0, start = c(1997,1),frequency = 4)


##檢測平穩性

ADF法

adf.test(na.omit(cover.timeseries))
adf.test(na.omit(leaf.timeseries))
adf.test(na.omit(relative.timeseries))
adf.test(na.omit(biomass.timeseries))

adf.test(na.omit(pmm.timeseries))
adf.test(na.omit(oni.timeseries))
adf.test(na.omit(nino4.timeseries))
adf.test(na.omit(cp.timeseries))

adf.test(na.omit(temp.timeseries))
adf.test(na.omit(rain.timeseries))
adf.test(na.omit(salt.timeseries))
adf.test(na.omit(tempde.timeseries))
## 
##  Augmented Dickey-Fuller Test
## 
## data:  na.omit(cover.timeseries)
## Dickey-Fuller = -1.2938, Lag order = 4, p-value = 0.8655
## alternative hypothesis: stationary
## 
## 
##  Augmented Dickey-Fuller Test
## 
## data:  na.omit(leaf.timeseries)
## Dickey-Fuller = -2.7016, Lag order = 4, p-value = 0.2881
## alternative hypothesis: stationary
## 
## 
##  Augmented Dickey-Fuller Test
## 
## data:  na.omit(relative.timeseries)
## Dickey-Fuller = -1.9812, Lag order = 4, p-value = 0.5836
## alternative hypothesis: stationary
## 
## 
##  Augmented Dickey-Fuller Test
## 
## data:  na.omit(biomass.timeseries)
## Dickey-Fuller = -1.2938, Lag order = 4, p-value = 0.8655
## alternative hypothesis: stationary
## 
## 
##  Augmented Dickey-Fuller Test
## 
## data:  na.omit(pmm.timeseries)
## Dickey-Fuller = -2.5903, Lag order = 4, p-value = 0.3325
## alternative hypothesis: stationary
## 
## 
##  Augmented Dickey-Fuller Test
## 
## data:  na.omit(oni.timeseries)
## Dickey-Fuller = -3.6285, Lag order = 4, p-value = 0.03462
## alternative hypothesis: stationary
## 
## 
##  Augmented Dickey-Fuller Test
## 
## data:  na.omit(nino4.timeseries)
## Dickey-Fuller = -3.4569, Lag order = 4, p-value = 0.04983
## alternative hypothesis: stationary
## 
## 
##  Augmented Dickey-Fuller Test
## 
## data:  na.omit(cp.timeseries)
## Dickey-Fuller = -3.834, Lag order = 4, p-value = 0.02091
## alternative hypothesis: stationary
## 
## 
##  Augmented Dickey-Fuller Test
## 
## data:  na.omit(temp.timeseries)
## Dickey-Fuller = -2.2262, Lag order = 4, p-value = 0.483
## alternative hypothesis: stationary
## 
## 
##  Augmented Dickey-Fuller Test
## 
## data:  na.omit(rain.timeseries)
## Dickey-Fuller = -4.7815, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary
## 
## 
##  Augmented Dickey-Fuller Test
## 
## data:  na.omit(salt.timeseries)
## Dickey-Fuller = -2.4348, Lag order = 4, p-value = 0.3975
## alternative hypothesis: stationary
## 
## 
##  Augmented Dickey-Fuller Test
## 
## data:  na.omit(tempde.timeseries)
## Dickey-Fuller = -2.2284, Lag order = 4, p-value = 0.4821
## alternative hypothesis: stationary


KPSS法

kpss.test(na.omit(cover.timeseries))
kpss.test(na.omit(leaf.timeseries))
kpss.test(na.omit(relative.timeseries))
kpss.test(na.omit(biomass.timeseries))

kpss.test(na.omit(pmm.timeseries))
kpss.test(na.omit(oni.timeseries))
kpss.test(na.omit(nino4.timeseries))
kpss.test(na.omit(cp.timeseries))

kpss.test(na.omit(temp.timeseries))
kpss.test(na.omit(rain.timeseries))
kpss.test(na.omit(salt.timeseries))
kpss.test(na.omit(tempde.timeseries))
## 
##  KPSS Test for Level Stationarity
## 
## data:  na.omit(cover.timeseries)
## KPSS Level = 0.55658, Truncation lag parameter = 3, p-value = 0.02892
## 
## 
##  KPSS Test for Level Stationarity
## 
## data:  na.omit(leaf.timeseries)
## KPSS Level = 0.15951, Truncation lag parameter = 3, p-value = 0.1
## 
## 
##  KPSS Test for Level Stationarity
## 
## data:  na.omit(relative.timeseries)
## KPSS Level = 1.0795, Truncation lag parameter = 3, p-value = 0.01
## 
## 
##  KPSS Test for Level Stationarity
## 
## data:  na.omit(biomass.timeseries)
## KPSS Level = 0.55658, Truncation lag parameter = 3, p-value = 0.02892
## 
## 
##  KPSS Test for Level Stationarity
## 
## data:  na.omit(pmm.timeseries)
## KPSS Level = 0.8593, Truncation lag parameter = 3, p-value = 0.01
## 
## 
##  KPSS Test for Level Stationarity
## 
## data:  na.omit(oni.timeseries)
## KPSS Level = 0.11392, Truncation lag parameter = 3, p-value = 0.1
## 
## 
##  KPSS Test for Level Stationarity
## 
## data:  na.omit(nino4.timeseries)
## KPSS Level = 0.19939, Truncation lag parameter = 3, p-value = 0.1
## 
## 
##  KPSS Test for Level Stationarity
## 
## data:  na.omit(cp.timeseries)
## KPSS Level = 0.22521, Truncation lag parameter = 3, p-value = 0.1
## 
## 
##  KPSS Test for Level Stationarity
## 
## data:  na.omit(temp.timeseries)
## KPSS Level = 0.78001, Truncation lag parameter = 3, p-value = 0.01
## 
## 
##  KPSS Test for Level Stationarity
## 
## data:  na.omit(rain.timeseries)
## KPSS Level = 0.12981, Truncation lag parameter = 3, p-value = 0.1
## 
## 
##  KPSS Test for Level Stationarity
## 
## data:  na.omit(salt.timeseries)
## KPSS Level = 0.17933, Truncation lag parameter = 3, p-value = 0.1
## 
## 
##  KPSS Test for Level Stationarity
## 
## data:  na.omit(tempde.timeseries)
## KPSS Level = 0.16384, Truncation lag parameter = 3, p-value = 0.1


##檢測差分

ADF法

ndiffs(cover.timeseries, test="adf")
ndiffs(leaf.timeseries, test="adf")
ndiffs(relative.timeseries, test="adf")
ndiffs(biomass.timeseries, test="adf")
ndiffs(pmm.timeseries, test="adf")
ndiffs(oni.timeseries, test="adf")
ndiffs(nino4.timeseries, test="adf")
ndiffs(cp.timeseries, test="adf")
ndiffs(temp.timeseries, test="adf")
ndiffs(rain.timeseries, test="adf")
ndiffs(salt.timeseries, test="adf")
ndiffs(tempde.timeseries, test="adf")
## [1] 1
## [1] 0
## [1] 0
## [1] 1
## [1] 0
## [1] 0
## [1] 0
## [1] 0
## [1] 0
## [1] 0
## [1] 0
## [1] 0


KPSS法

ndiffs(cover.timeseries, test="kpss")
ndiffs(leaf.timeseries, test="kpss")
ndiffs(relative.timeseries, test="kpss")
ndiffs(biomass.timeseries, test="kpss")
ndiffs(pmm.timeseries, test="kpss")
ndiffs(oni.timeseries, test="kpss")
ndiffs(nino4.timeseries, test="kpss")
ndiffs(cp.timeseries, test="kpss")
ndiffs(temp.timeseries, test="kpss")
ndiffs(rain.timeseries, test="kpss")
ndiffs(salt.timeseries, test="kpss")
ndiffs(tempde.timeseries, test="kpss")
## [1] 1
## [1] 0
## [1] 1
## [1] 1
## [1] 1
## [1] 0
## [1] 0
## [1] 0
## [1] 1
## [1] 0
## [1] 0
## [1] 0

cover, biomass必做差分;relative, PMM, temp可做可不做

##檢測季節性差分

nsdiffs(cover.timeseries)
nsdiffs(leaf.timeseries)
nsdiffs(relative.timeseries)
nsdiffs(biomass.timeseries)
nsdiffs(pmm.timeseries)
nsdiffs(oni.timeseries)
nsdiffs(nino4.timeseries)
nsdiffs(cp.timeseries)
nsdiffs(temp.timeseries)
nsdiffs(rain.timeseries)
nsdiffs(salt.timeseries)
nsdiffs(tempde.timeseries)
## [1] 0
## [1] 0
## [1] 0
## [1] 0
## [1] 0
## [1] 0
## [1] 0
## [1] 0
## [1] 0
## [1] 0
## [1] 0
## [1] 0



##差分比較

coverdiff <- diff(cover.timeseries)

par(mfrow=c(2,2))
plot.ts(cover.timeseries,main="原始cover")
plot.ts(coverdiff,main="一階差分cover")
acf(cover.timeseries, na = na.omit, ,main="原始cover"); acf(coverdiff, na = na.omit, main="一階差分cover")

coverdiff <- diff(cover.timeseries)
par(mfrow=c(2,2))
plot.ts(cover.timeseries,main="原始cover")
plot.ts(coverdiff,main="一階差分cover")
acf(cover.timeseries, na = na.omit, ,main="原始cover");acf(coverdiff, na = na.omit, main="一階差分cover")

relativediff <- diff(relative.timeseries)

#par(mfrow=c(2,2))
#plot.ts(relative.timeseries,main="原始relative")
#plot.ts(relativediff,main="一階差分relative")
#acf(relative.timeseries, na = na.omit, ,main="原始relative"); acf(relativediff, na = na.omit, main="一階差分relative")
pmmdiff <- diff(pmm.timeseries)
#par(mfrow=c(2,2))
#plot.ts(pmm.timeseries,main="原始pmm")
#plot.ts(pmmdiff,main="一階差分pmm")
#acf(pmm.timeseries, na = na.omit, ,main="原始pmm");acf(pmmdiff, na = na.omit, main="一階差分pmm")
tempdiff <- diff(temp.timeseries)

#par(mfrow=c(2,2))
#plot.ts(temp.timeseries,main="原始temp")
#plot.ts(tempdiff,main="一階差分temp")
#(temp.timeseries, na = na.omit, ,main="原始temp");acf(tempdiff, na = na.omit, main="一階差分temp")



##建立校準時間序列

因為差分的原理為第一項與第二項之差為新序列之第一項,因此計算a與b之格蘭傑因果時,其結果會與現實差了一個時間單位。然而當此差距會包含0時,將導致一盲區無法被計算(沒有自回歸),因此需要建立一個向後推移一單位之時間序列。

cover1.timeseries <- ts(season$aboveground距平1, start = c(1997,1),frequency = 4)
cover1.timeseries <- ts(season$cover距平1, start = c(1997,1),frequency = 4)
relative1.timeseries <- ts(season$relative距平1, start = c(1997,1),frequency = 4)
pmm1.timeseries <- ts(season$PMM1, start = c(1997,1),frequency = 4)
temp1.timeseries <- ts(season$TX距平1, start = c(1997,1),frequency = 4)



##建立校準差分序列

cover1diff <- diff(cover1.timeseries)
relative1diff <- diff(relative1.timeseries)
cover1diff <- diff(cover1.timeseries)
pmm1diff <- diff(pmm1.timeseries)
temp1diff <- diff(temp1.timeseries)



##cover與PMM的滯後選擇 & Granger Causality Test

cover差分與PMM

bio_pmm <- cbind(pmm.timeseries, coverdiff)
VARselect(y=na.omit(bio_pmm), lag.max = 8, type = c("const"))
## $selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      5      1      1      5 
## 
## $criteria
##                 1          2          3          4          5          6
## AIC(n) -1.0142979 -1.0510566 -1.0587417 -1.0573471 -1.0670205 -0.9913257
## HQ(n)  -0.9402701 -0.9276770 -0.8860102 -0.8352637 -0.7955852 -0.6705385
## SC(n)  -0.8288989 -0.7420582 -0.6261439 -0.5011500 -0.3872240 -0.1879298
## FPE(n)  0.3626879  0.3497067  0.3472703  0.3481844  0.3455014  0.3737179
##                  7          8
## AIC(n) -0.90715411 -0.8607084
## HQ(n)  -0.53701509 -0.4412175
## SC(n)   0.01984113  0.1898862
## FPE(n)  0.40810781  0.4296932

AIC best at 5, BIC(SC) best at 1.

grangertest(coverdiff ~ pmm.timeseries, order = 1, data = bio_pmm)
grangertest(pmm.timeseries ~ coverdiff, order = 1, data = bio_pmm)
grangertest(coverdiff ~ pmm.timeseries, order = 5, data = bio_pmm)
grangertest(pmm.timeseries ~ coverdiff, order = 5, data = bio_pmm)
## Granger causality test
## 
## Model 1: coverdiff ~ Lags(coverdiff, 1:1) + Lags(pmm.timeseries, 1:1)
## Model 2: coverdiff ~ Lags(coverdiff, 1:1)
##   Res.Df Df     F Pr(>F)
## 1     79                
## 2     80 -1 0.164 0.6866
## Granger causality test
## 
## Model 1: pmm.timeseries ~ Lags(pmm.timeseries, 1:1) + Lags(coverdiff, 1:1)
## Model 2: pmm.timeseries ~ Lags(pmm.timeseries, 1:1)
##   Res.Df Df      F Pr(>F)  
## 1     79                   
## 2     80 -1 3.8065 0.0546 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Granger causality test
## 
## Model 1: coverdiff ~ Lags(coverdiff, 1:5) + Lags(pmm.timeseries, 1:5)
## Model 2: coverdiff ~ Lags(coverdiff, 1:5)
##   Res.Df Df      F Pr(>F)
## 1     67                 
## 2     72 -5 0.7377 0.5978
## Granger causality test
## 
## Model 1: pmm.timeseries ~ Lags(pmm.timeseries, 1:5) + Lags(coverdiff, 1:5)
## Model 2: pmm.timeseries ~ Lags(pmm.timeseries, 1:5)
##   Res.Df Df      F Pr(>F)  
## 1     67                   
## 2     72 -5 2.1196 0.0737 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

cover lead PMM by 1 and 5. (實質上cover與PMM同時/領先4季)

cover1差分與PMM

bio1_pmm <- cbind(pmm.timeseries, cover1diff)
VARselect(y=na.omit(bio1_pmm), lag.max = 8, type = c("const"))
## $selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      4      1      1      4 
## 
## $criteria
##                1         2         3         4         5         6         7
## AIC(n)  4.329266  4.342998  4.250210  4.248838  4.311245  4.304390  4.375736
## HQ(n)   4.403293  4.466378  4.422941  4.470922  4.582680  4.625177  4.745875
## SC(n)   4.514665  4.651996  4.682807  4.805036  4.991042  5.107786  5.302731
## FPE(n) 75.895018 76.968340 70.196567 70.186961 74.851557 74.549391 80.371998
##                8
## AIC(n)  4.405685
## HQ(n)   4.825176
## SC(n)   5.456280
## FPE(n) 83.238458

AIC best at 4, BIC(SC) best at 1.

grangertest(cover1diff ~ pmm.timeseries, order = 1, data = bio1_pmm)
grangertest(pmm.timeseries ~ cover1diff, order = 1, data = bio1_pmm)
grangertest(cover1diff ~ pmm.timeseries, order = 4, data = bio1_pmm)
grangertest(pmm.timeseries ~ cover1diff, order = 4, data = bio1_pmm)
## Granger causality test
## 
## Model 1: cover1diff ~ Lags(cover1diff, 1:1) + Lags(pmm.timeseries, 1:1)
## Model 2: cover1diff ~ Lags(cover1diff, 1:1)
##   Res.Df Df      F Pr(>F)
## 1     79                 
## 2     80 -1 0.4728 0.4937
## Granger causality test
## 
## Model 1: pmm.timeseries ~ Lags(pmm.timeseries, 1:1) + Lags(cover1diff, 1:1)
## Model 2: pmm.timeseries ~ Lags(pmm.timeseries, 1:1)
##   Res.Df Df      F  Pr(>F)  
## 1     79                    
## 2     80 -1 2.8061 0.09786 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Granger causality test
## 
## Model 1: cover1diff ~ Lags(cover1diff, 1:4) + Lags(pmm.timeseries, 1:4)
## Model 2: cover1diff ~ Lags(cover1diff, 1:4)
##   Res.Df Df      F Pr(>F)
## 1     70                 
## 2     74 -4 0.5734 0.6828
## Granger causality test
## 
## Model 1: pmm.timeseries ~ Lags(pmm.timeseries, 1:4) + Lags(cover1diff, 1:4)
## Model 2: pmm.timeseries ~ Lags(pmm.timeseries, 1:4)
##   Res.Df Df      F Pr(>F)
## 1     70                 
## 2     74 -4 1.9923 0.1051

cover lead PMM by 1. (實質上cover與PMM同時)

##cover差分與PMM差分

bio_pmmdiff <- cbind(pmmdiff, coverdiff)
VARselect(y=na.omit(bio_pmmdiff), lag.max = 8, type = c("const"))
## $selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      2      2      2      2 
## 
## $criteria
##                 1          2          3          4          5          6
## AIC(n) -0.8650512 -1.0477954 -0.9991713 -0.9772263 -0.9519731 -0.8752591
## HQ(n)  -0.7910234 -0.9244157 -0.8264398 -0.7551429 -0.6805378 -0.5544719
## SC(n)  -0.6796522 -0.7387969 -0.5665735 -0.4210292 -0.2721765 -0.0718632
## FPE(n)  0.4210659  0.3508490  0.3685859  0.3772292  0.3876272  0.4197116
##                 7          8
## AIC(n) -0.7878636 -0.8624107
## HQ(n)  -0.4177246 -0.4429198
## SC(n)   0.1391317  0.1881839
## FPE(n)  0.4598139  0.4289623

AIC, BIC(SC) best at 2.

grangertest(coverdiff ~ pmmdiff, order = 2, data = bio_pmmdiff)
grangertest(pmmdiff ~ coverdiff, order = 2, data = bio_pmmdiff)
## Granger causality test
## 
## Model 1: coverdiff ~ Lags(coverdiff, 1:2) + Lags(pmmdiff, 1:2)
## Model 2: coverdiff ~ Lags(coverdiff, 1:2)
##   Res.Df Df     F Pr(>F)
## 1     76                
## 2     78 -2 0.056 0.9455
## Granger causality test
## 
## Model 1: pmmdiff ~ Lags(pmmdiff, 1:2) + Lags(coverdiff, 1:2)
## Model 2: pmmdiff ~ Lags(pmmdiff, 1:2)
##   Res.Df Df      F Pr(>F)
## 1     76                 
## 2     78 -2 0.9167 0.4042

No significant result.

##cover與CPENSO的滯後選擇 & Granger

cover差分與CPENSO

bio_cp <- cbind(cp.timeseries, coverdiff)
VARselect(y=na.omit(bio_cp), lag.max = 8, type = c("const"))
## $selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      3      2      2      3 
## 
## $criteria
##                  1           2          3           4          5           6
## AIC(n) -3.32683904 -3.47991151 -3.5065544 -3.47740155 -3.4540930 -3.42831736
## HQ(n)  -3.25131008 -3.35402991 -3.3303201 -3.25081466 -3.1771534 -3.10102520
## SC(n)  -3.13711686 -3.16370788 -3.0638693 -2.90823502 -2.7584450 -2.60618793
## FPE(n)  0.03590989  0.03082394  0.0300371  0.03096887  0.0317688  0.03270223
##                  7           8
## AIC(n) -3.33540562 -3.26958404
## HQ(n)  -2.95776082 -2.84158660
## SC(n)  -2.38679473 -2.19449171
## FPE(n)  0.03604363  0.03871911

AIC best at 3, BIC(SC) best at 1.

grangertest(coverdiff ~ cp.timeseries, order = 1, data = bio_cp)
grangertest(cp.timeseries ~ coverdiff, order = 1, data = bio_cp)
grangertest(coverdiff ~ cp.timeseries, order = 3, data = bio_cp)
grangertest(cp.timeseries ~ coverdiff, order = 3, data = bio_cp)
## Granger causality test
## 
## Model 1: coverdiff ~ Lags(coverdiff, 1:1) + Lags(cp.timeseries, 1:1)
## Model 2: coverdiff ~ Lags(coverdiff, 1:1)
##   Res.Df Df      F Pr(>F)
## 1     76                 
## 2     77 -1 0.1825 0.6704
## Granger causality test
## 
## Model 1: cp.timeseries ~ Lags(cp.timeseries, 1:1) + Lags(coverdiff, 1:1)
## Model 2: cp.timeseries ~ Lags(cp.timeseries, 1:1)
##   Res.Df Df     F Pr(>F)
## 1     76                
## 2     77 -1 0.128 0.7215
## Granger causality test
## 
## Model 1: coverdiff ~ Lags(coverdiff, 1:3) + Lags(cp.timeseries, 1:3)
## Model 2: coverdiff ~ Lags(coverdiff, 1:3)
##   Res.Df Df      F Pr(>F)
## 1     70                 
## 2     73 -3 1.1799 0.3237
## Granger causality test
## 
## Model 1: cp.timeseries ~ Lags(cp.timeseries, 1:3) + Lags(coverdiff, 1:3)
## Model 2: cp.timeseries ~ Lags(cp.timeseries, 1:3)
##   Res.Df Df      F Pr(>F)
## 1     70                 
## 2     73 -3 1.5869 0.2003

No significant result.

cover1差分與CPENSO

bio1_cp <- cbind(cp.timeseries, cover1diff)
VARselect(y=na.omit(bio1_cp), lag.max = 8, type = c("const"))
## $selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      3      2      2      3 
## 
## $criteria
##               1        2        3        4        5        6        7        8
## AIC(n) 2.000314 1.836317 1.800780 1.808594 1.875617 1.951404 2.057411 2.084002
## HQ(n)  2.075338 1.961356 1.975835 2.033665 2.150704 2.276506 2.432529 2.509136
## SC(n)  2.188571 2.150078 2.240046 2.373365 2.565893 2.767184 2.998696 3.150791
## FPE(n) 7.392060 6.276084 6.061529 6.117268 6.555093 7.092803 7.919126 8.177701

AIC best at 3, BIC(SC) best at 2.

grangertest(cover1diff ~ cp.timeseries, order = 2, data = bio1_cp)
grangertest(cp.timeseries ~ cover1diff, order = 2, data = bio1_cp)
grangertest(cover1diff ~ cp.timeseries, order = 3, data = bio1_cp)
grangertest(cp.timeseries ~ cover1diff, order = 3, data = bio1_cp)
## Granger causality test
## 
## Model 1: cover1diff ~ Lags(cover1diff, 1:2) + Lags(cp.timeseries, 1:2)
## Model 2: cover1diff ~ Lags(cover1diff, 1:2)
##   Res.Df Df      F Pr(>F)
## 1     74                 
## 2     76 -2 1.3046 0.2774
## Granger causality test
## 
## Model 1: cp.timeseries ~ Lags(cp.timeseries, 1:2) + Lags(cover1diff, 1:2)
## Model 2: cp.timeseries ~ Lags(cp.timeseries, 1:2)
##   Res.Df Df      F Pr(>F)
## 1     74                 
## 2     76 -2 0.8809 0.4187
## Granger causality test
## 
## Model 1: cover1diff ~ Lags(cover1diff, 1:3) + Lags(cp.timeseries, 1:3)
## Model 2: cover1diff ~ Lags(cover1diff, 1:3)
##   Res.Df Df      F Pr(>F)
## 1     71                 
## 2     74 -3 1.1311 0.3424
## Granger causality test
## 
## Model 1: cp.timeseries ~ Lags(cp.timeseries, 1:3) + Lags(cover1diff, 1:3)
## Model 2: cp.timeseries ~ Lags(cp.timeseries, 1:3)
##   Res.Df Df      F Pr(>F)
## 1     71                 
## 2     74 -3 0.4525 0.7163

No significant result.

##cover與Salt的滯後選擇 & Granger Causality Test

cover差分與Salt

bio_salt <- cbind(salt.timeseries, coverdiff)
VARselect(y=na.omit(bio_salt), lag.max = 8, type = c("const"))
## $selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      2      2      1      2 
## 
## $criteria
##                 1          2          3          4          5          6
## AIC(n) -2.0013017 -2.1037585 -2.0286503 -2.0782186 -2.0553771 -2.0522108
## HQ(n)  -1.9272739 -1.9803788 -1.8559188 -1.8561352 -1.7839418 -1.7314236
## SC(n)  -1.8159026 -1.7947601 -1.5960526 -1.5220214 -1.3755806 -1.2488149
## FPE(n)  0.1351708  0.1220454  0.1316563  0.1254442  0.1285914  0.1293622
##                 7          8
## AIC(n) -2.0603694 -2.0663946
## HQ(n)  -1.6902304 -1.6469037
## SC(n)  -1.1333742 -1.0158000
## FPE(n)  0.1288071  0.1286873

AIC, BIC(SC) best at 2.

grangertest(coverdiff ~ salt.timeseries, order = 1, data = bio_salt)
grangertest(salt.timeseries ~ coverdiff, order = 1, data = bio_salt)
grangertest(coverdiff ~ salt.timeseries, order = 2, data = bio_salt)
grangertest(salt.timeseries ~ coverdiff, order = 2, data = bio_salt)
## Granger causality test
## 
## Model 1: coverdiff ~ Lags(coverdiff, 1:1) + Lags(salt.timeseries, 1:1)
## Model 2: coverdiff ~ Lags(coverdiff, 1:1)
##   Res.Df Df      F Pr(>F)
## 1     79                 
## 2     80 -1 0.9089 0.3433
## Granger causality test
## 
## Model 1: salt.timeseries ~ Lags(salt.timeseries, 1:1) + Lags(coverdiff, 1:1)
## Model 2: salt.timeseries ~ Lags(salt.timeseries, 1:1)
##   Res.Df Df      F Pr(>F)
## 1     79                 
## 2     80 -1 0.0019 0.9653
## Granger causality test
## 
## Model 1: coverdiff ~ Lags(coverdiff, 1:2) + Lags(salt.timeseries, 1:2)
## Model 2: coverdiff ~ Lags(coverdiff, 1:2)
##   Res.Df Df      F Pr(>F)
## 1     76                 
## 2     78 -2 0.8818 0.4182
## Granger causality test
## 
## Model 1: salt.timeseries ~ Lags(salt.timeseries, 1:2) + Lags(coverdiff, 1:2)
## Model 2: salt.timeseries ~ Lags(salt.timeseries, 1:2)
##   Res.Df Df     F Pr(>F)
## 1     76                
## 2     78 -2 0.931 0.3986

No significant result.

cover1差分與Salt

bio1_salt <- cbind(salt.timeseries, cover1diff)
VARselect(y=na.omit(bio1_salt), lag.max = 8, type = c("const"))
## $selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      3      1      1      3 
## 
## $criteria
##                1         2         3         4         5         6         7
## AIC(n)  3.379094  3.332888  3.290836  3.342437  3.411806  3.446901  3.485658
## HQ(n)   3.453122  3.456268  3.463568  3.564521  3.683241  3.767688  3.855797
## SC(n)   3.564493  3.641887  3.723434  3.898635  4.091602  4.250297  4.412653
## FPE(n) 29.346687 28.030253 26.894611 28.353812 30.449437 31.625786 33.002668
##                8
## AIC(n)  3.495164
## HQ(n)   3.914655
## SC(n)   4.545759
## FPE(n) 33.488041

AIC best at 1, BIC(SC) best at 3.

grangertest(cover1diff ~ salt.timeseries, order = 1, data = bio1_salt)
grangertest(salt.timeseries ~ cover1diff, order = 1, data = bio1_salt)
grangertest(cover1diff ~ salt.timeseries, order = 3, data = bio1_salt)
grangertest(salt.timeseries ~ cover1diff, order = 3, data = bio1_salt)
## Granger causality test
## 
## Model 1: cover1diff ~ Lags(cover1diff, 1:1) + Lags(salt.timeseries, 1:1)
## Model 2: cover1diff ~ Lags(cover1diff, 1:1)
##   Res.Df Df      F Pr(>F)
## 1     79                 
## 2     80 -1 0.5277 0.4697
## Granger causality test
## 
## Model 1: salt.timeseries ~ Lags(salt.timeseries, 1:1) + Lags(cover1diff, 1:1)
## Model 2: salt.timeseries ~ Lags(salt.timeseries, 1:1)
##   Res.Df Df      F Pr(>F)
## 1     79                 
## 2     80 -1 0.5716 0.4519
## Granger causality test
## 
## Model 1: cover1diff ~ Lags(cover1diff, 1:3) + Lags(salt.timeseries, 1:3)
## Model 2: cover1diff ~ Lags(cover1diff, 1:3)
##   Res.Df Df      F Pr(>F)
## 1     73                 
## 2     76 -3 1.5766 0.2024
## Granger causality test
## 
## Model 1: salt.timeseries ~ Lags(salt.timeseries, 1:3) + Lags(cover1diff, 1:3)
## Model 2: salt.timeseries ~ Lags(salt.timeseries, 1:3)
##   Res.Df Df      F Pr(>F)
## 1     73                 
## 2     76 -3 1.2911 0.2839

No significant result.

##leaf與CPENSO的滯後選擇 & Granger

leaf與CPENSO

leaf_cp <- cbind(leaf.timeseries, cp.timeseries)
VARselect(y=na.omit(leaf_cp), lag.max = 8, type = c("const"))
## $selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      2      2      1      2 
## 
## $criteria
##                  1           2           3           4           5           6
## AIC(n) -3.32560421 -3.40520545 -3.38987453 -3.36077394 -3.29852737 -3.32029344
## HQ(n)  -3.25058057 -3.28016605 -3.21481936 -3.13570302 -3.02344068 -2.99519099
## SC(n)  -3.13734727 -3.09144389 -2.95060833 -2.79600312 -2.60825192 -2.50451336
## FPE(n)  0.03595412  0.03321427  0.03375278  0.03479602  0.03710876  0.03642075
##                  7           8
## AIC(n) -3.22247168 -3.18527212
## HQ(n)  -2.84735347 -2.76013815
## SC(n)  -2.28118698 -2.11848279
## FPE(n)  0.04033233  0.04209346

AIC best at 3, BIC(SC) best at 1.

grangertest(leaf.timeseries ~ cp.timeseries, order = 1, data = leaf_cp)
grangertest(cp.timeseries ~ leaf.timeseries, order = 1, data = leaf_cp)
grangertest(leaf.timeseries ~ cp.timeseries, order = 2, data = leaf_cp)
grangertest(cp.timeseries ~ leaf.timeseries, order = 2, data = leaf_cp)
## Granger causality test
## 
## Model 1: leaf.timeseries ~ Lags(leaf.timeseries, 1:1) + Lags(cp.timeseries, 1:1)
## Model 2: leaf.timeseries ~ Lags(leaf.timeseries, 1:1)
##   Res.Df Df      F Pr(>F)
## 1     77                 
## 2     78 -1 0.8517  0.359
## Granger causality test
## 
## Model 1: cp.timeseries ~ Lags(cp.timeseries, 1:1) + Lags(leaf.timeseries, 1:1)
## Model 2: cp.timeseries ~ Lags(cp.timeseries, 1:1)
##   Res.Df Df      F Pr(>F)
## 1     77                 
## 2     78 -1 0.9212 0.3402
## Granger causality test
## 
## Model 1: leaf.timeseries ~ Lags(leaf.timeseries, 1:2) + Lags(cp.timeseries, 1:2)
## Model 2: leaf.timeseries ~ Lags(leaf.timeseries, 1:2)
##   Res.Df Df      F Pr(>F)
## 1     74                 
## 2     76 -2 1.0416  0.358
## Granger causality test
## 
## Model 1: cp.timeseries ~ Lags(cp.timeseries, 1:2) + Lags(leaf.timeseries, 1:2)
## Model 2: cp.timeseries ~ Lags(cp.timeseries, 1:2)
##   Res.Df Df      F Pr(>F)
## 1     74                 
## 2     76 -2 0.3126 0.7325

No significant result.

##leaf與Temp的滯後選擇 & Granger

leaf與Tempde

tempde_leaf <- cbind(tempde.timeseries, leaf.timeseries)
VARselect(y=na.omit(tempde_leaf), lag.max = 8, type = c("const"))
## $selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      1      1      1      1 
## 
## $criteria
##                  1           2           3           4          5           6
## AIC(n) -3.72270958 -3.66976608 -3.65499839 -3.63299038 -3.6795957 -3.69345965
## HQ(n)  -3.64917218 -3.54720374 -3.48341111 -3.41237817 -3.4099586 -3.37479757
## SC(n)  -3.53870432 -3.36309064 -3.22565277 -3.08097459 -3.0049097 -2.89610351
## FPE(n)  0.02417038  0.02549213  0.02588863  0.02649611  0.0253367  0.02505538
##                  7           8
## AIC(n) -3.60122193 -3.54621964
## HQ(n)  -3.23353491 -3.12950768
## SC(n)  -2.68119561 -2.50352314
## FPE(n)  0.02757835  0.02928053

AIC best at 3, BIC best at 1.

grangertest(leaf.timeseries ~ tempde.timeseries, order = 1, data = tempde_leaf)
grangertest(tempde.timeseries ~ leaf.timeseries, order = 1, data = tempde_leaf)
grangertest(leaf.timeseries ~ tempde.timeseries, order = 3, data = tempde_leaf)
grangertest(tempde.timeseries ~ leaf.timeseries, order = 3, data = tempde_leaf)
## Granger causality test
## 
## Model 1: leaf.timeseries ~ Lags(leaf.timeseries, 1:1) + Lags(tempde.timeseries, 1:1)
## Model 2: leaf.timeseries ~ Lags(leaf.timeseries, 1:1)
##   Res.Df Df     F Pr(>F)
## 1     80                
## 2     81 -1 0.005 0.9438
## Granger causality test
## 
## Model 1: tempde.timeseries ~ Lags(tempde.timeseries, 1:1) + Lags(leaf.timeseries, 1:1)
## Model 2: tempde.timeseries ~ Lags(tempde.timeseries, 1:1)
##   Res.Df Df      F  Pr(>F)  
## 1     80                    
## 2     81 -1 3.7659 0.05583 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Granger causality test
## 
## Model 1: leaf.timeseries ~ Lags(leaf.timeseries, 1:3) + Lags(tempde.timeseries, 1:3)
## Model 2: leaf.timeseries ~ Lags(leaf.timeseries, 1:3)
##   Res.Df Df      F Pr(>F)
## 1     74                 
## 2     77 -3 0.6092 0.6111
## Granger causality test
## 
## Model 1: tempde.timeseries ~ Lags(tempde.timeseries, 1:3) + Lags(leaf.timeseries, 1:3)
## Model 2: tempde.timeseries ~ Lags(tempde.timeseries, 1:3)
##   Res.Df Df      F Pr(>F)
## 1     74                 
## 2     77 -3 1.1966  0.317

leaf lead temp by 1.

##leaf與Salt的滯後選擇 & Granger Causality Test

leaf與salt

salt_leaf <- cbind(salt.timeseries, leaf.timeseries)
VARselect(y=na.omit(salt_leaf), lag.max = 8, type = c("const"))
## $selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      3      1      1      3 
## 
## $criteria
##                 1          2         3          4          5          6
## AIC(n) -1.9692776 -1.9230704 -1.989260 -1.9476322 -1.8709809 -1.8400644
## HQ(n)  -1.8957402 -1.8005080 -1.817673 -1.7270200 -1.6013438 -1.5214023
## SC(n)  -1.7852723 -1.6163949 -1.559914 -1.3956164 -1.1962950 -1.0427082
## FPE(n)  0.1395691  0.1462132  0.136940  0.1429302  0.1546043  0.1598896
##                 7          8
## AIC(n) -1.7674580 -1.7436445
## HQ(n)  -1.3997710 -1.3269326
## SC(n)  -0.8474317 -0.7009480
## FPE(n)  0.1725686  0.1775936

AIC best at 5, BIC best at 1.

grangertest(leaf.timeseries ~ salt.timeseries, order = 1, data = salt_leaf)
grangertest(salt.timeseries ~ leaf.timeseries, order = 1, data = salt_leaf)
grangertest(leaf.timeseries ~ salt.timeseries, order = 5, data = salt_leaf)
grangertest(salt.timeseries ~ leaf.timeseries, order = 5, data = salt_leaf)
## Granger causality test
## 
## Model 1: leaf.timeseries ~ Lags(leaf.timeseries, 1:1) + Lags(salt.timeseries, 1:1)
## Model 2: leaf.timeseries ~ Lags(leaf.timeseries, 1:1)
##   Res.Df Df      F Pr(>F)
## 1     80                 
## 2     81 -1 0.7125 0.4011
## Granger causality test
## 
## Model 1: salt.timeseries ~ Lags(salt.timeseries, 1:1) + Lags(leaf.timeseries, 1:1)
## Model 2: salt.timeseries ~ Lags(salt.timeseries, 1:1)
##   Res.Df Df      F Pr(>F)
## 1     80                 
## 2     81 -1 2.3372 0.1303
## Granger causality test
## 
## Model 1: leaf.timeseries ~ Lags(leaf.timeseries, 1:5) + Lags(salt.timeseries, 1:5)
## Model 2: leaf.timeseries ~ Lags(leaf.timeseries, 1:5)
##   Res.Df Df      F Pr(>F)
## 1     68                 
## 2     73 -5 0.7927 0.5587
## Granger causality test
## 
## Model 1: salt.timeseries ~ Lags(salt.timeseries, 1:5) + Lags(leaf.timeseries, 1:5)
## Model 2: salt.timeseries ~ Lags(salt.timeseries, 1:5)
##   Res.Df Df      F Pr(>F)
## 1     68                 
## 2     73 -5 1.8067 0.1232

No significant result.

##Temp與PMM的滯後選擇 & Granger

Temp與PMM

temp_pmm <- cbind(temp.timeseries, pmm.timeseries)
VARselect(y=na.omit(temp_pmm), lag.max = 8, type = c("const"))
## $selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      3      3      3      3 
## 
## $criteria
##                1         2           3          4          5         6
## AIC(n) 0.2589601 0.1756239 -0.01794214 0.03566198 0.06342978 0.1163690
## HQ(n)  0.3270095 0.2890395  0.14083963 0.23980997 0.31294398 0.4112494
## SC(n)  0.4278695 0.4571395  0.37617963 0.54238997 0.68276398 0.8483094
## FPE(n) 1.2956506 1.1922818  0.98287977 1.03779402 1.06828936 1.1283142
##                7         8
## AIC(n) 0.1906305 0.1807359
## HQ(n)  0.5308772 0.5663488
## SC(n)  1.0351772 1.1378888
## FPE(n) 1.2181721 1.2099386

AIC, BIC(SC) best at 3.

grangertest(temp.timeseries ~ pmm.timeseries, order = 3, data = temp_pmm)
grangertest(pmm.timeseries ~ temp.timeseries, order = 3, data = temp_pmm)
## Granger causality test
## 
## Model 1: temp.timeseries ~ Lags(temp.timeseries, 1:3) + Lags(pmm.timeseries, 1:3)
## Model 2: temp.timeseries ~ Lags(temp.timeseries, 1:3)
##   Res.Df Df     F  Pr(>F)  
## 1     86                   
## 2     89 -3 3.098 0.03099 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Granger causality test
## 
## Model 1: pmm.timeseries ~ Lags(pmm.timeseries, 1:3) + Lags(temp.timeseries, 1:3)
## Model 2: pmm.timeseries ~ Lags(pmm.timeseries, 1:3)
##   Res.Df Df      F  Pr(>F)  
## 1     86                    
## 2     89 -3 3.9452 0.01093 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

PMM 領先 Temp 3季; Temp 領先 PMM 3季.

Temp差分與PMM

tempdiff_pmm <- cbind(tempdiff, pmm.timeseries)
VARselect(y=na.omit(tempdiff_pmm), lag.max = 8, type = c("const"))
## $selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      4      3      3      4 
## 
## $criteria
##                1         2         3          4         5         6         7
## AIC(n) 0.5346566 0.3465812 0.1300685 0.08766693 0.1421680 0.2154293 0.2888470
## HQ(n)  0.6031355 0.4607128 0.2898527 0.29310386 0.3932576 0.5121715 0.6312419
## SC(n)  0.7047192 0.6300189 0.5268813 0.59785482 0.7657310 0.9523673 1.1391601
## FPE(n) 1.7069553 1.4145830 1.1397007 1.09324743 1.1559115 1.2459999 1.3442083
##                8
## AIC(n) 0.3104210
## HQ(n)  0.6984685
## SC(n)  1.2741092
## FPE(n) 1.3779586

AIC best at 4, BIC(SC) best at 3.

grangertest(tempdiff ~ pmm.timeseries, order = 3, data = tempdiff_pmm)
grangertest(pmm.timeseries ~ tempdiff, order = 3, data = tempdiff_pmm)
grangertest(tempdiff ~ pmm.timeseries, order = 4, data = tempdiff_pmm)
grangertest(pmm.timeseries ~ tempdiff, order = 4, data = tempdiff_pmm)
## Granger causality test
## 
## Model 1: tempdiff ~ Lags(tempdiff, 1:3) + Lags(pmm.timeseries, 1:3)
## Model 2: tempdiff ~ Lags(tempdiff, 1:3)
##   Res.Df Df      F  Pr(>F)  
## 1     85                    
## 2     88 -3 2.7238 0.04927 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Granger causality test
## 
## Model 1: pmm.timeseries ~ Lags(pmm.timeseries, 1:3) + Lags(tempdiff, 1:3)
## Model 2: pmm.timeseries ~ Lags(pmm.timeseries, 1:3)
##   Res.Df Df     F  Pr(>F)  
## 1     85                   
## 2     88 -3 2.817 0.04392 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Granger causality test
## 
## Model 1: tempdiff ~ Lags(tempdiff, 1:4) + Lags(pmm.timeseries, 1:4)
## Model 2: tempdiff ~ Lags(tempdiff, 1:4)
##   Res.Df Df      F Pr(>F)
## 1     82                 
## 2     86 -4 1.0453  0.389
## Granger causality test
## 
## Model 1: pmm.timeseries ~ Lags(pmm.timeseries, 1:4) + Lags(tempdiff, 1:4)
## Model 2: pmm.timeseries ~ Lags(pmm.timeseries, 1:4)
##   Res.Df Df      F  Pr(>F)  
## 1     82                    
## 2     86 -4 2.3089 0.06478 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

PMM lead Temp by 3; Temp lead PMM by 3 and 4. (實質上PMM領先Temp 4季;Temp領先PMM 2和3季)

Temp與PMM差分

temp_pmmdiff <- cbind(temp.timeseries, pmmdiff)
VARselect(y=na.omit(temp_pmmdiff), lag.max = 8, type = c("const"))
## $selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      3      3      2      3 
## 
## $criteria
##                1         2         3         4         5         6         7
## AIC(n) 0.4747269 0.2168313 0.1543781 0.1744728 0.1818065 0.2263362 0.2233073
## HQ(n)  0.5432059 0.3309629 0.3141624 0.3799097 0.4328961 0.5230785 0.5657021
## SC(n)  0.6447896 0.5002690 0.5511909 0.6846607 0.8053695 0.9632743 1.0736204
## FPE(n) 1.6076631 1.2424496 1.1677459 1.1923885 1.2026503 1.2596643 1.2589342
##                8
## AIC(n) 0.2540264
## HQ(n)  0.6420739
## SC(n)  1.2177146
## FPE(n) 1.3023997

AIC best at 3, BIC(SC) best at 2.

grangertest(temp.timeseries ~ pmmdiff, order = 2, data = temp_pmmdiff)
grangertest(pmmdiff ~ temp.timeseries, order = 2, data = temp_pmmdiff)
grangertest(temp.timeseries ~ pmmdiff, order = 3, data = temp_pmmdiff)
grangertest(pmmdiff ~ temp.timeseries, order = 3, data = temp_pmmdiff)
## Granger causality test
## 
## Model 1: temp.timeseries ~ Lags(temp.timeseries, 1:2) + Lags(pmmdiff, 1:2)
## Model 2: temp.timeseries ~ Lags(temp.timeseries, 1:2)
##   Res.Df Df      F  Pr(>F)  
## 1     88                    
## 2     90 -2 3.8896 0.02406 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Granger causality test
## 
## Model 1: pmmdiff ~ Lags(pmmdiff, 1:2) + Lags(temp.timeseries, 1:2)
## Model 2: pmmdiff ~ Lags(pmmdiff, 1:2)
##   Res.Df Df      F Pr(>F)
## 1     88                 
## 2     90 -2 0.1097 0.8962
## Granger causality test
## 
## Model 1: temp.timeseries ~ Lags(temp.timeseries, 1:3) + Lags(pmmdiff, 1:3)
## Model 2: temp.timeseries ~ Lags(temp.timeseries, 1:3)
##   Res.Df Df      F  Pr(>F)  
## 1     85                    
## 2     88 -3 2.2637 0.08689 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Granger causality test
## 
## Model 1: pmmdiff ~ Lags(pmmdiff, 1:3) + Lags(temp.timeseries, 1:3)
## Model 2: pmmdiff ~ Lags(pmmdiff, 1:3)
##   Res.Df Df      F Pr(>F)
## 1     85                 
## 2     88 -3 1.6811 0.1772

PMM lead Temp by 2 and 3. (實質上PMM領先Temp 3和4季)

Temp差分與PMM差分

tempdiff_pmmdiff <- cbind(tempdiff, pmmdiff)
VARselect(y=na.omit(tempdiff_pmmdiff), lag.max = 8, type = c("const"))
## $selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      4      3      2      4 
## 
## $criteria
##                1         2         3         4         5         6         7
## AIC(n) 0.6492568 0.2127901 0.1605420 0.1354522 0.1936255 0.2625013 0.2667114
## HQ(n)  0.7177358 0.3269217 0.3203263 0.3408892 0.4447150 0.5592436 0.6091063
## SC(n)  0.8193194 0.4962278 0.5573548 0.6456401 0.8171884 0.9994394 1.1170246
## FPE(n) 1.9142225 1.2374387 1.1749660 1.1467569 1.2169486 1.3060540 1.3147804
##                8
## AIC(n) 0.3140577
## HQ(n)  0.7021052
## SC(n)  1.2777459
## FPE(n) 1.3829789

AIC best at 4, BIC(SC) best at 2.

grangertest(tempdiff ~ pmmdiff, order = 2, data = tempdiff_pmmdiff)
grangertest(pmmdiff ~ tempdiff, order = 2, data = tempdiff_pmmdiff)
grangertest(tempdiff ~ pmmdiff, order = 4, data = tempdiff_pmmdiff)
grangertest(pmmdiff ~ tempdiff, order = 4, data = tempdiff_pmmdiff)
## Granger causality test
## 
## Model 1: tempdiff ~ Lags(tempdiff, 1:2) + Lags(pmmdiff, 1:2)
## Model 2: tempdiff ~ Lags(tempdiff, 1:2)
##   Res.Df Df      F  Pr(>F)  
## 1     88                    
## 2     90 -2 3.7053 0.02851 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Granger causality test
## 
## Model 1: pmmdiff ~ Lags(pmmdiff, 1:2) + Lags(tempdiff, 1:2)
## Model 2: pmmdiff ~ Lags(pmmdiff, 1:2)
##   Res.Df Df      F Pr(>F)
## 1     88                 
## 2     90 -2 2.0579 0.1338
## Granger causality test
## 
## Model 1: tempdiff ~ Lags(tempdiff, 1:4) + Lags(pmmdiff, 1:4)
## Model 2: tempdiff ~ Lags(tempdiff, 1:4)
##   Res.Df Df      F Pr(>F)
## 1     82                 
## 2     86 -4 1.0528 0.3853
## Granger causality test
## 
## Model 1: pmmdiff ~ Lags(pmmdiff, 1:4) + Lags(tempdiff, 1:4)
## Model 2: pmmdiff ~ Lags(pmmdiff, 1:4)
##   Res.Df Df      F  Pr(>F)  
## 1     82                    
## 2     86 -4 2.3456 0.06135 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

PMM 領先 Temp 2季; Temp 領先 PMM by 4季.

##Salt與PMM的滯後選擇 & Granger

Salt與PMM

salt_pmm <- cbind(salt.timeseries, pmm.timeseries)
VARselect(y=na.omit(salt_pmm), lag.max = 8, type = c("const"))
## $selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      4      1      1      4 
## 
## $criteria
##               1        2        3        4        5        6        7        8
## AIC(n) 1.793458 1.775214 1.764027 1.733818 1.804333 1.883607 1.932778 1.961071
## HQ(n)  1.866995 1.897776 1.935615 1.954430 2.073970 2.202269 2.300465 2.377782
## SC(n)  1.977463 2.081889 2.193373 2.285834 2.479019 2.680963 2.852804 3.003767
## FPE(n) 6.010691 5.903791 5.842007 5.674889 6.100847 6.622013 6.981580 7.217133

AIC best at 3, BIC(SC) best at 1.

grangertest(salt.timeseries ~ pmm.timeseries, order = 1, data = salt_pmm)
grangertest(pmm.timeseries ~ salt.timeseries, order = 1, data = salt_pmm)
grangertest(salt.timeseries ~ pmm.timeseries, order = 3, data = salt_pmm)
grangertest(pmm.timeseries ~ salt.timeseries, order = 3, data = salt_pmm)
## Granger causality test
## 
## Model 1: salt.timeseries ~ Lags(salt.timeseries, 1:1) + Lags(pmm.timeseries, 1:1)
## Model 2: salt.timeseries ~ Lags(salt.timeseries, 1:1)
##   Res.Df Df      F   Pr(>F)   
## 1     80                      
## 2     81 -1 8.6202 0.004339 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Granger causality test
## 
## Model 1: pmm.timeseries ~ Lags(pmm.timeseries, 1:1) + Lags(salt.timeseries, 1:1)
## Model 2: pmm.timeseries ~ Lags(pmm.timeseries, 1:1)
##   Res.Df Df      F Pr(>F)
## 1     80                 
## 2     81 -1 0.3919 0.5331
## Granger causality test
## 
## Model 1: salt.timeseries ~ Lags(salt.timeseries, 1:3) + Lags(pmm.timeseries, 1:3)
## Model 2: salt.timeseries ~ Lags(salt.timeseries, 1:3)
##   Res.Df Df      F Pr(>F)
## 1     74                 
## 2     77 -3 2.0146 0.1192
## Granger causality test
## 
## Model 1: pmm.timeseries ~ Lags(pmm.timeseries, 1:3) + Lags(salt.timeseries, 1:3)
## Model 2: pmm.timeseries ~ Lags(pmm.timeseries, 1:3)
##   Res.Df Df      F Pr(>F)
## 1     74                 
## 2     77 -3 0.4773  0.699

PMM 領先 salt 1季.