library("readxl")
library(bayesvl)
setwd("~/Desktop/ClimateChange-Business")
dat <- read_excel("StrategyandStock.xlsx")
colnames(dat) 
keeps <- c("AverageStockPrice", "EnvironmentalExpendituresInvestments","LogCO2", "LogNetIncome", "EnvironmentalProducts","RenewableEnergyUse")
dat <- dat[keeps]
dat<-na.omit(dat)
data_num <- as.data.frame(apply(dat, 2, as.numeric))
sapply(data_num, class)
attach(data_num)
#stock~Income&interactions with environmental efforts - uninformative prior
model1 <- bayesvl()

model1 <- bvl_addNode(model1, "AverageStockPrice", "norm")
model1 <- bvl_addNode(model1, "LogNetIncome", "norm")
model1 <- bvl_addNode(model1, "EnvironmentalProducts", "binom")
model1 <- bvl_addNode(model1, "EnvironmentalExpendituresInvestments", "binom")
model1 <- bvl_addNode(model1, "RenewableEnergyUse", "binom")
model1 <- bvl_addNode(model1, "LogCO2", "norm")

model1 <- bvl_addNode(model1, "Income_CO2", "trans")
model1 <- bvl_addArc(model1, "LogNetIncome", "Income_CO2", "*") 
model1 <- bvl_addArc(model1, "LogCO2", "Income_CO2", "*") 

model1 <- bvl_addNode(model1, "Income_EnvironmentalProduct", "trans")
model1 <- bvl_addArc(model1, "LogNetIncome", "Income_EnvironmentalProduct", "*") 
model1 <- bvl_addArc(model1, "EnvironmentalProducts", "Income_EnvironmentalProduct", "*") 

model1 <- bvl_addNode(model1, "Income_RenewableEnergy", "trans")
model1 <- bvl_addArc(model1, "LogNetIncome", "Income_RenewableEnergy", "*") 
model1 <- bvl_addArc(model1, "RenewableEnergyUse", "Income_RenewableEnergy", "*")

model1 <- bvl_addNode(model1, "Income_EnvironmentInvestment", "trans")
model1 <- bvl_addArc(model1, "LogNetIncome", "Income_EnvironmentInvestment", "*") 
model1 <- bvl_addArc(model1, "EnvironmentalExpendituresInvestments", "Income_EnvironmentInvestment", "*")

model1 <- bvl_addArc(model1, "LogNetIncome",  "AverageStockPrice","slope")
model1 <- bvl_addArc(model1, "EnvironmentalProducts",  "AverageStockPrice","slope")
model1 <- bvl_addArc(model1, "EnvironmentalExpendituresInvestments",  "AverageStockPrice","slope")
model1 <- bvl_addArc(model1, "RenewableEnergyUse",  "AverageStockPrice","slope")
model1 <- bvl_addArc(model1, "Income_EnvironmentalProduct",  "AverageStockPrice","slope")
model1 <- bvl_addArc(model1, "Income_EnvironmentInvestment",  "AverageStockPrice","slope")
model1 <- bvl_addArc(model1, "Income_RenewableEnergy",  "AverageStockPrice","slope")
model1 <- bvl_addArc(model1, "LogCO2",  "AverageStockPrice","slope")
model1 <- bvl_addArc(model1, "Income_CO2",  "AverageStockPrice","slope")

model_string <- bvl_model2Stan(model1)
cat(model_string)

bvl_bnPlot(model1)
# detect the number of cores of your CPU
options(mc.cores = parallel::detectCores())

# fit the model using appropriate technical parameters
model1 <- bvl_modelFit(model1, data_num, warmup = 2000, iter = 5000, chains = 4)


summary(model1)
bvl_plotTrace(model1)
bvl_plotAcfs(model1,3,4,params = NULL)
bvl_plotGelmans(model1,3,4,params = NULL)
bvl_plotDensity(model1,c("b_LogNetIncome_AverageStockPrice","b_EnvironmentalProducts_AverageStockPrice",
                         "b_EnvironmentalExpendituresInvestments_AverageStockPrice","b_RenewableEnergyUse_AverageStockPrice","b_Income_EnvironmentInvestment_AverageStockPrice","b_Income_RenewableEnergy_AverageStockPrice","b_LogCO2_AverageStockPrice","b_Income_CO2_AverageStockPrice"))+theme_bw()
