#################################################################################
# Regression Analysis                                                           #
#################################################################################

# Table 1 Main Paper
# Table 1 Supplementary Material
# Figure 8 Main Paper
# Figure 9 Main Paper
# Figure 10 Main Paper
# Figure 11 Main Paper

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

# remove 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                                                                       #
#################################################################################

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)
)

# Table 1: Difference between aggregated AIC Values and the minimum aggregated
# AIC Value for different Tradecorrelation Structures.
AIC_TAB<-round(AIC_TAB,1)
colnames(AIC_TAB)<-c("TC1","TC2","TC3","OLS")
stargazer(AIC_TAB,type="text")

#################################################################################
# Table 1 Supplementary Material                                                #
#################################################################################

# Supplementary Material:
# Table 1*: Difference between aggregated AIC Values and the minimum aggregated
# AIC Value for different Tradecorrelation Structures in the SEM and the SAR Model.                                                                     #



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")

# Now construct the AIC weights
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)")))

#################################################################################
# Figure 8                                                                      #
#################################################################################

# Figure 8: Annual Akaike Weights for different Specifcations (Raw Values
# in Gray, Smoothed in Black).

pdf("2) Estimating the Models/figure8.pdf",height = 7.5,width = 15)
par(mfrow=c(2,2),mar=c(4,8,4,4))
for (i in 1:4){
  plot(W[,i],type="l",ylim=c(0,1),col="lightgray",ylab="Akaike weights",xlab="t",xaxt="n",cex.axis=2,cex.lab=2,main=names[i],cex.main=2)
  lines(smooth.spline(W[,i],nknots = 65,spar=0.4),ylim=c(0,1),type="l",lwd=3,cex.axis=2)
  axis(1,1:65,1952:2016,cex.axis=2)
}
dev.off()


# open containers for plotting the coefficients
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()

# Extract the model output

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



# Now prepare the data for plotting

ts<-rho_errorsar_TC3
pv<-2*pnorm(-abs(ts/prho_errorsar_TC3),0,1)
Kiu<-rho_errorsar_TC3-2*prho_errorsar_TC3
Kio<-rho_errorsar_TC3+2*prho_errorsar_TC3

year<-1952:2016

pv[pv>0.1]<-3
pv[pv>0.05 & pv<0.1]<-2
pv[pv<0.05]<-1

data<-data.frame(year,ts, Kiu,Kio, pv)
colnames(data)<-c("Year","Rho","KIu","KIo", 'pi')
p1<-qplot(Year, Rho, data=data) + geom_line(color='black')+ggtitle(expression(rho))
p1<-p1+geom_ribbon(aes(ymin=KIu,ymax=KIo), alpha=0.3)
p1<-p1+scale_colour_manual(values=c('1'='springgreen4', '2'='orange1', '3'='red'))
p1<-p1+theme(legend.position='none', axis.text=element_text(size=10),axis.title=element_text(size=11),
             plot.title=element_text(size=14, face='bold'))+ geom_point(aes(colour=factor(pi)))

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


ts<-errorsar_b1
pv<-pvalerrorsar_b1
Kiu<-errorsar_b1-2*perrorsar_b1
Kio<-errorsar_b1+2*perrorsar_b1

year<-1952:2016

pv[pv>0.1]<-3
pv[pv>0.05 & pv<0.1]<-2
pv[pv<0.05]<-1

data<-data.frame(year,ts, Kiu,Kio, pv)
colnames(data)<-c("Year","Beta","KIu","KIo", 'pi')
p1<-qplot(Year, Beta, data=data) + geom_line(color='black')+ggtitle("Constant")
p1<-p1+geom_ribbon(aes(ymin=KIu,ymax=KIo), alpha=0.3)
p1<-p1+scale_colour_manual(values=c('1'='springgreen4', '2'='orange1', '3'='red'))
p1<-p1+theme(legend.position='none', axis.text=element_text(size=10),axis.title=element_text(size=11),
             plot.title=element_text(size=14, face='bold'))+ geom_point(aes(colour=factor(pi)))
p_2<-p1+geom_abline(intercept =0,slope=0)

ts<-errorsar_b2
pv<-pvalerrorsar_b2
Kiu<-errorsar_b2-2*perrorsar_b2
Kio<-errorsar_b2+2*perrorsar_b2

year<-1952:2016

pv[pv>0.1]<-3
pv[pv>0.05 & pv<0.1]<-2
pv[pv<0.05]<-1

data<-data.frame(year,ts, Kiu,Kio, pv)
colnames(data)<-c("Year","Beta","KIu","KIo", 'pi')
p1<-qplot(Year, Beta, data=data) + geom_line(color='black')+ggtitle("log(GDP), sender")
p1<-p1+geom_ribbon(aes(ymin=KIu,ymax=KIo), alpha=0.3)
p1<-p1+scale_colour_manual(values=c('1'='springgreen4', '2'='orange1', '3'='red'))
p1<-p1+theme(legend.position='none', axis.text=element_text(size=10),axis.title=element_text(size=11),
             plot.title=element_text(size=14, face='bold'))+ geom_point(aes(colour=factor(pi)))
