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

# Figure 1 Supplementary Material
# Figure 2 Supplementary Material

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

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


# load packages
library(spdep)
library(pspline)
library(mvtnorm)
library(stargazer)
library(maxLik)
library(psych)
library(ggplot2)
library(reshape)
library(gridExtra)


# Set working Directory
setwd("")

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

# open containers for the results
res_errorsar_TC3<-list()
res_own<-list()

for (t in 1952:2016){
  
  # Extract the Response of interest
  Y_use<-log(X[which(X[,2]==t),which(colnames(X)=="EX_ij")])
  # Extract the covariates of interest
  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")])
  
  # Extract the weighting matrix of interest
  W_all<-TC1[[t-1951]]
  # Exclude self-loops
  diag(W_all)<-0
  # Row-normalize W_i
  group<-diag(1/rowSums(W_all))
  group[group== Inf]<-0
  W_all<-group%*%W_all
  
  # Include all three weighting matrices
  
  W_9<-TC1[[t-1951]]
  diag(W_9)<-0
  group<-diag(1/rowSums(W_9))
  group[group== Inf]<-0
  W_9<-group%*%W_9
  
  W_10<-TC2[[t-1951]]
  diag(W_10)<-0
  group<-diag(1/rowSums(W_10))
  group[group== Inf]<-0
  W_10<-group%*%W_10
  
  W_11<-TC3[[t-1951]]
  diag(W_11)<-0
  group<-diag(1/rowSums(W_11))
  group[group== Inf]<-0
  W_11<-group%*%W_11
  
  # Find starting values by employing the errorSAR model
  
  errorsarlm_TC1<-errorsarlm(Y_use~-1+X_use,listw = mat2listw(W_all),zero.policy=TRUE)
  
  # Defin the log-likelihood as a function of rho and sigma^2
  
  logl_1<-function(input1){
    beta<-coef(errorsarlm_TC1)[-1]
    rho<-input1[1:3]
    sigma2<-input1[4]
    A<-diag(rep(1,length(Y_use)))-rho[1]*W_9-rho[2]*W_10-rho[3]*W_11
    
    res<- -(length(Y_use)/2)*log(pi*sigma2)+ log(sum(eval(A)))-(1/(2*sigma2))*(t(A%*%(Y_use-X_use%*%beta))%*%A%*%(Y_use-X_use%*%beta))
    return(res)
  }
  

  # maximize with respect to rho  
  res_1<-maxLik(logl_1,start=c(0,0,0,errorsarlm_TC1$s2))
  
  # Reweight the results in order to finde beta
  A_est<-diag(rep(1,length(Y_use)))-res_1$estimate[1]*W_9-res_1$estimate[2]*W_10-res_1$estimate[3]*W_11
  W_est<-res_1$estimate[1]*W_9+res_1$estimate[2]*W_10+res_1$estimate[3]*W_11
  A_XX<-t(X_use)%*%X_use-t(X_use)%*%W_est%*%X_use-t(X_use)%*%t(W_est)%*%X_use+t(X_use)%*%t(W_est)%*%W_est%*%X_use
  A_XY<-t(X_use)%*%Y_use-t(X_use)%*%W_est%*%Y_use-t(X_use)%*%t(W_est)%*%Y_use+t(X_use)%*%t(W_est)%*%W_est%*%Y_use
  beta_new<-solve(A_XX)%*%A_XY
  
  # Define the log-likelihood as a function of beta, rho and sigma^2
  logl_2<-function(input1){
    beta<-input1[1:6]
    rho<-input1[7:9]
    sigma2<-input1[10]
    A<-diag(rep(1,length(Y_use)))-rho[1]*W_9-rho[2]*W_10-rho[3]*W_11
    
    res<- -(length(Y_use)/2)*log(pi*sigma2)+ log(sum(eval(A)))-(1/(2*sigma2))*(t(A%*%(Y_use-X_use%*%beta))%*%A%*%(Y_use-X_use%*%beta))
    return(res)
  }
  
  # Define the gradient of the log-likelihood as a function of beta, rho and sigma^2
  gradient_logl_2<-function(input1){
    beta<-input1[1:6]
    rho<-input1[7:9]
    sigma2<-input1[10]
    A<-diag(rep(1,length(Y_use)))-rho[1]*W_9-rho[2]*W_10-rho[3]*W_11
    Q_1<--W_9-t(W_9)+2*rho[1]*t(W_9)%*%W_9+rho[2]*(t(W_9)%*%W_10+t(W_10)%*%W_9)+rho[3]*(t(W_9)%*%W_11+t(W_11)%*%W_9)
    Q_2<--W_10-t(W_10)+2*rho[2]*t(W_10)%*%W_10+rho[1]*(t(W_9)%*%W_10+t(W_10)%*%W_9)+rho[3]*(t(W_10)%*%W_11+t(W_11)%*%W_10)
    Q_3<--W_11-t(W_11)+2*rho[3]*t(W_11)%*%W_11+rho[1]*(t(W_11)%*%W_9+t(W_11)%*%W_9)+rho[2]*(t(W_10)%*%W_11+t(W_11)%*%W_10)
    comp1<-(1/sigma2)*(t(A%*%X_use)%*%(A%*%Y_use-A%*%X_use%*%beta))
    comp_rho1<--tr(solve(A)%*%W_9)-(1/sigma2)*( t(Y_use)%*%Q_1%*%Y_use-2*t(beta)%*%t(X_use)%*%Q_1%*%Y_use+t(beta)%*%t(X_use)%*%Q_1%*%X_use%*%beta )
    comp_rho2<--tr(solve(A)%*%W_10)-(1/sigma2)*( t(Y_use)%*%Q_2%*%Y_use-2*t(beta)%*%t(X_use)%*%Q_2%*%Y_use+t(beta)%*%t(X_use)%*%Q_2%*%X_use%*%beta )
    comp_rho3<--tr(solve(A)%*%W_11)-(1/sigma2)*( t(Y_use)%*%Q_3%*%Y_use-2*t(beta)%*%t(X_use)%*%Q_3%*%Y_use+t(beta)%*%t(X_use)%*%Q_3%*%X_use%*%beta )
    comp_sigma<--length(Y_use)/(2*sigma2)+1/(2*sigma2^2)*(t((A%*%(Y_use-X_use%*%beta)))%*%(A%*%(Y_use-X_use%*%beta)))
    
    return<-c(comp1,comp_rho1,comp_rho2,comp_rho3,comp_sigma)
  }
  
  # Define the Hessian of the log-likelihood as a function of beta, rho and sigma^2
  hessian_logl_2<-function(input1){
    beta<-input1[1:6]
    rho<-input1[7:9]
    sigma2<-input1[10]
    A<-diag(rep(1,length(Y_use)))-rho[1]*W_9-rho[2]*W_10-rho[3]*W_11
    Sigma<-solve(A)%*%solve(t(A))*sigma2
    
    Q_1<--W_9-t(W_9)+2*rho[1]*t(W_9)%*%W_9+rho[2]*(t(W_9)%*%W_10+t(W_10)%*%W_9)+rho[3]*(t(W_9)%*%W_11+t(W_11)%*%W_9)
    Q_2<--W_10-t(W_10)+2*rho[2]*t(W_10)%*%W_10+rho[1]*(t(W_9)%*%W_10+t(W_10)%*%W_9)+rho[3]*(t(W_10)%*%W_11+t(W_11)%*%W_10)
    Q_3<--W_11-t(W_11)+2*rho[3]*t(W_11)%*%W_11+rho[1]*(t(W_11)%*%W_9+t(W_11)%*%W_9)+rho[2]*(t(W_10)%*%W_11+t(W_11)%*%W_10)
   
    #beta
    dell_beta_beta<- - 1/(2*sigma2)*(t(A%*%X_use)%*%(A%*%X_use))
    dell_beta_rho1<-rep(0,6)
    dell_beta_rho2<-rep(0,6)
    dell_beta_rho3<-rep(0,6)
    dell_beta_sigma<-rep(0,6)
    H_beta<-cbind(dell_beta_beta,dell_beta_rho1,dell_beta_rho2,dell_beta_rho3,dell_beta_sigma)
    
    # rho
    dell_rho1_rho1<--tr(solve(A)%*%W_9%*%solve(A)%*%W_9) - 1/(sigma2)*(tr(t(W_9)%*%W_9%*%Sigma))
    dell_rho1_rho2<--tr(solve(A)%*%W_9%*%solve(A)%*%W_10) - 1/(2*sigma2)*(tr((t(W_9)%*%W_10+t(W_10)%*%W_9)%*%Sigma))
    dell_rho1_rho3<--tr(solve(A)%*%W_9%*%solve(A)%*%W_11) -1/(2*sigma2)*(tr((t(W_9)%*%W_11+t(W_11)%*%W_9)%*%Sigma))
    dell_rho2_rho1<-dell_rho1_rho2
    dell_rho2_rho2<--tr(solve(A)%*%W_10%*%solve(A)%*%W_10) - 1/(sigma2)*(tr(t(W_10)%*%W_10%*%Sigma))
    dell_rho2_rho3<--tr(solve(A)%*%W_10%*%solve(A)%*%W_11) - 1/(2*sigma2)*(tr((t(W_11)%*%W_10+t(W_10)%*%W_11)%*%Sigma))
    dell_rho3_rho1<-dell_rho1_rho3
    dell_rho3_rho2<-dell_rho2_rho3
    dell_rho3_rho3<--tr(solve(A)%*%W_11%*%solve(A)%*%W_11) - 1/(sigma2)*(tr(t(W_11)%*%W_11%*%Sigma))
    dell_rho1_sigma<--(1/(2*sigma2^2))*( tr(Q_1%*%Sigma) )
    dell_rho2_sigma<--(1/(2*sigma2^2))*( tr(Q_2%*%Sigma) )
    dell_rho3_sigma<--(1/(2*sigma2^2))*( tr(Q_3%*%Sigma))
    
    H_rho1<-c(dell_beta_rho1,dell_rho1_rho1,dell_rho1_rho2,dell_rho1_rho3,dell_rho1_sigma)
    H_rho2<-c(dell_beta_rho2,dell_rho2_rho1,dell_rho2_rho2,dell_rho2_rho3,dell_rho2_sigma)
    H_rho3<-c(dell_beta_rho3,dell_rho3_rho1,dell_rho3_rho2,dell_rho3_rho3,dell_rho3_sigma)
    
    
    # sigma
    dell_sigma_beta<-dell_beta_sigma
    dell_sigma_rho1<-dell_rho1_sigma
    dell_sigma_rho2<-dell_rho2_sigma
    dell_sigma_rho3<-dell_rho3_sigma
    dell_sigma_sigma<- -length(Y_use)/(2*sigma2^2)
    H_sigma<-c(dell_sigma_beta,dell_sigma_rho1,dell_sigma_rho2,dell_sigma_rho3,dell_sigma_sigma)
    
    hessian<-rbind(H_beta,H_rho1,H_rho2,H_rho3,H_sigma)
    return(hessian)
    }
  
  # Optimize the log-likelihood
  res_2<-maxLik(logl_2,start=c(beta_new,res_1$estimate),grad=gradient_logl_2,hess=hessian_logl_2)
  # Evaluate the Hessian
  hessian<-hessian_logl_2(res_2$estimate)
  # Check whether the Hessian is invertible and whether the 
  # Fisher Information has non-negative entries
  print(diag(solve(-hessian)))

  # Save the results of the procedure
 res_own[[t-1951]]<-cbind(res_2$estimate,sqrt(diag(solve(-hessian))))
 res_errorsar_TC3[[t-1951]]<-errorsarlm_TC1
 
    print(t)
  
}