bvl_plotIntervals(model1,c("b_LogNetIncome_AverageStockPrice","b_EnvironmentalProducts_AverageStockPrice",
                           "b_EnvironmentalExpendituresInvestments_AverageStockPrice","b_RenewableEnergyUse_AverageStockPrice","b_Income_EnvironmentalProduct_AverageStockPrice","b_Income_EnvironmentInvestment_AverageStockPrice","b_Income_RenewableEnergy_AverageStockPrice","b_LogCO2_AverageStockPrice","b_Income_CO2_AverageStockPrice"))+theme_bw()+geom_vline(xintercept = 0,color = "red", size=1)
bvl_plotParams(model1,credMass = 0.89,params = NULL)
loo1<-bvl_stanLoo(model1)
plot(loo1)
#stock~Income&interactions with environmental factors - Informative priors
model1a <- bayesvl()

model1a <- bvl_addNode(model1a, "AverageStockPrice", "norm")
model1a <- bvl_addNode(model1a, "LogNetIncome", "norm")
model1a <- bvl_addNode(model1a, "EnvironmentalProducts", "binom")
model1a <- bvl_addNode(model1a, "EnvironmentalExpendituresInvestments", "binom")
model1a <- bvl_addNode(model1a, "RenewableEnergyUse", "binom")
model1a <- bvl_addNode(model1a, "LogCO2", "norm")

model1a <- bvl_addNode(model1a, "Income_CO2", "trans")
model1a <- bvl_addArc(model1a, "LogNetIncome", "Income_CO2", "*") 
model1a <- bvl_addArc(model1a, "LogCO2", "Income_CO2", "*") 

model1a <- bvl_addNode(model1a, "Income_EnvironmentalProduct", "trans")
model1a <- bvl_addArc(model1a, "LogNetIncome", "Income_EnvironmentalProduct", "*") 
model1a <- bvl_addArc(model1a, "EnvironmentalProducts", "Income_EnvironmentalProduct", "*") 

model1a <- bvl_addNode(model1a, "Income_RenewableEnergy", "trans")
model1a <- bvl_addArc(model1a, "LogNetIncome", "Income_RenewableEnergy", "*") 
model1a <- bvl_addArc(model1a, "RenewableEnergyUse", "Income_RenewableEnergy", "*")

model1a <- bvl_addNode(model1a, "Income_EnvironmentInvestment", "trans")
model1a <- bvl_addArc(model1a, "LogNetIncome", "Income_EnvironmentInvestment", "*") 
model1a <- bvl_addArc(model1a, "EnvironmentalExpendituresInvestments", "Income_EnvironmentInvestment", "*")

model1a <- bvl_addArc(model1a, "LogNetIncome",  "AverageStockPrice","slope")
model1a <- bvl_addArc(model1a, "EnvironmentalProducts",  "AverageStockPrice","slope", priors = c("b_EnvironmentalProducts_AverageStockPrice ~ normal (0.5,0.5)"))
model1a <- bvl_addArc(model1a, "EnvironmentalExpendituresInvestments",  "AverageStockPrice","slope", priors = c("b_EnvironmentalExpendituresInvestments_AverageStockPrice ~ normal (0.5,0.5)"))
model1a <- bvl_addArc(model1a, "RenewableEnergyUse",  "AverageStockPrice","slope", priors = c("b_RenewableEnergyUse_AverageStockPrice ~ normal (0.5,0.5)"))
model1a <- bvl_addArc(model1a, "Income_EnvironmentalProduct",  "AverageStockPrice","slope")
model1a <- bvl_addArc(model1a, "Income_EnvironmentInvestment",  "AverageStockPrice","slope")
model1a <- bvl_addArc(model1a, "Income_RenewableEnergy",  "AverageStockPrice","slope")
model1a <- bvl_addArc(model1a, "LogCO2",  "AverageStockPrice","slope", priors = c("b_LogCO2_AverageStockPrice ~ normal (-8,8)"))
model1a <- bvl_addArc(model1a, "Income_CO2",  "AverageStockPrice","slope")

