#################################################################################
# Simulation                                                                    #
#################################################################################

# Table 2 Main Paper

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

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

# load packages
library(spdep)
library(MASS)
library(stargazer)
# Set working Directory
setwd("")

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


# Define a function that creates the tradecorrelation structures
tradecorr_create<-function(ajdM,Ind_structure="str1",Returntype="binary"){
  
  #################################################################################################
  # This function uses an Adjacency Matrix as Input and Creates several Outputs that represent ####
  # the tradecorrelation structure.         #######################################################
  #################################################################################################
  # ajdM ##########################################################################################
  # - Symmetric Adjacency Matrices that represent a graph. It is possible that there exist binary #
  #   positive valued entries                                                                     #
  #                                                                                               #
  # tradcorr structure ############################################################################
  # - Argument that speciefies how the tradecorrelation should look:                              #
  #   * str1: N_1 (Y_ij depends on all other exports of i and on the export of j to i)            #
  #   * str2: N_2 (Y_ij Y_ij depends on all other imports of j and a reciprocial trade (j,i))     #
  #   * str3: N_3 (Y_ij depends on all other imports and exports of i and j)                      #
  # Returntype ####################################################################################
  # - Argument that specifies the Return of the function                                          #
  #   * list: Retunrs a list that shows the cond. independence structure                          #
  #   * binary: Returs a binary matrix representing the independence structure                    #
  #   * weighted: Returns a matrix with rowsums normed to one                                     #
  #                                                                                               #
  #################################################################################################
  
  
  
  big_N<-list()
  if (Ind_structure=="str1"){
    for (i in 1:dim(A)[1]) {
      for (j in 1:dim(A)[1]){
        if (A[i,j]>0) {
          
          name<-paste("Y",i,j,sep="_")
          big_N[[name]]<-paste("Y",i,j,sep="_")
          if (j<dim(A)[1]) {
            for (k in (j+1):dim(A)[1]){
              if (A[i,k]>0){
                
                big_N[[name]]<-append(big_N[[name]],paste("Y",i,k,sep="_"))
                
              }
            }
          }
          if (j>1) {
            for (l in 1:(j-1)){
              if (A[i,l]>0){
                
                big_N[[name]]<-append(big_N[[name]],paste("Y",i,l,sep="_"))
              }
            }
          }
          
          if (A[j,i]>0) {
            big_N[[name]]<-append(big_N[[name]],paste("Y",j,i,sep="_"))
          }
        }
      }
    }
    
    
  }
  
  
  # End N_1
  
  if (Ind_structure=="str2"){
    
    for (i in 1:dim(A)[1]) {
      for (j in 1:dim(A)[1]){
        if (A[i,j]>0) {
          
          name<-paste("Y",i,j,sep="_")
          big_N[[name]]<-paste("Y",i,j,sep="_")
          if (A[j,i]>0){
            big_N[[name]]<-append(big_N[[name]],paste("Y",j,i,sep="_"))
          }
          
          for (l in 1:dim(A)[1]){
            if (A[l,j]>0){
              big_N[[name]]<-append(big_N[[name]],paste("Y",l,j,sep="_"))
            }
          }
        }
      }
    }
    
  }
  
  
  # End N_2
  
  if (Ind_structure=="str3"){
    
    for (i in 1:dim(A)[1]) {
      for (j in 1:dim(A)[1]){
        if (A[i,j]>0) {
          
          name<-paste("Y",i,j,sep="_")
          big_N[[name]]<-paste("Y",i,j,sep="_")
          for (k in 1:dim(A)[1]){
            if (A[k,j]>0){
              big_N[[name]]<-append(big_N[[name]],paste("Y",k,j,sep="_"))
            }
          }
          for (l in 1:dim(A)[1]){
            if (A[j,l]>0){
              big_N[[name]]<-append(big_N[[name]],paste("Y",j,l,sep="_"))
            }
          }
          
          for (m in 1:dim(A)[1]){
            if (A[i,m]>0){
              big_N[[name]]<-append(big_N[[name]],paste("Y",i,m,sep="_"))
            }
          }
          for (n in 1:dim(A)[1]){
            if (A[n,i]>0){
              big_N[[name]]<-append(big_N[[name]],paste("Y",n,i,sep="_"))
            }
          }
        }
      }
    }
    
  }
  
  
  len<-length(big_N)
  
  W<-matrix(rep(0,len^2),ncol=len,nrow=len)
  
  for (i in 1:len) {
    for (j in 1:len) {
      if (names(big_N)[j] %in% big_N[[i]]){
        
        W[i,j]<-1}
    }
  }
  
  for (i in 1:len) {
    for (j in 1:len) {
      if (W[i,j]==1){
        
        W[j,i]<-1
      }
    }
  }
  
  diag(W)<-0
  W_weight<-W
  rs<-rowSums(W)
  for (i in 1:len){
    for (j in 1:len) {
      if (W[i,j]==1){
        W_weight[i,j]<-1/rs[i]
      }
    }
  }
  
  
  if (Returntype=="list"){
    return(big_N)
  }
  
  if (Returntype=="binary"){
    return(W)
  }
  
  if (Returntype=="weighted"){
    return(W_weight)
  }
  
}

