#################################################################################
#################################################################################
# Code for the Supplementary Material                                           #
#################################################################################

# Figure 3 Supplementary Material


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

# rm the old stuff
rm(list=ls())


# load packages
library(spdep)
library(pspline)
library(stargazer)
library(ggplot2)
library(car)

# Set working Directory
setwd("")

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

#################################################################################
# Run the regressions for all tradecorrelation structures                       #
#################################################################################

## TC_1
res_errorsar_TC1<-list()
res_lagsar_TC1<-list()
res_lm<-list()

for (t in 1952:2016){
  print(t)
  Y_use<-log(X[which(X[,2]==t),which(colnames(X)=="EX_ij")])
  X_use<-cbind(1,
               X[which(X[,2]==t),which(colnames(X)=="lGDP_i")],
               X[which(X[,2]==t),which(colnames(X)=="lGDP_j")],
               X[which(X[,2]==t),which(colnames(X)=="DA_ij")],
               X[which(X[,2]==t),which(colnames(X)=="POL_ij")],
               X[which(X[,2]==t),which(colnames(X)=="D_ij")])
  
  W_all<-TC1[[t-1951]]
  diag(W_all)<-0
  group<-diag(1/rowSums(W_all))
  group[group== Inf]<-0
  W_all<-group%*%W_all
  res_errorsar_TC1[[t-1951]]<-errorsarlm(Y_use~-1+X_use,listw = mat2listw(W_all),zero.policy=TRUE)
  res_lagsar_TC1[[t-1951]]<-lagsarlm(Y_use~-1+X_use,listw = mat2listw(W_all),zero.policy=TRUE)
  res_lm[[t-1951]]<-lm(Y_use~-1+X_use)
  
}


# Calculate the AIC
AIC_lagsar_TC1<-c()
AIC_errorsar_TC1<-c()
AIC_lm<-c()
for (i in 1:65){
  AIC_errorsar_TC1<-c(AIC_errorsar_TC1,AIC(res_errorsar_TC1[[i]]))
  AIC_lagsar_TC1<-c(AIC_lagsar_TC1,AIC(res_lagsar_TC1[[i]]))
  AIC_lm<-c(AIC_lm,AIC(res_lm[[i]]))
  print(i)
}

## TC_2
res_errorsar_TC2<-list()
res_lagsar_TC2<-list()

for (t in 1952:2016){
  print(t)
  Y_use<-log(X[which(X[,2]==t),which(colnames(X)=="EX_ij")])
  X_use<-cbind(1,
               X[which(X[,2]==t),which(colnames(X)=="lGDP_i")],
               X[which(X[,2]==t),which(colnames(X)=="lGDP_j")],
               X[which(X[,2]==t),which(colnames(X)=="DA_ij")],
               X[which(X[,2]==t),which(colnames(X)=="POL_ij")],
               X[which(X[,2]==t),which(colnames(X)=="D_ij")])
  
  W_all<-TC2[[t-1951]]
  diag(W_all)<-0
  group<-diag(1/rowSums(W_all))
  group[group== Inf]<-0
  W_all<-group%*%W_all
  res_lagsar_TC2[[t-1951]]<-lagsarlm(Y_use~-1+X_use,listw = mat2listw(W_all),zero.policy=TRUE)
  res_errorsar_TC2[[t-1951]]<-errorsarlm(Y_use~-1+X_use,listw = mat2listw(W_all),zero.policy=TRUE)
}


# Calculate the AIC
AIC_errorsar_TC2<-c()
AIC_lagsar_TC2<-c()

for (i in 1:65){
  AIC_errorsar_TC2<-c(AIC_errorsar_TC2,AIC(res_errorsar_TC2[[i]]))
  AIC_lagsar_TC2<-c(AIC_lagsar_TC2,AIC(res_lagsar_TC2[[i]]))
}

## TC3
res_errorsar_TC3<-list()
res_lagsar_TC3<-list()

