-
Notifications
You must be signed in to change notification settings - Fork 1
/
AggregateGHG_CODE.R
102 lines (86 loc) · 4.94 KB
/
AggregateGHG_CODE.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
# AGGREGATED GHG CODE #
#############################################################################
req_packages<-c("devtools","glmnet","igraph","readxl","grid","gridExtra","Rmisc","ggplot2","sandwich","bootUR","mFilter","psych","parallel","scales","MASS")
install.packages(req_packages)
lapply(req_packages, require, character.only=TRUE)
devtools::install_github("Marga8/HDGCvar")
library(HDGCvar)
#load dataset and plot it
# AggregateGHG <- readRDS("own_path/AggregateGHG.rds")
DATA_R_clima2<-AggregateGHG[,-1] #-1 to remove the dates column
Dates<-c(1871:2018)
dataset_x_ggplot<-cbind(Dates, DATA_R_clima2)
plots_series<-list()
plots_series[[1]]<-ggplot(data = dataset_x_ggplot, aes(x = Dates, y = T))+ geom_line(color = "#00AFBB", size = 0.5) +
xlab("Years") + ylab("Temperature Anomalies")
plots_series[[2]]<-ggplot(data = dataset_x_ggplot, aes(x = Dates, y = G))+ geom_line(color = "#00AFBB", size = 0.5) +
xlab("Years") + ylab("GHG Concentration")
plots_series[[3]]<-ggplot(data = dataset_x_ggplot, aes(x = Dates, y = S))+ geom_line(color = "#00AFBB", size = 0.5) +
xlab("Years") + ylab("Solar Activity")
plots_series[[4]]<-ggplot(data = dataset_x_ggplot, aes(x = Dates, y = V))+ geom_line(color = "#00AFBB", size = 0.5) +
xlab("Years") + ylab("Stratospheric Aerosols")
plots_series[[5]]<-ggplot(data = dataset_x_ggplot, aes(x = Dates, y = Y))+ geom_line(color = "#00AFBB", size = 0.5) +
xlab("Years") + ylab("GDP (log 2010 USD)")
plots_series[[6]]<-ggplot(data = dataset_x_ggplot, aes(x = Dates, y = A))+ geom_line(color = "#00AFBB", size = 0.5) +
xlab("Years") + ylab("Tropospheric Aerosols")
plots_series[[7]]<-ggplot(data = dataset_x_ggplot, aes(x = Dates, y = N))+ geom_line(color = "#00AFBB", size = 0.5) +
xlab("Years") + ylab("El Nino SOI")
plots_series[[8]]<-ggplot(data = dataset_x_ggplot, aes(x = Dates, y = O))+ geom_line(color = "#00AFBB", size = 0.5) +
xlab("Years") + ylab("Ocean Heat Content")
grid.arrange(grobs=plots_series, ncol=2)
#lag-length selection (VAR order)
selected_lag<-HDGCvar::lags_upbound_BIC(DATA_R_clima2,p_max=10)
print(selected_lag)
#compute network
network<-HDGCvar::HDGC_VAR_all(DATA_R_clima2, p = selected_lag, d = 2, bound = 0.5 * nrow(DATA_R_clima2),
parallel = TRUE, n_cores = NULL)
#plot network
HDGCvar::Plot_GC_all(network, Stat_type="FS_cor",alpha=0.10, multip_corr=list(T,"none"),
directed=T, layout=layout.circle, main="Climate Network",edge.arrow.size=.2,
vertex.size=16, vertex.color=c("lightblue"), vertex.frame.color="blue",
vertex.label.size=4,vertex.label.color="black",vertex.label.cex=0.8,
vertex.label.dist=c(0,0,0,0,0,0,0,0), edge.curved=0,cluster=list(F,5,"black",0.8,1,0))
## HEAT MAP table for p-values ##
##########################################################################
pvals<-as.matrix(network[["tests"]][,,2,2])
pv_max <- 0.15
trunc_pvals <- pmin(pvals, pv_max)
series <-colnames(DATA_R_clima2)
GC_df <- data.frame(pvalue = c(trunc_pvals),
series_to = factor(rep(series, ncol(DATA_R_clima2)), levels = series),
series_from = factor(rep(series, each = ncol(DATA_R_clima2)), levels = series))
ggplot(GC_df, aes(x = series_from, y = series_to, fill = pvalue)) +
geom_tile(colour = "grey50") +
scale_y_discrete(limits = rev(levels(GC_df$series_from))) +
labs(x = "Granger causality from", y = "Granger causality to", fill = "p-value") +
scale_fill_gradient2(low = "darkblue", mid = "blue", high = "white", midpoint = pv_max / 2) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
axis.text.x = element_text(angle = 0, size = 9),
axis.text.y = element_text(size = 9),
legend.position = "none")
#highlight paths
##########################################################################
input<-as.matrix(network[["tests"]][,,2,2])#p_values
input<-t(input)
input[input < 0.1] <- 1 #put =1 values < alpha
input[is.na(input)] <- 0 #put =0 the diagonal
input[input != 1] <- 0 #put =0 values > alpha
my_graph=graph_from_adjacency_matrix(input, mode='directed',diag=F,add.rownames = TRUE )
V(my_graph)$label = rownames(input)
# Preferences
FROM_V <- 4 #put =3 for GDP(Y) or adjust as preferred
TO_V <- 6
# Calculate all simple paths from FROM_V to TO_V as list of vertecy sequences
graph_outlet <- all_simple_paths(my_graph,from=FROM_V,to=TO_V)
graph_outlet
V(my_graph)$color <- "white"
V(my_graph)$color[unique(unlist(graph_outlet))] <- "lightblue"
V(my_graph)$color[c(FROM_V,TO_V)] <- "yellow"
# Colour each of the paths
E(my_graph)$color <- "gray"
E(my_graph)$width<- 1
lapply(graph_outlet, function(x) E(my_graph, path=x)$color <<- "red")
lapply(graph_outlet, function(x) E(my_graph, path=x)$width <<- 3)
# Plot all paths and mark the group of vertecies through which paths flow
plot(my_graph,layout=layout.circle,main="Climate Network",edge.arrow.size=.5,
vertex.size=c(16,16,16,25,16,25,16,16),vertex.label.size=4)