model_string <- bvl_model2Stan(model1a)
cat(model_string)

# detect the number of cores of your CPU
options(mc.cores = parallel::detectCores())

# fit the model using appropriate technical parameters
model1a <- bvl_modelFit(model1a, data_num, warmup = 2000, iter = 5000, chains = 4)


summary(model1a)

#stock~income&environmental factors - uninformative priors
model2 <- bayesvl()

model2 <- bvl_addNode(model2, "AverageStockPrice", "norm")
model2 <- bvl_addNode(model2, "EnvironmentalProducts", "binom")
model2 <- bvl_addNode(model2, "EnvironmentalExpendituresInvestments", "binom")
model2 <- bvl_addNode(model2, "RenewableEnergyUse", "binom")
model2 <- bvl_addNode(model2, "LogCO2", "norm")
model2 <- bvl_addNode(model2, "LogNetIncome", "norm")

model2 <- bvl_addArc(model2, "EnvironmentalProducts",  "AverageStockPrice","slope")
model2 <- bvl_addArc(model2, "EnvironmentalExpendituresInvestments",  "AverageStockPrice","slope")
model2 <- bvl_addArc(model2, "RenewableEnergyUse",  "AverageStockPrice","slope")
model2 <- bvl_addArc(model2, "LogCO2",  "AverageStockPrice","slope")
model2 <- bvl_addArc(model2, "LogNetIncome",  "AverageStockPrice","slope")


model_string <- bvl_model2Stan(model2)
cat(model_string)

bvl_bnPlot(model2)
# detect the number of cores of your CPU
options(mc.cores = parallel::detectCores())

# fit the model using appropriate technical parameters
model2 <- bvl_modelFit(model2, data_num, warmup = 2000, iter = 5000, chains = 4)


summary(model2)
bvl_plotTrace(model2)
bvl_plotAcfs(model2,3,2,params = NULL)
bvl_plotGelmans(model2,3,2,params = NULL)
bvl_plotDensity(model2,c("b_EnvironmentalProducts_AverageStockPrice","b_EnvironmentalExpendituresInvestments_AverageStockPrice",
                         "b_RenewableEnergyUse_AverageStockPrice","b_LogCO2_AverageStockPrice","b_LogNetIncome_AverageStockPrice"))+theme_bw()
bvl_plotIntervals(model2,c("b_EnvironmentalProducts_AverageStockPrice","b_EnvironmentalExpendituresInvestments_AverageStockPrice",
                           "b_RenewableEnergyUse_AverageStockPrice","b_LogCO2_AverageStockPrice","b_LogNetIncome_AverageStockPrice"))+theme_bw()+geom_vline(xintercept = 0,color = "red", size=1)
bvl_plotParams(model2,2,3,credMass = 0.89,params = NULL)
loo2<-bvl_stanLoo(model2)
plot(loo2)

#stock~income&environmental factors - informative priors
model2a <- bayesvl()

model2a <- bvl_addNode(model2a, "AverageStockPrice", "norm")
model2a <- bvl_addNode(model2a, "EnvironmentalProducts", "binom")
model2a <- bvl_addNode(model2a, "EnvironmentalExpendituresInvestments", "binom")
model2a <- bvl_addNode(model2a, "RenewableEnergyUse", "binom")
model2a <- bvl_addNode(model2a, "LogCO2", "norm")
model2a <- bvl_addNode(model2a, "LogNetIncome", "norm")

model2a <- bvl_addArc(model2a, "EnvironmentalProducts",  "AverageStockPrice","slope", priors = c("b_EnvironmentalProducts_AverageStockPrice ~ normal (0.5,0.5)"))
model2a <- bvl_addArc(model2a, "EnvironmentalExpendituresInvestments",  "AverageStockPrice","slope", priors = c("b_EnvironmentalExpendituresInvestments_AverageStockPrice ~ normal (0.5,0.5)"))
model2a <- bvl_addArc(model2a, "RenewableEnergyUse",  "AverageStockPrice","slope", priors = c("b_RenewableEnergyUse_AverageStockPrice ~ normal (0.5,0.5)"))
model2a <- bvl_addArc(model2a, "LogCO2",  "AverageStockPrice","slope", priors = c("b_LogCO2_AverageStockPrice ~ normal (-8,8)"))
model2a <- bvl_addArc(model2a, "LogNetIncome",  "AverageStockPrice","slope")


