#################################################################################
#################################################################################
# Creating descriptive figures                                                  #
#################################################################################

# Figure 1 Main Paper
# Figure 2 Main Paper
# Figure 3 Main Paper
# Figure 4 Main Paper
# Figure 5 Main Paper
# Figure 6 Main Paper
# Figure 7 Main Paper

#################################################################################
# Head                                                                          #
#################################################################################

# remove the old Stuff
rm(list=ls())

# load some packages
library("statnet")
library("ggplot2")
library("igraph")

# Set working Directory
setwd("")

# Load the Data
load("Data.RData")


# Aggregate the export volume in each year for the
# countries included in the model
tiv_all<-c()
for (i in 1:67){
  tiv_all<-c(tiv_all, sum(relevant[[i]]))
}


# Calculate the export volume of the United States ...
share_us<-c()
# ... the Soviet Union ...
share_su<-c()
# ... and Russia.
share_ru<-c()

for (i in 1:67){
  select<-relevant[[i]]
  share_us<-c(share_us,rowSums(select)[which(names=="United States")])
  share_su<-c(share_su,rowSums(select)[which(names=="Soviet Union")])
  share_ru<-c(share_ru,rowSums(select)[which(names=="Russia")])
}

# Calculate the remaining export volume
rest<-tiv_all-share_us-share_su-share_ru

# Combine the export volume of the Soviet Union and Russia
share_su_ru<-share_su+share_ru

# Prepare plotting
d <- data.frame(year=rep(1950:2016,each=3),countries=rep(c("Soviet Union/Russia","United States","other"),67),TIV=c(rbind(share_su_ru,share_us,rest)))
d$countries<-factor(d$countries,c("other","United States","Soviet Union/Russia"))
p<-ggplot(d, aes(x=year,y=TIV,group=countries,fill=countries)) + geom_area( position = "stack")

p1<-p+  scale_fill_grey(start=0, end=0.85)
p1<-p1+  theme_bw()+ theme(legend.position="top")
Density<-c()
for ( i in 1:67){
  net_bin<-amk[[i]][nEX[,i]==1,nEX[,i]==1]
  net_bin[net_bin>0]<-1
  Density<-c(Density,sum(net_bin)/(dim(net_bin)[1]*(dim(net_bin)[1]-1)))
}


d <- data.frame(year=rep(1950:2016,each=1),Density=c(Density))

p2<-ggplot(d, aes(x=year,y=Density)) +geom_line()
p2<-p2+ theme_bw()

#################################################################################
# Figure 1                                                                      #
#################################################################################

# Create Figure 1: Development of the exported TIV for 1950-2016 for the United
# States, Soviet Union/Russia and all other Countries.
pdf("1) Drawing descriptive Figures/figure1.pdf",height =3,width = 6)
p1
dev.off()


#################################################################################
# Figure 3                                                                      #
#################################################################################

# Create Figure 3: Development of the Network Density (Share of realized Trade
# Relations as compared to all possible Trade Relations) for 1950-2016.
pdf("1) Drawing descriptive Figures/figure3.pdf",height =3,width = 6)
p2
dev.off()



# Extract the outdegree ...
outdeg<-c()
# ... the exports ...
export<-c()
# ... the indegree and ...
indeg<-c()
# ... the imports.
import<-c()
for (i in 1:67){
  select<-amk[[i]][1:224,1:224]
  select_bin<-select
  select_bin[select_bin>0]<-1
  
  outdeg<-c(outdeg,rowSums(select_bin))
  indeg<-c(indeg,colSums(select_bin))
  
  
  import<-c(import,colSums(select))
  export<-c(export,rowSums(select))
  
}


# Extract means... 
export_means<-c()
import_means<-c()
# ... minimum values ... 
export_min<-c()
import_min<-c()
# ... and maximum values.
export_max<-c()
import_max<-c()

for (i in 1: max(outdeg)){
  if (length(which(outdeg==i))>0){
  export_means<-c(export_means,mean(export[which(outdeg==i)]))
  export_min<-c(export_min,min(export[which(outdeg==i)]))
  export_max<-c(export_max,max(export[which(outdeg==i)]))
  }
  
  if (length(which(outdeg==i))==0){
    export_means<-c(export_means,0)
    export_min<-c(export_min,0)
    export_max<-c(export_max,0)
  }
  
  if (length(which(indeg==i))>0){
  import_means<-c(import_means,mean(import[which(indeg==i)]))
  import_min<-c(import_min,min(import[which(indeg==i)]))
  import_max<-c(import_max,max(import[which(indeg==i)]))
  }
  
  if (length(which(indeg==i))==0){
    import_means<-c(import_means,0)
    import_min<-c(import_min,0)
    import_max<-c(import_max,0)
  }
  
}

#################################################################################
# Figure 2                                                                      #
#################################################################################

# Create Figure 2: Log-Log Plot of TIV Export vs. Outdegree (left) and TIV Import
# vs. Indegree (right). The solid Line in the Middle gives the Mean of Exports
# or Imports associated with the respective In- or Outdegree. Whiskers
# indicate minimum and maximum Values.

pdf("1) Drawing descriptive Figures/figure2.pdf",height =8,width = 16)