for (t in 1952:2016){
  print(t)
  Y_use<-log(X[which(X[,2]==t),which(colnames(X)=="EX_ij")])
  X_use<-cbind(1,
               X[which(X[,2]==t),which(colnames(X)=="lGDP_i")],
               X[which(X[,2]==t),which(colnames(X)=="lGDP_j")],
               X[which(X[,2]==t),which(colnames(X)=="DA_ij")],
               X[which(X[,2]==t),which(colnames(X)=="POL_ij")],
               X[which(X[,2]==t),which(colnames(X)=="D_ij")])
  
  W_all<-TC3[[t-1951]]
  diag(W_all)<-0
  group<-diag(1/rowSums(W_all))
  group[group== Inf]<-0
  W_all<-group%*%W_all
  res_errorsar_TC3[[t-1951]]<-errorsarlm(Y_use~-1+X_use,listw = mat2listw(W_all),zero.policy=TRUE)
  res_lagsar_TC3[[t-1951]]<-lagsarlm(Y_use~-1+X_use,listw = mat2listw(W_all),zero.policy=TRUE)
}


# Calculate the AIC
AIC_errorsar_TC3<-c()
AIC_lagsar_TC3<-c()

for (i in 1:65){
  AIC_errorsar_TC3<-c(AIC_errorsar_TC3,AIC(res_errorsar_TC3[[i]]))
  AIC_lagsar_TC3<-c(AIC_lagsar_TC3,AIC(res_lagsar_TC3[[i]]))
}



### Table 1: Difference between aggregated AIC Values and the minimum aggregated AIC Value for different Tradecorrelation Structures

AIC_TAB<-cbind(sum(AIC_errorsar_TC1)-sum(AIC_errorsar_TC3),
               
               sum(AIC_errorsar_TC2)-sum(AIC_errorsar_TC3),
               
               sum(AIC_errorsar_TC3)-sum(AIC_errorsar_TC3),
               
               sum(AIC_lm)-sum(AIC_errorsar_TC3)
)


AIC_TAB<-round(AIC_TAB,1)
colnames(AIC_TAB)<-c("TC1","TC2","TC3","OLS")
stargazer(AIC_TAB,type="latex")


AIC_TAB2<-rbind(cbind(sum(AIC_errorsar_TC1)-sum(AIC_errorsar_TC3),
               
               sum(AIC_errorsar_TC2)-sum(AIC_errorsar_TC3),
               
               sum(AIC_errorsar_TC3)-sum(AIC_errorsar_TC3),sum(AIC_lm)-sum(AIC_errorsar_TC3)),
               
               cbind(sum(AIC_lagsar_TC1)-sum(AIC_errorsar_TC3),
               sum(AIC_lagsar_TC2)-sum(AIC_errorsar_TC3),
               sum(AIC_lagsar_TC3)-sum(AIC_errorsar_TC3),sum(AIC_lm)-sum(AIC_errorsar_TC3))
               

)


AIC_TAB2<-round(AIC_TAB2,1)
colnames(AIC_TAB2)<-c("TC1","TC2","TC3","OLS")
stargazer(AIC_TAB2,type="text")


X<-cbind(AIC_errorsar_TC1,AIC_errorsar_TC2,AIC_errorsar_TC3,AIC_lm)
colnames(X)<-c("TC1","TC2","TC3","OLS")
min<-c()

for (i in 1:65){
  min<-c(min, which(X[i,]==min(X[i,])))
}
minval<-c()

for (i in 1:65){
  minval<-c(minval, X[i,which(X[i,]==min(X[i,]))])
}
Delta_i<- X-minval

norm<-c()
for (i in 1:65){
  norm<-c(norm,sum(exp(-0.5*Delta_i[i,])))
}

W<-c()

for (i in 1:65){
  col<-c()
  for (j in 1:4){
    col<-c(col,exp(-0.5*Delta_i[i,j])/norm[i])
  }
  W<-rbind(W,col)
}