model_string <- bvl_model2Stan(model2a)
cat(model_string)

bvl_bnPlot(model2a)
# detect the number of cores of your CPU
options(mc.cores = parallel::detectCores())

# fit the model using appropriate technical parameters
model2a <- bvl_modelFit(model2a, data_num, warmup = 2000, iter = 5000, chains = 4)
summary(model2a)

#stock~environment factors
model3 <- bayesvl()

model3 <- bvl_addNode(model3, "AverageStockPrice", "norm")
model3 <- bvl_addNode(model3, "EnvironmentalProducts", "binom")
model3 <- bvl_addNode(model3, "EnvironmentalExpendituresInvestments", "binom")
model3 <- bvl_addNode(model3, "RenewableEnergyUse", "binom")
model3 <- bvl_addNode(model3, "LogCO2", "norm")

model3 <- bvl_addArc(model3, "EnvironmentalProducts",  "AverageStockPrice","slope")
model3 <- bvl_addArc(model3, "EnvironmentalExpendituresInvestments",  "AverageStockPrice","slope")
model3 <- bvl_addArc(model3, "RenewableEnergyUse",  "AverageStockPrice","slope")
model3 <- bvl_addArc(model3, "LogCO2",  "AverageStockPrice","slope")


model_string <- bvl_model2Stan(model3)
cat(model_string)

bvl_bnPlot(model3)
# detect the number of cores of your CPU
options(mc.cores = parallel::detectCores())

# fit the model using appropriate technical parameters
model3 <- bvl_modelFit(model3, data_num, warmup = 2000, iter = 5000, chains = 4)


summary(model3)
bvl_plotTrace(model3)
bvl_plotAcfs(model3,3,2,params = NULL)
bvl_plotGelmans(model3,3,2,params = NULL)
bvl_plotDensity(model3,c("b_EnvironmentalProducts_AverageStockPrice","b_EnvironmentalExpendituresInvestments_AverageStockPrice",
                         "b_RenewableEnergyUse_AverageStockPrice","b_LogCO2_AverageStockPrice"))+theme_bw()
bvl_plotIntervals(model3,c("b_EnvironmentalProducts_AverageStockPrice","b_EnvironmentalExpendituresInvestments_AverageStockPrice",
                           "b_RenewableEnergyUse_AverageStockPrice","b_LogCO2_AverageStockPrice"),
)+theme_bw()+geom_vline(xintercept = 0,color = "red", size=1)
bvl_plotParams(model3,2,3,credMass = 0.89,params = NULL)
loo3<-bvl_stanLoo(model3)
plot(loo3)

#stock~environmental factors - informative priors
model3a <- bayesvl()

model3a <- bvl_addNode(model3a, "AverageStockPrice", "norm")
model3a <- bvl_addNode(model3a, "EnvironmentalProducts", "binom")
model3a <- bvl_addNode(model3a, "EnvironmentalExpendituresInvestments", "binom")
model3a <- bvl_addNode(model3a, "RenewableEnergyUse", "binom")
model3a <- bvl_addNode(model3a, "LogCO2", "norm")

model3a <- bvl_addArc(model3a, "EnvironmentalProducts",  "AverageStockPrice","slope", priors = c("b_EnvironmentalProducts_AverageStockPrice ~ normal (0.5,0.5)"))
model3a <- bvl_addArc(model3a, "EnvironmentalExpendituresInvestments",  "AverageStockPrice","slope", priors = c("b_EnvironmentalExpendituresInvestments_AverageStockPrice ~ normal (0.5,0.5)"))
model3a <- bvl_addArc(model3a, "RenewableEnergyUse",  "AverageStockPrice","slope", priors = c("b_RenewableEnergyUse_AverageStockPrice ~ normal (0.5,0.5)"))
model3a <- bvl_addArc(model3a, "LogCO2",  "AverageStockPrice","slope", priors = c("b_LogCO2_AverageStockPrice ~ normal (-8,8)"))