#################################################################################
# Set the Parameters                                                            #
#################################################################################

years <-1952:2016
N<-100
k=3 
constant<-5
beta<-c(constant,rep(1,k))
rho<-0.5
sigma2<-1

#################################################################################
# Open some containers                                                          #
#################################################################################

res_11<-c()
res_12<-c()
res_13<-c()
res_14<-c()
res_21<-c()
res_22<-c()
res_23<-c()
res_24<-c()
res_31<-c()
res_32<-c()
res_33<-c()
res_34<-c()
res_41<-c()
res_42<-c()
res_43<-c()
res_44<-c()

aic_11<-c()
aic_12<-c()
aic_13<-c()
aic_14<-c()
aic_21<-c()
aic_22<-c()
aic_23<-c()
aic_24<-c()
aic_31<-c()
aic_32<-c()
aic_33<-c()
aic_34<-c()
aic_41<-c()
aic_42<-c()
aic_43<-c()
aic_44<-c()

for (rep in 1:N){
  print(rep)
  for (t in years){
    tau<-t-1949
    # starting network is the real network from the respective year
    A<-amk[[tau]][nEX[,tau]==1,nEX[,tau]==1]
    # binarize the network
    A[A!=0]<-1
    
    
    X<-cbind(1,matrix(c(rnorm(sum(A)*k,0,1)),nrow=sum(A),ncol=k))
    
    ##1
    W_all<-tradecorr_create(A,Ind_structure = "str1")
    diag(W_all)<-0
    group<-diag(1/rowSums(W_all))
    group[group== Inf]<-0
    W<-group%*%W_all
    
    S=sigma2*solve(diag(rep(1,sum(A)))-rho*W )%*%solve(diag(rep(1,sum(A)))-rho*t(W) )
    Y<-mvrnorm(n=1,mu=X%*%beta,Sigma=S)
    
    #11
    W_all<-tradecorr_create(A,Ind_structure = "str1")
    diag(W_all)<-0
    group<-diag(1/rowSums(W_all))
    group[group== Inf]<-0
    W<-group%*%W_all
    
    est<-errorsarlm(Y~X,listw = mat2listw(W),zero.policy=TRUE)
    res_11<-cbind(res_11,coef(est))
    aic_11<-c(aic_11,AIC(est))
    #12
    W_all<-tradecorr_create(A,Ind_structure = "str2")
    diag(W_all)<-0
    group<-diag(1/rowSums(W_all))
    group[group== Inf]<-0
    W<-group%*%W_all
    
    est<-errorsarlm(Y~X,listw = mat2listw(W),zero.policy=TRUE)
    res_12<-cbind(res_12,coef(est))
    aic_12<-c(aic_12,AIC(est))
    
    #13
    W_all<-tradecorr_create(A,Ind_structure = "str3")
    diag(W_all)<-0
    group<-diag(1/rowSums(W_all))
    group[group== Inf]<-0
    W<-group%*%W_all
    
    est<-errorsarlm(Y~X,listw = mat2listw(W),zero.policy=TRUE)
    res_13<-cbind(res_13,coef(est))
    aic_13<-c(aic_13,AIC(est))
    
    #14
    est<-lm(Y~X)
    aic_14<-c(aic_14,AIC(est))
    
    
    ##2
    W_all<-tradecorr_create(A,Ind_structure = "str2")
    diag(W_all)<-0
    group<-diag(1/rowSums(W_all))
    group[group== Inf]<-0
    W<-group%*%W_all
    
    S=sigma2*solve(diag(rep(1,sum(A)))-rho*W )%*%solve(diag(rep(1,sum(A)))-rho*t(W) )
    Y<-mvrnorm(n=1,mu=X%*%beta,Sigma=S)
    
    #21
    W_all<-tradecorr_create(A,Ind_structure = "str1")
    diag(W_all)<-0
    group<-diag(1/rowSums(W_all))
    group[group== Inf]<-0
    W<-group%*%W_all
    
    est<-errorsarlm(Y~X,listw = mat2listw(W),zero.policy=TRUE)
    res_21<-cbind(res_21,coef(est))
    aic_21<-c(aic_21,AIC(est))
    
    #22
    W_all<-tradecorr_create(A,Ind_structure = "str2")
    diag(W_all)<-0
    group<-diag(1/rowSums(W_all))
    group[group== Inf]<-0
    W<-group%*%W_all
    
    est<-errorsarlm(Y~X,listw = mat2listw(W),zero.policy=TRUE)
    res_22<-cbind(res_22,coef(est))
    aic_22<-c(aic_22,AIC(est))
    
    #23
    W_all<-tradecorr_create(A,Ind_structure = "str3")
    diag(W_all)<-0
    group<-diag(1/rowSums(W_all))
    group[group== Inf]<-0
    W<-group%*%W_all
    
    est<-errorsarlm(Y~X,listw = mat2listw(W),zero.policy=TRUE)
    res_23<-cbind(res_23,coef(est))
    aic_23<-c(aic_23,AIC(est))
    
    #24
    
    est<-lm(Y~X)
    aic_24<-c(aic_24,AIC(est))
    
    
    ##3
    W_all<-tradecorr_create(A,Ind_structure = "str3")
    diag(W_all)<-0
    group<-diag(1/rowSums(W_all))
    group[group== Inf]<-0
    W<-group%*%W_all
    
    S=sigma2*solve(diag(rep(1,sum(A)))-rho*W )%*%solve(diag(rep(1,sum(A)))-rho*t(W) )
    Y<-mvrnorm(n=1,mu=X%*%beta,Sigma=S)
    
    #31
    W_all<-tradecorr_create(A,Ind_structure = "str1")
    diag(W_all)<-0
    group<-diag(1/rowSums(W_all))
    group[group== Inf]<-0
    W<-group%*%W_all
    
    est<-errorsarlm(Y~X,listw = mat2listw(W),zero.policy=TRUE)
    res_31<-cbind(res_31,coef(est))
    aic_31<-c(aic_31,AIC(est))
    
    #32
    W_all<-tradecorr_create(A,Ind_structure = "str2")
    diag(W_all)<-0
    group<-diag(1/rowSums(W_all))
    group[group== Inf]<-0
    W<-group%*%W_all
    
    est<-errorsarlm(Y~X,listw = mat2listw(W),zero.policy=TRUE)
    res_32<-cbind(res_32,coef(est))
    aic_32<-c(aic_32,AIC(est))
    
    #33
    W_all<-tradecorr_create(A,Ind_structure = "str3")
    diag(W_all)<-0
    group<-diag(1/rowSums(W_all))
    group[group== Inf]<-0
    W<-group%*%W_all
    
    est<-errorsarlm(Y~X,listw = mat2listw(W),zero.policy=TRUE)
    res_33<-cbind(res_33,coef(est))
    aic_33<-c(aic_33,AIC(est))
    
    #34
    est<-lm(Y~X)
    aic_34<-c(aic_34,AIC(est))
    
    
    
    ##4
    W_all<-tradecorr_create(A,Ind_structure = "str3")
    
    S=sigma2*diag(rep(1,dim(W_all)[1]))
    Y<-mvrnorm(n=1,mu=X%*%beta,Sigma=S)
    
    #41
    W_all<-tradecorr_create(A,Ind_structure = "str1")
    diag(W_all)<-0
    group<-diag(1/rowSums(W_all))
    group[group== Inf]<-0
    W<-group%*%W_all
    
    est<-errorsarlm(Y~X,listw = mat2listw(W),zero.policy=TRUE)
    res_41<-cbind(res_41,coef(est))
    aic_41<-c(aic_41,AIC(est))
    
    #42
    W_all<-tradecorr_create(A,Ind_structure = "str2")
    diag(W_all)<-0
    group<-diag(1/rowSums(W_all))
    group[group== Inf]<-0
    W<-group%*%W_all
    
    est<-errorsarlm(Y~X,listw = mat2listw(W),zero.policy=TRUE)
    res_42<-cbind(res_42,coef(est))
    aic_42<-c(aic_42,AIC(est))
    
    #43
    W_all<-tradecorr_create(A,Ind_structure = "str3")
    diag(W_all)<-0
    group<-diag(1/rowSums(W_all))
    group[group== Inf]<-0
    W<-group%*%W_all
    
    est<-errorsarlm(Y~X,listw = mat2listw(W),zero.policy=TRUE)
    res_43<-cbind(res_43,coef(est))
    aic_43<-c(aic_43,AIC(est))
    
    #44
    est<-lm(Y~X)
    aic_44<-c(aic_44,AIC(est))
    
    
    
    
    
    
    
    print(t)
  }
  
}