p_3<-p1+geom_abline(intercept =0,slope=0)+ylim(-0.1,1.7)


ts<-errorsar_b3
pv<-pvalerrorsar_b3
Kiu<-errorsar_b3-2*perrorsar_b3
Kio<-errorsar_b3+2*perrorsar_b3

year<-1952:2016

pv[pv>0.1]<-3
pv[pv>0.05 & pv<0.1]<-2
pv[pv<0.05]<-1

data<-data.frame(year,ts, Kiu,Kio, pv)
colnames(data)<-c("Year","Beta","KIu","KIo", 'pi')
p1<-qplot(Year, Beta, data=data) + geom_line(color='black')+ggtitle("log(GDP), receiver")
p1<-p1+geom_ribbon(aes(ymin=KIu,ymax=KIo), alpha=0.3)
p1<-p1+scale_colour_manual(values=c('1'='springgreen4', '2'='orange1', '3'='red'))
p1<-p1+theme(legend.position='none', axis.text=element_text(size=10),axis.title=element_text(size=11),
             plot.title=element_text(size=14, face='bold'))+ geom_point(aes(colour=factor(pi)))
p_4<-p1+geom_abline(intercept =0,slope=0)+ylim(-0.1,1.7)

ts<-errorsar_b4
pv<-pvalerrorsar_b4
Kiu<-errorsar_b4-2*perrorsar_b4
Kio<-errorsar_b4+2*perrorsar_b4

year<-1952:2016

pv[pv>0.1]<-3
pv[pv>0.05 & pv<0.1]<-2
pv[pv<0.05]<-1

data<-data.frame(year,ts, Kiu,Kio, pv)
colnames(data)<-c("Year","Beta","KIu","KIo", 'pi')
p1<-qplot(Year, Beta, data=data) + geom_line(color='black')+ggtitle("Formal Alliance, Sender-Receiver")
p1<-p1+geom_ribbon(aes(ymin=KIu,ymax=KIo), alpha=0.3)
p1<-p1+scale_colour_manual(values=c('1'='springgreen4', '2'='orange1', '3'='red'))
p1<-p1+theme(legend.position='none', axis.text=element_text(size=10),axis.title=element_text(size=11),
             plot.title=element_text(size=14, face='bold'))+ geom_point(aes(colour=factor(pi)))
p_5<-p1+geom_abline(intercept =0,slope=0)


ts<-errorsar_b5
pv<-pvalerrorsar_b5
Kiu<-errorsar_b5-2*perrorsar_b5
Kio<-errorsar_b5+2*perrorsar_b5

year<-1952:2016

pv[pv>0.1]<-3
pv[pv>0.05 & pv<0.1]<-2
pv[pv<0.05]<-1

data<-data.frame(year,ts, Kiu,Kio, pv)
colnames(data)<-c("Year","Beta","KIu","KIo", 'pi')
p1<-qplot(Year, Beta, data=data) + geom_line(color='black')+ggtitle("Absolute Difference, Polity Score")
p1<-p1+geom_ribbon(aes(ymin=KIu,ymax=KIo), alpha=0.3)
p1<-p1+scale_colour_manual(values=c('1'='springgreen4', '2'='orange1', '3'='red'))
p1<-p1+theme(legend.position='none', axis.text=element_text(size=10),axis.title=element_text(size=11),
             plot.title=element_text(size=14, face='bold'))+ geom_point(aes(colour=factor(pi)))
p_6<-p1+geom_abline(intercept =0,slope=0)


ts<-errorsar_b6
pv<-pvalerrorsar_b6
Kiu<-errorsar_b6-2*perrorsar_b6
Kio<-errorsar_b6+2*perrorsar_b6

year<-1952:2016

pv[pv>0.1]<-3
pv[pv>0.05 & pv<0.1]<-2
pv[pv<0.05]<-1

data<-data.frame(year,ts, Kiu,Kio, pv)
colnames(data)<-c("Year","Beta","KIu","KIo", 'pi')
p1<-qplot(Year, Beta, data=data) + geom_line(color='black')+ggtitle("log(Distance), Sender-Receiver")
p1<-p1+geom_ribbon(aes(ymin=KIu,ymax=KIo), alpha=0.3)
p1<-p1+scale_colour_manual(values=c('1'='springgreen4', '2'='orange1', '3'='red'))
p1<-p1+theme(legend.position='none', axis.text=element_text(size=10),axis.title=element_text(size=11),
             plot.title=element_text(size=14, face='bold'))+ geom_point(aes(colour=factor(pi)))
p_7<-p1+geom_abline(intercept =0,slope=0)

#################################################################################
# Figure 9                                                                      #
#################################################################################

# Figure 9: Time Series of Coefcients from the Model with Tradecorrelation
# Structure N_sender-receiver. p-values by colours (green p < 0:05, orange p 
# <0:1, red p > 0:1). Shaded Regions give +- 2 standard deviations.

pdf("2) Estimating the Models/figure9.pdf", width = 12,height=12)
multiplot(p_3,p_7,p_5,p_4,p_6,p_1,cols=2)
dev.off()


rm(list=ls())
