############################################################################################
#############           TIME-SERIES ANALYSES WITH GRANGER'S METHOD         #################
############################################################################################

#### 1. Packages
library(tseries)
library(Hmisc)
library(psych)
library(corrplot)
library(lmtest)
library(vars)
library(urca)


#### 2. Files and data
setwd("...")
data <- read.csv2("Pelletier2021_BrSDivR.csv")
as.data.frame(data)
data$DivR[is.na(data$DivR)]<- 0.2224972 #removing NA and adding DivR mean
data$DivR <- as.numeric(data$DivR)
data$BrS <- as.numeric(data$BrS)
names(data)[1] <- "Yr"
str(data)
sapply(data, mode)
attach(data)

#### Coefficient of variation (CV)
sapply(data, function(x) sd(x) / mean(x) * 100)


#### 3. Plot time series
par(mfrow=c(2,1))
plot(Yr, BrS, ty ="l", main = "Breeding success")
abline(h = mean(BrS), lty = 3)
abline(lm(BrS~Yr), col="gray", lty=5)

plot(Yr, DivR, ty ="l", main = "Divorce rate")
abline(h = mean(DivR), lty = 3)
abline(lm(DivR~Yr), col="gray", lty=5)

dev.off()


#### 4. AUtocorrelation function
par(mfrow=c(2,1))
acf(BrS, lag.max = 5, main = "Breeding success", ylab = "Autocorrelation")
acf(DivR, lag.max = 5, main = "Divorce rate", ylab = "Autocorrelation")
dev.off()


#### 5. Detrend time series (we can use also diff() function in 'tseries' package,
####  but diff() remove decreases sample of 1...
reg1 = lm(BrS~Yr)
summary(reg1) # STATIONARY
R1 = reg1$residuals
BrS_R = rstandard(lm(BrS~Yr))

reg2 = lm(DivR~Yr)
summary(reg2) # STATIONARY
R2 = reg2$residuals
DivR_R = rstandard(lm(DivR~Yr))


#### 6. Descriptive statistics
psych::describe(data)


#### 7. Stationarity tests
adf1 = summary(ur.df(BrS, type = "trend", selectlag = "AIC"))
print(adf1) # STATIONARY CONFIRMED
kpss.test(BrS, null = "Trend")

adf2 = summary(ur.df(DivR, type = "trend", selectlag = "AIC"))
print(adf2) # STATIONARY CONFIRMED
kpss.test(DivR, null = "Trend")


#### 8. Cross-correlation function
par(mfrow=c(2,1))
ccf(BrS, DivR, lag.max = 5, main = "Breeding success & Divorce rate", ylab = "Cross-correlation",lwd = 3)
ccf(DivR, BrS, lag.max = 5, main = "Divorce rate & Breeding success", ylab = "Cross-correlation",lwd = 3)
dev.off()


#### 9. Bivariate Vector autoregressive regression (VAR) analyses

### 9.1 Combine variables two-by-two
dat1 = cbind(BrS, DivR)

### 9.2 Lag order selection (with the lowest AIC value)
VARselect(dat1, lag.max = 2, type = "const") #lag = 1

### 9.3 VAR estimation
x1 = VAR(dat1, p = 1, type = "const")
summary(x1)
x1r <- restrict(x1)
summary(x1r)

### 9.4 Diagnostic testing (only for large number of lags and series with ~100 obs.)
serial.test(x1r, lags.pt = 3, type = "PT.adjusted")

### 9.5 Granger test
causality(x1r, cause = "BrS") #BrS Granger-cause DivR data (with restrict())
causality(x1r, cause = "DivR") #DivR do not Granger-cause BrS data (with restrict())


### 9.6 Impulse response analysis
plot(irf(x1r, cumulative = FALSE, impulse = 'BrS', response = 'DivR',  n.ahead = 5, boot = TRUE, runs = 1000, ci = 0.95))

IRA1r <- irf(x1r, cumulative = FALSE, impulse = 'BrS', response = 'DivR',  n.ahead = 5, boot = TRUE, runs = 1000, ci = 0.95)
IRA1r

### 9.7 Forecast error variance decomposition
plot(fevd(x1r, n.ahead = 3), nsteps = 3)

fevd1r <- fevd(x1r, n.ahead = 3)
fevd1r