#################################################################################
# Figure ?? - Not Included in the Paper showing the KDE Plot                                         #
#################################################################################

pdf("4) Simulations/kernel_density.pdf")
par(mfrow=c(4,3))
plot(density(res_11[1,]),ylab="DGP1",xlim=c(-1,1),main="Model 1",xlab="")
abline(v=rho,col="red")
abline(v=mean(res_11[1,]),col="gray",lty=2)
abline(v=quantile(res_11[1,],0.05),col="gray",lty=1)
abline(v=quantile(res_11[1,],0.95),col="gray",lty=1)

plot(density(res_12[1,]),ylab="",xlim=c(-1,1),main="Model 2",xlab="")
abline(v=rho,col="red")
abline(v=mean(res_12[1,]),col="gray",lty=2)
abline(v=quantile(res_12[1,],0.05),col="gray",lty=1)
abline(v=quantile(res_12[1,],0.95),col="gray",lty=1)

plot(density(res_13[1,]),ylab="",xlim=c(-1,1),main="Model 3",xlab="")
abline(v=rho,col="red")
abline(v=mean(res_13[1,]),col="gray",lty=2)
abline(v=quantile(res_13[1,],0.05),col="gray",lty=1)
abline(v=quantile(res_13[1,],0.95),col="gray",lty=1)