names<-c(expression(paste("Tradcorrelation, Sender (",N["sender"],")")),expression(paste("Tradcorrelation, Receiver (",N["receiver"],")")),expression(paste("Tradcorrelation, Sender-Receiver (",N["sender-receiver"],")")),expression(paste("OLS (",rho,"=0)")))


# Now we re-estimate the Stuff with the Akaike-Weights applied

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

res_errorsar_weighted<-list()

for (t in 1952:2016){
  
  Y_use<-log(X[which(X[,2]==t),which(colnames(X)=="EX_ij")])
  X_use<-cbind(1,
               X[which(X[,2]==t),which(colnames(X)=="lGDP_i")],
               X[which(X[,2]==t),which(colnames(X)=="lGDP_j")],
               X[which(X[,2]==t),which(colnames(X)=="DA_ij")],
               X[which(X[,2]==t),which(colnames(X)=="POL_ij")],
               log(1+X[which(X[,2]==t),which(colnames(X)=="D_ij")]))

  
  
  W_all<-W[t-1951,1]*TC1[[t-1951]]+W[t-1951,2]*TC2[[t-1951]]+W[t-1951,3]*TC3[[t-1951]]
  diag(W_all)<-0
  group<-diag(1/rowSums(W_all))
  group[group== Inf]<-0
  W_all<-group%*%W_all
  
  
  res_errorsar_weighted[[t-1951]]<-errorsarlm(Y_use~-1+X_use,listw = mat2listw(W_all),zero.policy=TRUE)
print(t)
  
}




#################################################################################
#################################################################################
# Plot the Coefficients                                                         #
#################################################################################


# open containers 
rho_errorsar_TC3<-c()
prho_errorsar_TC3<-c()
pvalrho_errorsar_TC3<-c()

errorsar_b1<-c()
perrorsar_b1<-c()
pvalerrorsar_b1<-c()
errorsar_b2<-c()
perrorsar_b2<-c()
pvalerrorsar_b2<-c()
errorsar_b3<-c()
perrorsar_b3<-c()
pvalerrorsar_b3<-c()
errorsar_b4<-c()
perrorsar_b4<-c()
pvalerrorsar_b4<-c()
errorsar_b5<-c()
perrorsar_b5<-c()
pvalerrorsar_b5<-c()
errorsar_b6<-c()
perrorsar_b6<-c()
pvalerrorsar_b6<-c()
errorsar_b7<-c()
perrorsar_b7<-c()
pvalerrorsar_b7<-c()
errorsar_b8<-c()
perrorsar_b8<-c()
pvalerrorsar_b8<-c()