par(mfrow=c(1,2),mar=c(5,5,5,5))
plot(1:67,export_means,log="xy",type="b",ylim=c(0.01,20000),xlab="Outdegree",ylab="Aggregated Export",cex.main=3,cex.lab=2,cex.axis=2)
arrows(1:67, export_min, 1:67, export_max, length=0.05, angle=90, code=3,col="gray")
lines(1:67,export_min,col="gray")
lines(1:67,export_max,col="gray")


plot(1:16,import_means[1:16],log="xy",type="b",ylim=c(0.01,20000),xlab="Indegree",ylab="Aggregated Import",cex.main=3,cex.lab=2,cex.axis=2)
arrows(1:16, import_min[1:16], 1:16, import_max[1:16], length=0.05, angle=90, code=3,col="gray")
lines(1:16,import_min[1:16],col="gray")
lines(1:16,import_max[1:16],col="gray")

dev.off()
# Do not be afraid of the warning messages, they are
# simply related to the fact that we do not observe
# all outdegree/indegree combinations


# We look at the network in 1952
t=1952
# Scale the nodes
cex<-1+log(1+amk[[t-1949]][1:224,1:224])
# Scale the edges
net_50<-amk[[t-1949]][1:224,1:224]*0.2
cex_nodes<-2+log(1+rowSums(net_50))*0.3
# Give Names
colnames(net_50)<-laenderliste[1:224,5]
rownames(net_50)<-laenderliste[1:224,5]
# Define it as a network
net_50<-network(net_50)

#################################################################################
# Figure 4                                                                      #
#################################################################################

# Create Figure 4: Directed valued Network of Arms Trade in 1952. Edges (Nodes)
# are scaled proportional to the logarithmic Exports (logarithmic aggregated
# Exports) measured in TIV.
pdf("1) Drawing descriptive Figures/figure4.pdf",height =20,width = 10)
plot(net_50,displayisolates=F,displaylabels=T,label.cex=0.8,main=paste("Arms trade network in ",t),edge.lwd=cex,vertex.col="gray",edge.col="gray",vertex.cex=cex_nodes,label.pos=5,interactive=F)
dev.off()

# This is a force-directed algorithm and hence the network
# representation changes with each plot

# Define the log-TIV values for transfers in t=1952
Y_use<-log(X[which(X[,2]%in%t),which(colnames(X)%in%"EX_ij")])
# Define the Sending and Receiving Set for t=1952
R_use<-R[which(X[,2]%in%t),]
# Gives names
R_use<-paste(laenderliste[R_use[,1],5],"-",laenderliste[R_use[,2],5],sep="")

## Encode the different dependence structures

# TC1
W_1<-TC1[[t-1951]]
# TC2
W_2<-TC2[[t-1951]]
# TC 3
W_3<-TC3[[t-1951]]

#################################################################################
# Figure 5                                                                      #
#################################################################################

# Create Figure 5: Tradecorrelation Structure N_sender.
# Colouring according to logarithmic TIV, ranging from Yellow (low) to Red (high).

# Remove transfers that are not connected
rS<-rowSums(W_1)
W_kill<-W_1[-which(rS==0),-which(rS==0)]
Y_use1<-Y_use[-which(rS==0)]
R_use1<-R_use[-which(rS==0)]
g<-graph_from_adjacency_matrix(W_kill,mode="undirected")
col_Y<-colorRampPalette(c('yellow', 'red'))(length(Y_use1))[rank(Y_use1)]
tkplot(g, vertex.color=col_Y,vertex.size=10,vertex.label.family = "sans",vertex.label.color="black",layout=layout.kamada.kawai(g),vertex.label=R_use1)

#################################################################################
# Figure 6                                                                      #
#################################################################################

# Create Figure 6: Tradecorrelation Structure N_receiver.
# Colouring according to logarithmic TIV values, ranging from Yellow (low) to Red (high).

rS<-rowSums(W_2)
W_kill<-W_2[-which(rS==0),-which(rS==0)]
Y_use2<-Y_use[-which(rS==0)]
R_use2<-R_use[-which(rS==0)]
g<-graph_from_adjacency_matrix(W_kill,mode="undirected")
col_Y<-colorRampPalette(c('yellow', 'red'))(length(Y_use2))[rank(Y_use2)]
tkplot(g, vertex.color=col_Y,vertex.size=10,vertex.label.family = "sans",vertex.label.color="black",layout=layout.fruchterman.reingold(g),vertex.label=R_use2)

#################################################################################
# Figure 7                                                                      #
#################################################################################

# Figure 7: Tradecorrelation Structure Nsender-receiver. Colouring according
# to logarithmic TIV, ranging from Yellow (low) to Red (high).
rS<-rowSums(W_3)
R_use3<-R_use
W_kill_norm<-diag(1/rS)%*%W_3
g<-graph_from_adjacency_matrix(W_3,mode="undirected")
col_Y<-colorRampPalette(c('yellow', 'red'))(length(Y_use))[rank(Y_use)]
tkplot(g, vertex.color=col_Y,vertex.size=10,vertex.label.family = "sans",vertex.label.color="black",layout=layout.kamada.kawai(g),vertex.label=R_use3)

# remove the old Stuff
rm(list=ls())