model_string <- bvl_model2Stan(model3a)
cat(model_string)


# detect the number of cores of your CPU
options(mc.cores = parallel::detectCores())

# fit the model using appropriate technical parameters
model3a <- bvl_modelFit(model3a, data_num, warmup = 2000, iter = 5000, chains = 4)
summary(model3a)

#stock~Income
model4 <- bayesvl()

model4 <- bvl_addNode(model4, "LogNetIncome", "norm")
model4 <- bvl_addNode(model4, "AverageStockPrice", "norm")

model4 <- bvl_addArc(model4, "LogNetIncome",  "AverageStockPrice","slope")

model_string <- bvl_model2Stan(model4)
cat(model_string)

bvl_bnPlot(model4)
# detect the number of cores of your CPU
options(mc.cores = parallel::detectCores())

# fit the model using appropriate technical parameters
model4 <- bvl_modelFit(model4, data_num, warmup = 2000, iter = 5000, chains = 4)

summary(model4)
bvl_plotTrace(model4)
bvl_plotAcfs(model4)
bvl_plotGelmans(model4,params = NULL)
bvl_plotTrace(model4)
bvl_plotAcfs(model4)
bvl_plotGelmans(model4,params = NULL)
bvl_plotDensity(model4,c("b_LogNetIncome_AverageStockPrice")
)+theme_bw()
bvl_plotIntervals(model4,c("b_LogNetIncome_AverageStockPrice"),
)+theme_bw()+geom_vline(xintercept = 0,color = "red", size=1)
bvl_plotParams(model4,credMass = 0.89,params = NULL)
loo1<-bvl_stanLoo(model4)
plot(loo4)


# Weight 
library("loo")
loo1 <- bvl_stanLoo(model1)
loo2 <- bvl_stanLoo(model2)
loo3 <- bvl_stanLoo(model3)
loo4 <- bvl_stanLoo(model4)
print(loo1)
print(loo2)
print(loo3)
print(loo4)

log_lik_1 <- extract_log_lik(model1@stanfit, parameter_name = "log_lik_AverageStockPrice", merge_chains = FALSE)
r_eff1 <- relative_eff(exp(log_lik_1))
loo_1 <- loo(log_lik_1, r_eff = r_eff1, cores = 2)

log_lik_2 <- extract_log_lik(model2@stanfit, parameter_name = "log_lik_AverageStockPrice", merge_chains = FALSE)
r_eff2 <- relative_eff(exp(log_lik_2))
loo_2 <- loo(log_lik_2, r_eff = r_eff2, cores = 2)

log_lik_3 <- extract_log_lik(model3@stanfit, parameter_name = "log_lik_AverageStockPrice", merge_chains = FALSE)
r_eff3 <- relative_eff(exp(log_lik_3))
loo_3 <- loo(log_lik_3, r_eff = r_eff3, cores = 2)

log_lik_4 <- extract_log_lik(model4@stanfit, parameter_name = "log_lik_AverageStockPrice", merge_chains = FALSE)
r_eff4 <- relative_eff(exp(log_lik_4))
loo_4 <- loo(log_lik_4, r_eff = r_eff4, cores = 2)

loo_list <- list(model1=loo_1, model2=loo_2, model3=loo_3, model4=loo_4)

stacking_wts <-loo_model_weights(loo_list) # stacking weight
pbma_BB_wts<-loo_model_weights(loo_list, method = "pseudobma") # pseudo-BMA+ weight
pbma_wts <-loo_model_weights(loo_list, method = "pseudobma", BB = FALSE) # pseudo-BMA weight

(waic1 <- waic(log_lik_1))
(waic2 <- waic(log_lik_2))
(waic3 <- waic(log_lik_3))
(waic4 <- waic(log_lik_4))
waics <- c(
  waic1$estimates["elpd_waic", 1],
  waic2$estimates["elpd_waic", 1],
  waic3$estimates["elpd_waic", 1],
  waic4$estimates["elpd_waic", 1])

waic_wts <- exp(waics) / sum(exp(waics))
round(cbind(waic_wts, pbma_wts, pbma_BB_wts, stacking_wts),4 )