rm(list = ls())

setwd("[enter directory in which results will be stored]")

library(dplyr)
library(ggplot2)

# Random seed
set.seed(30)

# Functions used below

# Function for Dawid-Sebastiani score
dss <- function(m, v, y){
  -mean(dnorm(y, mean = m, sd = sqrt(v), log = TRUE))
}
# Function for logarithmic score
logs <- function(m1, v1, nu1, m2, v2, nu2, y, design = 1){
  if (design == 1){
    d1 <- dnorm(y, mean = m1, sd = sqrt(v1))
    d2 <- dnorm(y, mean = m2, sd = sqrt(v2))  
  } else {
    aux1 <- v1*((nu1-2)/nu1)
    aux2 <- v2*((nu2-2)/nu2)
    d1 <- dt((y-m1)/sqrt(aux1), df = nu1)/sqrt(aux1)
    d2 <- dt((y-m2)/sqrt(aux2), df = nu2)/sqrt(aux2)
  }
  -mean(log(w*d1+(1-w)*d2))
}

# Identifier for simulation design
design <- 1 # 1 is uneqal MSEs, 2 is scaled t distribution
if (design == 1){
  svnm <- "unequal_MSEs_"
} else {
  svnm <- "t_distributed_forecasts_"
}

# Setup parameters
# Sample size in each iteration
n <- 1e4 
# Number of iterations
n_rep <- 1e3
# Grid of weights considered
w_grid <- seq(from = 0, to = 1, by = .005)

# Parameters for simulation 
if (design == 1){
  s2_x2 <- 1.5 
} else {
  s2_x2 <- 1
}
nu1 <- nu2 <- 5
v1 <- s2_x2 + 1
v2 <- 2
w_opt <- 1-s2_x2/(1+s2_x2)

# Initialize matrices to store results
res_cor <- res_dss_lp <- res_dss_clp <- res_dss_vulp <- 
  res_logs_lp <- res_logs_clp <- res_logs_vulp <- res_v_lp <- 
  res_v_clp <- res_v_vulp <- 
  matrix(NA, n_rep, length(w_grid))

# Loop across simulation iterations
for (jj in 1:n_rep){
  
  if (jj %% 1e3 == 0){
    print(paste("Now running iteration", jj))
  }
  
  # Draw data
  if (design == 1){
    dat <- matrix(rnorm(3*n), n, 3)
    x1 <- dat[,1]
    x2 <- sqrt(s2_x2) * dat[,2]  
  } else {
    dat <- cbind(matrix(rt(2*n, df = 5), n, 2), 
                 rnorm(n))
    x1 <- sqrt(3/5)*dat[,1]
    x2 <- sqrt(3/5)*dat[,2]
  }
  u <- dat[,3]
  y <- x1+x2+u
  
  # Loop across weight choices
  for (w in w_grid){
   
    ind <- which(w_grid == w)
    
    # Mean, average variance, disagreement D, squared error S
    m <- w*x1 + (1-w)*x2
    av <- w*v1 + (1-w)*v2
    d <- w*(1-w)*(x1-x2)^2
    s <- (y-m)^2

    # Correlation between D and S
    if (!w %in% c(0,1)){
      res_cor[jj, ind] <- cor(d,s)  
    }

    # Variance
    res_v_lp[jj, ind] <- av + mean(d)
    res_v_clp[jj, ind] <- av
    res_v_vulp[jj, ind] <- av - mean(d)
    
    # Scores
    res_dss_lp[jj, ind] <- dss(m, av + d, y) 
    res_dss_clp[jj, ind] <- dss(m, av, y)
    res_dss_vulp[jj, ind] <- dss(m, av - mean(d), y)    
    res_logs_lp[jj, ind] <- logs(x1, v1, nu1, x2, v2, nu2, y, 
                                 design = design)
    res_logs_clp[jj, ind] <- logs(m, v1, nu1, m, v2, nu2, y, 
                                  design = design)
    res_logs_vulp[jj, ind] <- logs(m, v1-mean(d), nu1, m, v2-mean(d), nu2, y, 
                                   design = design)
    
  }
}

# Collect results, prepare figures
df1 <- data.frame(w = w_grid, measure = "dss", pool = "LP", 
                  value = colMeans(res_dss_lp))
df2 <- data.frame(w = w_grid, measure = "dss", pool = "CLP", 
                  value = colMeans(res_dss_clp))
df3 <- data.frame(w = w_grid, measure = "dss", pool = "VULP", 
                  value = colMeans(res_dss_vulp))
df4 <- data.frame(w = w_grid, measure = "logs", pool = "LP", 
                  value = colMeans(res_logs_lp))
df5 <- data.frame(w = w_grid, measure = "logs", pool = "CLP", 
                  value = colMeans(res_logs_clp))
df6 <- data.frame(w = w_grid, measure = "logs", pool = "VULP", 
                  value = colMeans(res_logs_vulp))
df7 <- data.frame(w = w_grid, measure = "v", pool = "LP", 
                  value = colMeans(res_v_lp))
df8 <- data.frame(w = w_grid, measure = "v", pool = "CLP", 
                  value = colMeans(res_v_clp))
df9 <- data.frame(w = w_grid, measure = "v", pool = "VULP", 
                  value = colMeans(res_v_vulp))
df <- rbind(df1, df2, df3, df4, df5, df6, df7, df8, df9)
df_cor <- data.frame(x = w_grid, y = colMeans(res_cor))

# Save aggregate results 
mc_results <- list(df = df, df_cor = df_cor)
save(mc_results, 
     file = paste0("MC_results_", substr(svnm, 1, nchar(svnm) - 1), ".RData"))