# Simulation Codes for obtaining Table 5 in the paper 
# ``About predictions in spatial autoregressive models:
#   Optimal and almost optimal strategies'' (2016) 
#
# Number of cores: simulations done with 8 cores 
#################################################################

EM_par <- function(j, path.dir, param.RData, nb.simul = 1000)
{
 # initialization of the parameters
 load(paste(path.dir, param.RData, sep=""))
 # import the EM function
 source(paste(path.dir, "function/fEMSAR.R", sep=""))

 # we fix rho=0.35 and configuration 2 with 54 out-of sample 
 n <- nrow(x) 
 rho <- 0.35
            
 # we only focus on two situations
 # sigma = 1, 3 
 # W = W1, W2

 # the different configurations
 if(j==1)
 {outsamp <- list.param$config[[4]] 
  sigma <- 1
  W.listw <- list.param$W[[1]]
  W <- spdep::listw2mat(W.listw)
 }else
  {if(j==2)
   {outsamp <- list.param$config[[4]] 
    sigma <- 1
    W.listw <- list.param$W[[2]]
    W <- spdep::listw2mat(W.listw)
   }else
    if(j==3)
    {outsamp <- list.param$config[[4]] 
     sigma <- 3
     W.listw <- list.param$W[[1]]
     W <- spdep::listw2mat(W.listw)
    }else
    if(j==4)
    {outsamp <- list.param$config[[4]] 
     sigma <- 3
     W.listw <- list.param$W[[2]]
     W <- spdep::listw2mat(W.listw)
    }else
     if(j==5)
     {outsamp <- list.param$config[[6]] 
      sigma <- 1
      W.listw <- list.param$W[[1]]
      W <- spdep::listw2mat(W.listw)
     }else
      {if(j==6)
       {outsamp <- list.param$config[[6]] 
        sigma <- 1
        W.listw <- list.param$W[[2]]
        W <- spdep::listw2mat(W.listw)
        }else
        if(j==7)
         {outsamp <- list.param$config[[6]] 
          sigma <- 3
          W.listw <- list.param$W[[1]]
          W <- spdep::listw2mat(W.listw)
        }else
        {outsamp <- list.param$config[[6]] 
         sigma <- 3
         W.listw <- list.param$W[[2]]
         W <- spdep::listw2mat(W.listw)
        }
       }
      }


 no <- length(outsamp)

 # jacobian matrix
 A <- diag(1,n) - rho*W

  # We build Ws
  Wss <- W[-outsamp,-outsamp]
  sub.COL.nb.2 <- spdep::mat2listw(Wss, style="W")

 # initialisaion
 res.beta <- matrix(0, nb.simul, 4)    
 res.beta.EM <- matrix(0, nb.simul, 4)  
 
 res.rho <- matrix(0, nb.simul, 1)    
 res.rho.EM <- matrix(0, nb.simul, 1)   
 
 res.sigma <- matrix(0, nb.simul, 1)    
 res.sigma.EM <- matrix(0, nb.simul, 1)  
   
 for (l in 1:nb.simul)
  {
   set.seed(l)
   Y0 <- x%*%beta0 + rnorm(n, 0, sigma)    # on tire une erreur différente
   y <- solve(A, Y0)                      # y = (ys yo)
   don <- as.data.frame(x[,-1])
   don <- data.frame(don, y)

##########
 # estimate a SAR model on the S sample
   sar.res <- spdep::lagsarlm(y~., data=don[-outsamp,],
    sub.COL.nb.2 , method="eigen")
   
 # coefficient beta
   res.beta[l,] <- sar.res$coefficients
 # rho
   res.rho[l] <- sar.res$rho
 # Sigma2 
   res.sigma[l] <- sar.res$s2
   
 # EM algorithm
   Ys <- y[-outsamp]
   rho.E <- 0.2
   beta.E <- c(0, 0, 0, 0)
   sigma.E <- var(Ys)
   theta <- fEMSAR(Ys, x, W, rho.E, sigma.E, beta.E, outsamp)

  # beta coefficient
   res.beta.EM[l,] <- theta$beta
 # rho
   res.rho.EM[l] <- theta$rho
 # Sigma2 
   res.sigma.EM[l] <- theta$sigma2        
  }

 return(res=list(res.beta, res.beta.EM, 
            res.rho, res.rho.EM, 
            res.sigma, res.sigma.EM))

}

# parallel computing
path.dir <- "Z:/Thibault Pro/research/statistique spatiale/prédiction Michel Christine/article spatial economic analysis/codes/"
library("snowfall")
sfInit(parallel=TRUE, cpus=8)
result.EM <- sfClusterApplyLB(1:8, EM_par, path.dir = path.dir, param.RData="parameters/param-out.RData", nb.simul = 1000)
sfStop()

save(result.EM, file=paste(path.dir, "result/res-EM.RData", sep=""))