plot(density(res_21[1,]),ylab="DGP2",xlim=c(-1,1),main="",xlab="")
abline(v=rho,col="red")
abline(v=mean(res_21[1,]),col="gray",lty=2)
abline(v=quantile(res_21[1,],0.05),col="gray",lty=1)
abline(v=quantile(res_21[1,],0.95),col="gray",lty=1)

plot(density(res_22[1,]),ylab="",xlim=c(-1,1),main="",xlab="")
abline(v=rho,col="red")
abline(v=mean(res_22[1,]),col="gray",lty=2)
abline(v=quantile(res_22[1,],0.05),col="gray",lty=1)
abline(v=quantile(res_22[1,],0.95),col="gray",lty=1)

plot(density(res_23[1,]),ylab="",xlim=c(-1,1),main="",xlab="")
abline(v=rho,col="red")
abline(v=mean(res_23[1,]),col="gray",lty=2)
abline(v=quantile(res_23[1,],0.05),col="gray",lty=1)
abline(v=quantile(res_23[1,],0.95),col="gray",lty=1)

plot(density(res_31[1,]),ylab="DGP3",xlim=c(-1,1),main="",xlab="")
abline(v=rho,col="red")
abline(v=mean(res_31[1,]),col="gray",lty=2)
abline(v=quantile(res_31[1,],0.05),col="gray",lty=1)
abline(v=quantile(res_31[1,],0.95),col="gray",lty=1)

plot(density(res_32[1,]),ylab="",xlim=c(-1,1),main="",xlab="")
abline(v=rho,col="red")
abline(v=mean(res_32[1,]),col="gray",lty=2)
abline(v=quantile(res_32[1,],0.05),col="gray",lty=1)
abline(v=quantile(res_32[1,],0.95),col="gray",lty=1)