for (t in 1952:2016){
  
  
  rho_errorsar_TC3<-c(rho_errorsar_TC3,res_errorsar_weighted[[t-1951]]$lambda)
  prho_errorsar_TC3<-c(prho_errorsar_TC3,summary(res_errorsar_weighted[[t-1951]])$ lambda.se)
  pvalrho_errorsar_TC3<-c(pvalrho_errorsar_TC3,summary(res_errorsar_weighted[[t-1951]])$ LR$p.value)
  
  errorsar_b1<-c(errorsar_b1,summary(res_errorsar_weighted[[t-1951]])$Coef[1,1])
  perrorsar_b1<-c(perrorsar_b1,summary(res_errorsar_weighted[[t-1951]])$Coef[1,2])
  pvalerrorsar_b1<-c(pvalerrorsar_b1,summary(res_errorsar_weighted[[t-1951]])$Coef[1,4])
  
  errorsar_b2<-c(errorsar_b2,summary(res_errorsar_weighted[[t-1951]])$Coef[2,1])
  perrorsar_b2<-c(perrorsar_b2,summary(res_errorsar_weighted[[t-1951]])$Coef[2,2])
  pvalerrorsar_b2<-c(pvalerrorsar_b2,summary(res_errorsar_weighted[[t-1951]])$Coef[2,4])
  
  errorsar_b3<-c(errorsar_b3,summary(res_errorsar_weighted[[t-1951]])$Coef[3,1])
  perrorsar_b3<-c(perrorsar_b3,summary(res_errorsar_weighted[[t-1951]])$Coef[3,2])
  pvalerrorsar_b3<-c(pvalerrorsar_b3,summary(res_errorsar_weighted[[t-1951]])$Coef[3,4])
  
  errorsar_b4<-c(errorsar_b4,summary(res_errorsar_weighted[[t-1951]])$Coef[4,1])
  perrorsar_b4<-c(perrorsar_b4,summary(res_errorsar_weighted[[t-1951]])$Coef[4,2])
  pvalerrorsar_b4<-c(pvalerrorsar_b4,summary(res_errorsar_weighted[[t-1951]])$Coef[4,4])
  
  errorsar_b5<-c(errorsar_b5,summary(res_errorsar_weighted[[t-1951]])$Coef[5,1])
  perrorsar_b5<-c(perrorsar_b5,summary(res_errorsar_weighted[[t-1951]])$Coef[5,2])
  pvalerrorsar_b5<-c(pvalerrorsar_b5,summary(res_errorsar_weighted[[t-1951]])$Coef[5,4])
  
  errorsar_b6<-c(errorsar_b6,summary(res_errorsar_weighted[[t-1951]])$Coef[6,1])
  perrorsar_b6<-c(perrorsar_b6,summary(res_errorsar_weighted[[t-1951]])$Coef[6,2])
  pvalerrorsar_b6<-c(pvalerrorsar_b6,summary(res_errorsar_weighted[[t-1951]])$Coef[6,4])

}



# Now prepare the data

ts<-rho_errorsar_TC3
year<-1952:2016
data<-data.frame(year,ts)
colnames(data)<-c("Year","Rho")
p1<-qplot(Year, Rho, data=data) + geom_line(color='black')+ggtitle(expression(rho))
p_1<-p1+geom_abline(intercept =0,slope=0)
###


ts<-errorsar_b1
data<-data.frame(year,ts)
colnames(data)<-c("Year","Beta")
p1<-qplot(Year, Beta, data=data) + geom_line(color='black')+ggtitle("Constant")

p_2<-p1+geom_abline(intercept =0,slope=0)

ts<-errorsar_b2
data<-data.frame(year,ts)
colnames(data)<-c("Year","Beta")
p1<-qplot(Year, Beta, data=data) + geom_line(color='black')+ggtitle("log(GDP), sender")
p_3<-p1+geom_abline(intercept =0,slope=0)+ylim(-0.1,1.7)


ts<-errorsar_b3
data<-data.frame(year,ts)
colnames(data)<-c("Year","Beta")
p1<-qplot(Year, Beta, data=data) + geom_line(color='black')+ggtitle("log(GDP), receiver")
p_4<-p1+geom_abline(intercept =0,slope=0)+ylim(-0.1,1.7)

ts<-errorsar_b4
data<-data.frame(year,ts)
colnames(data)<-c("Year","Beta")
p1<-qplot(Year, Beta, data=data) + geom_line(color='black')+ggtitle("Formal Alliance, Sender-Receiver")
p_5<-p1+geom_abline(intercept =0,slope=0)


ts<-errorsar_b5
data<-data.frame(year,ts)
colnames(data)<-c("Year","Beta")
p1<-qplot(Year, Beta, data=data) + geom_line(color='black')+ggtitle("Absolute Difference, Polity Score")
p_6<-p1+geom_abline(intercept =0,slope=0)


ts<-errorsar_b6

data<-data.frame(year,ts)
colnames(data)<-c("Year","Beta")
p1<-qplot(Year, Beta, data=data) + geom_line(color='black')+ggtitle("log(Distance), Sender-Receiver")
p_7<-p1+geom_abline(intercept =0,slope=0)


#################################################################################
# Figure 3 - Supplementary Material                                             #
#################################################################################

pdf("3) Supplementary Material/figure3_SM.pdf", width = 12,height=12)
multiplot(p_3,p_7,p_5,p_4,p_6,p_1,cols=2)
dev.off()