# Save the results of the estimation procedure in containers

coef_own<-c()
coef_sar<-c()
sd_sar<-c()
p_sar<-c()
sd_own<-c()
p_own<-c()
for (t in 1:65){
  coef_own<-cbind(coef_own,res_own[[t]][,1])
  sd_own<-cbind(sd_own,res_own[[t]][,2])
  coef_sar<-cbind(coef_sar,coef(res_errorsar_TC3[[t]]))
  p_sar<-cbind(p_sar,summary(res_errorsar_TC3[[t]])$ LR$p.value)
  sd_sar<-cbind(sd_sar,summary(res_errorsar_TC3[[t]])$ lambda.se)
  p_own<-cbind(p_own,2*pnorm(-abs(res_own[[t]][,2]/sd_own[[t]])))
}


# Preparte the time series for plotting

ts<-as.numeric(coef_sar[1,])
pv<-(1-pnorm(ts/sd_sar[1,],0,1))
Kiu<-c(ts-2*sd_sar)
Kio<-c(ts+2*sd_sar)


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("Tradecorrelation 3")
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)+ylim(low=-1, high=1.2)
p_1
###


ts<-coef_own[7,]
pv<-2*pnorm(-abs(ts/sd_own[7,]))
Kiu<-c(ts-2*sd_own[7,])
Kio<-c(ts+2*sd_own[7,])

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_1","KIu","KIo", 'pi')
p1<-qplot(Year, Rho_1, data=data) + geom_line(color='black')+ggtitle("Tradecorrelation, 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_2<-p1+geom_abline(intercept =0,slope=0)+ylim(low=-1.4, high=2)
p_2

ts<-coef_own[8,]
pv<-2*pnorm(-abs(ts/sd_own[8,]))
Kiu<-c(ts-2*sd_own[8,])
Kio<-c(ts+2*sd_own[8,])

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_2","KIu","KIo", 'pi')
p1<-qplot(Year, Rho_2, data=data) + geom_line(color='black')+ggtitle("Tradecorrelation, 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_3<-p1+geom_abline(intercept =0,slope=0)+ylim(low=-1.4, high=2)
p_3

ts<-coef_own[9,]
pv<-2*pnorm(-abs(ts/sd_own[9,]))
Kiu<-c(ts-2*sd_own[9,])
Kio<-c(ts+2*sd_own[9,])

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_3","KIu","KIo", 'pi')
p1<-qplot(Year, Rho_3, data=data) + geom_line(color='black')+ggtitle("Tradecorrelation, 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_4<-p1+geom_abline(intercept =0,slope=0)+ylim(low=-1.4, high=2)
p_4


#################################################################################
# Figure 2 - Supplementary Material                                             #
#################################################################################


pdf("3) Supplementary Material/figure2_SM.pdf", width = 12,height=6)
multiplot(p_2,p_3,p_4,cols=3)
dev.off()



###### Now the coefficients


ts<-coef_own[1,]
pv<-2*pnorm(-abs(ts/sd_own[1,]))
Kiu<-c(ts-2*sd_own[1,])
Kio<-c(ts+2*sd_own[1,])

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_1<-p1+geom_abline(intercept =0,slope=0)#+ylim(low=-1.4, high=2)
p_1

ts<-coef_own[2,]
pv<-2*pnorm(-abs(ts/sd_own[2,]))
Kiu<-c(ts-2*sd_own[2,])
Kio<-c(ts+2*sd_own[2,])

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_2<-p1+geom_abline(intercept =0,slope=0)+ylim(low=-2, high=2)
p_2

ts<-coef_own[3,]
pv<-2*pnorm(-abs(ts/sd_own[3,]))
Kiu<-c(ts-2*sd_own[3,])
Kio<-c(ts+2*sd_own[3,])

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_3<-p1+geom_abline(intercept =0,slope=0)+ylim(low=-2, high=2)
p_3



ts<-coef_own[4,]
pv<-2*pnorm(-abs(ts/sd_own[4,]))
Kiu<-c(ts-2*sd_own[4,])
Kio<-c(ts+2*sd_own[4,])

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_4<-p1+geom_abline(intercept =0,slope=0)+ylim(low=-1.4, high=2)
p_4

ts<-coef_own[5,]
pv<-2*pnorm(-abs(ts/sd_own[5,]))
Kiu<-c(ts-2*sd_own[5,])
Kio<-c(ts+2*sd_own[5,])

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_5<-p1+geom_abline(intercept =0,slope=0)
p_5


ts<-coef_own[6,]
pv<-2*pnorm(-abs(ts/sd_own[6,]))
Kiu<-c(ts-2*sd_own[6,])
Kio<-c(ts+2*sd_own[6,])

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_6<-p1+geom_abline(intercept =0,slope=0)
p_6

#################################################################################
# Figure 1 - Supplementary Material                                             #
#################################################################################

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