plot(density(res_33[1,]),ylab="",xlim=c(-1,1),main="",xlab="")
abline(v=rho,col="red")
abline(v=mean(res_33[1,]),col="gray",lty=2)
abline(v=quantile(res_33[1,],0.05),col="gray",lty=1)
abline(v=quantile(res_33[1,],0.95),col="gray",lty=1)
plot(density(res_41[1,]),ylab="DGP4",xlim=c(-1,1),main="",xlab="")
abline(v=0,col="red")
abline(v=mean(res_41[1,]),col="gray",lty=2)
abline(v=quantile(res_41[1,],0.05),col="gray",lty=1)
abline(v=quantile(res_41[1,],0.95),col="gray",lty=1)

plot(density(res_42[1,]),ylab="",xlim=c(-1,1),main="",xlab="")
abline(v=0,col="red")
abline(v=mean(res_42[1,]),col="gray",lty=2)
abline(v=quantile(res_42[1,],0.05),col="gray",lty=1)
abline(v=quantile(res_42[1,],0.95),col="gray",lty=1)

plot(density(res_43[1,]),ylab="",xlim=c(-1,1),main="",xlab="")
abline(v=0,col="red")
abline(v=mean(res_43[1,]),col="gray",lty=2)
abline(v=quantile(res_43[1,],0.05),col="gray",lty=1)
abline(v=quantile(res_43[1,],0.95),col="gray",lty=1)
dev.off()


## AIC table

min_model<-c()
aic<-rbind(aic_11,aic_12,aic_13,aic_14)
min<-apply(aic,2,min)
for (i in 1:length(min)){
  min_model<-c(min_model,which(aic[,i]==min[i]))
}


DGP1<-cbind(  c(mean(aic_11),median(aic_11),length(which(min_model == 1))/length(min)),
              c(mean(aic_12),median(aic_12),length(which(min_model == 2))/length(min)),
              c(mean(aic_13),median(aic_13),length(which(min_model == 3))/length(min)),
              c(mean(aic_14),median(aic_14),length(which(min_model == 4))/length(min)))


min_model<-c()
aic<-rbind(aic_21,aic_22,aic_23,aic_24)
min<-apply(aic,2,min)
for (i in 1:length(min)){
  min_model<-c(min_model,which(aic[,i]==min[i]))
}


DGP2<-cbind(  c(mean(aic_21),median(aic_21),length(which(min_model == 1))/length(min)),
              c(mean(aic_22),median(aic_22),length(which(min_model == 2))/length(min)),
              c(mean(aic_23),median(aic_23),length(which(min_model == 3))/length(min)),
              c(mean(aic_24),median(aic_24),length(which(min_model == 4))/length(min)))

min_model<-c()
aic<-rbind(aic_31,aic_32,aic_33,aic_34)
min<-apply(aic,2,min)
for (i in 1:length(min)){
  min_model<-c(min_model,which(aic[,i]==min[i]))
}


DGP3<-cbind(  c(mean(aic_31),median(aic_31),length(which(min_model == 1))/length(min)),
              c(mean(aic_32),median(aic_32),length(which(min_model == 2))/length(min)),
              c(mean(aic_33),median(aic_33),length(which(min_model == 3))/length(min)),
              c(mean(aic_34),median(aic_34),length(which(min_model == 4))/length(min)))

min_model<-c()
aic<-rbind(aic_41,aic_42,aic_43,aic_44)
min<-apply(aic,2,min)
for (i in 1:length(min)){
  min_model<-c(min_model,which(aic[,i]==min[i]))
}


DGP4<-cbind(  c(mean(aic_41),median(aic_41),length(which(min_model == 1))/length(min)),
              c(mean(aic_42),median(aic_42),length(which(min_model == 2))/length(min)),
              c(mean(aic_43),median(aic_43),length(which(min_model == 3))/length(min)),
              c(mean(aic_44),median(aic_44),length(which(min_model == 4))/length(min)))



AIC_sim<-rbind(DGP1,DGP2,DGP3,DGP4)

#################################################################################
# Table 2                                                                       #
#################################################################################

stargazer(AIC_sim,type="text")

