rm(list = ls())
setwd("~/GitHub/OWN_LOCAL_PATH/Section_4")
set.seed(69)
source("RICBC_choice_sets_and_incentives.R")
source("RICBC_RI_SC_Choice.R")
library(bayesm)
library(ggplot2)
library(gridExtra)
beta.vec<-c(2,-1,1) #preferences: brand, price, discount
#min and max price for uniform draws
prmin=2
prmax=4
#lambdas<-seq(0.01,4.01,by=0.025)
lambdas<-seq(0.01,4.01,by=0.025) #sequence of information processing cost
rho<-c(2,0) #complex discount levels
rho.prob <- rep(1, length(rho))/length(rho)
# Simulation --------------------------------------------------------------
N = 1000
L = length(lambdas)
sim.list1<-array(list(),L)
for (i in 1:L){
sim.list1[[L]] <-array(list(), N)
}
betas<-matrix(c(rep(0,(length(beta.vec))*L)),ncol=length(beta.vec))
standarderrors<-matrix(c(rep(0,(length(beta.vec))*L)),ncol=length(beta.vec))
for(j in 1:L){
lambda<-lambdas[j]
t1 <- Sys.time()
for(i in 1:N){
X <- matrix(c(1,runif(1,prmin,prmax),
0,0), byrow = TRUE, ncol=2)
# Calculation of the states.
states.and.prior <- CreateStatesAndPrior(X, beta.vec, rho, rho.prob)
# RI Optimization with Shannon Costs
choice.probs.output <- CalcChoiceProbsUnderRIWithShannon(
Omega = states.and.prior[[1]],
mu = states.and.prior[[2]],
lambda = lambda,
max.iter = 10^7,
precision = 10^(-10))
state.dependent.choice.probs<-choice.probs.output[[2]]
uncond.choice.probs<-choice.probs.output[[1]]
eqm.mutual.info <- MutualInfo(prior = states.and.prior[[2]],
signal.dist = choice.probs.output[[1]],
posterior = choice.probs.output[[3]])
# Draw the state of the world.
curr.state <- sample(1:ncol(states.and.prior[[1]]), 1, replace = TRUE,
prob = states.and.prior[[2]])
#curr.state <- sample(1:ncol(states.and.prior[[1]]), 1, replace = TRUE,
#prob = rep(0.25,4))
# Combine current state and conditional choice probs to draw the choice
choice <- sample(1:nrow(X), 1, replace = TRUE,
prob = choice.probs.output[[2]][, curr.state])
# Output
choice.sim <- list(X, choice, curr.state, states.and.prior,
eqm.mutual.info, rho, rho.prob, # Regulation and its distribution
beta.vec, lambda,state.dependent.choice.probs,
uncond.choice.probs)
names(choice.sim) <- c("X", "choice", "curr.state", "states.and.prior",
"eqm.mutual.info", "rho", "rho.prob", "beta.vec", "lambda",
"state.dependent.choice.probs","uncond.choice.probs")
# -----------------------------------------------------------------------
# print(paste0("simulation number ", i))
sim.list1[[j]][[i]] <- choice.sim
}
t2 <- Sys.time()
t2 - t1
############################## Estimation #################################
###########################################################################
# obfuscated price bonuses at each state
num.alternatives = 2 #dim(sim.list1[[1]]$X)[1]
reg.states = expand.grid(rep(list(rho), num.alternatives - 1 ))
reg.states = rbind(t(reg.states),0)
## Model 2: MNL with both the observable price and the separate obfuscated price
## (bonus) as attributes
X = NULL
y = rep(0,N)
for(t in 1:N){
#add the bonus to the design
adj.design = cbind(sim.list1[[j]][[t]]$X, reg.states[,sim.list1[[j]][[t]]$curr.state])
colnames(adj.design)[3] = "bonus"
#storing individual data in y and X
X = rbind(X,adj.design)
y[t] = sim.list1[[j]][[t]]$choice
}
Data = list(X = X, y = y, p = num.alternatives)
Mcmc = list(R = 10000)
out_double_price1 = rmnlIndepMetrop(Data=Data, Mcmc=Mcmc)
betadraw <- colMeans(out_double_price1$betadraw)
names(betadraw) <- c("brand", "price", "discount")
#draws<-table(y)
SE<-c(sqrt(var(out_double_price1$betadraw[,1])),
sqrt(var(out_double_price1$betadraw[,2])),
sqrt(var(out_double_price1$betadraw[,3])))
#Bonus levels: -2,0
betas[j,]<-round(betadraw,2)
standarderrors[j,]<-round(SE,2)
}
colnames(betas)<-c("brand","price","discount")
#plotting results
g<-list("Lambda"=lambdas,"Beta_b"=betas[,1],"Beta_p"=betas[,2],
"Beta_d"=betas[,3])
df<-data.frame(g)
test1<-ggplot(data=df, aes(x=Lambda,y=Beta_b))+
geom_point( size=1) +
ggtitle("Brand")+
theme(axis.text=element_text(size=8),
axis.title=element_text(size=10),
legend.title = element_blank(),
legend.position = "top",
plot.title = element_text(hjust = 0.5, size=10, face = "bold"))+
#theme(plot.title = element_text(hjust = 0.5))+
labs(x="Information Processing Cost",y="Estimated Part Worth")+
coord_cartesian( ylim = c(0, 60))
test2<-ggplot(data=df, aes(x=Lambda,y=Beta_p))+
geom_point( size=1) +
ggtitle("Price")+
theme(axis.text=element_text(size=8),
axis.title=element_text(size=10),
legend.title = element_blank(),
legend.position = "top",
plot.title = element_text(hjust = 0.5, size=10, face = "bold"))+
#theme(plot.title = element_text(hjust = 0.5))+
labs(x="Information Processing Cost",y="Estimated Part Worth")+
coord_cartesian( ylim = c(0, -30))
test3<-ggplot(data=df, aes(x=Lambda,y=Beta_d))+
geom_point( size=1) +
ggtitle("Discount")+
theme(axis.text=element_text(size=8),
axis.title=element_text(size=10),
legend.title = element_blank(),
legend.position = "top",
plot.title = element_text(hjust = 0.5, size=10, face = "bold"))+
#theme(plot.title = element_text(hjust = 0.5))+
labs(x="Information Processing Cost",y="Estimated Part Worth")+
coord_cartesian( ylim = c(0, 30))
windows(height=3, width=9)
par(mfrow=c(1,3))
grid.arrange(test1,test2,test3,ncol=3)
#######################################################################
####Calculating Mc Fadden R?###########################################
#no of alternatives
nalt<-2
#no of lambda sims
NL<-length(sim.list1)
#no of simulated choces per lambda
N<-length(sim.list1[[1]])
probs<-array(rep(0,length(sim.list1[[1]])*nalt*length(sim.list1)),
dim = c(length(sim.list1),length(sim.list1[[1]]),nalt))
ys<-matrix(rep(0,length(sim.list1)*length(sim.list1[[1]])),
ncol=length(sim.list1[[1]]))
zwischenprobs<-matrix(rep(0,length(sim.list1)*length(sim.list1[[1]])),
ncol=length(sim.list1[[1]]))
loglkmod<-rep(0,length(sim.list1))
nullchoiceprobs<-matrix(rep(0,NL*nalt),ncol=nalt)
for (j in 1:length(probs[,1,1])){
for (i in 1:length(probs[1,,1])){
#Logit choice probs
probs[j,i,1]<-exp(betas[j,]%*%c(sim.list1[[j]][[i]]$X[1,],
sim.list1[[j]][[i]]$rho[sim.list1[[j]][[i]]$curr.state]))/
(exp(betas[j,]%*%c(sim.list1[[j]][[i]]$X[1,],
sim.list1[[j]][[i]]$rho[sim.list1[[j]][[i]]$curr.state]))+1)
probs[j,i,2]<- 1-probs[j,i,1]
#choice matrix, j rows is grid of lambdas and i columns are choice simulations
#per lambda
ys[j,i]<-sim.list1[[j]][[i]]$choice
#helpmatrix
zwischenprobs[j,i]<-log(probs[j,i,ys[j,i]])
}
loglkmod[j]<-sum(zwischenprobs[j,])
nullchoiceprobs[j,]<-c(sum(ys[j,]==1)/N, 1-sum(ys[j,]==1)/N)
}
loglknull<-rep(0,NL)
zwischennull<-matrix(rep(0,N*NL),ncol=N)
Mcfaddenrsquare<-rep(0,NL)
for (j in 1:NL){
for (i in 1:N){
zwischennull[j,i]<-log(nullchoiceprobs[j,ys[j,i]])
}
loglknull[j]<-sum(zwischennull[j,])
Mcfaddenrsquare[j]<-1-(loglkmod[j]/loglknull[j])
}
g2<-list("Lambda"=lambdas,"Mcfadd"=Mcfaddenrsquare)
df2<-data.frame(g)
test4<-ggplot(data=df2, aes(x=Lambda,y=Mcfaddenrsquare))+
geom_point( size=1) +
ggtitle("Model Fit")+
theme(axis.text=element_text(size=8),
axis.title=element_text(size=10),
legend.title = element_blank(),
legend.position = "top",
plot.title = element_text(hjust = 0.5, size=10, face = "bold"))+
#theme(plot.title = element_text(hjust = 0.5))+
labs(x="Information Processing Cost",y = expression ("MC Fadden"~R^2))+
coord_cartesian( ylim = c(0.5,1))
windows(height=3, width=5.5)
par(mfrow=c(1,1))
test4