rm(list = ls())
setwd("~/GitHub/OWN_LOCAL_PATH/Section_3/3.5")
library(bayesm)
library(Rcpp)
library(RcppArmadillo)
source('RI_DCM2_SIM_new.R')
Rcpp::sourceCpp('BA_stable.cpp')
Rcpp::sourceCpp('RI_MH_singlecomplex_cpp.cpp')
Rcpp::sourceCpp('RI_MH_singlecomplex_cpp_out.cpp')
###############################################
############### Simulation ###################
###############################################
# Simulation from RI model for discrete choice :
# - # of alternatives: 1 inside goods, 1 outside good
# - # of attributes: 1 price (simple), 1 discount (complex)
# - realization of the complex attribute is independent of other attribute levels
# - rho are vectors of attribute levels
#
# Generating data
#number of choices
N = 2000
#cost of information processing
lambda05 = 0.5 #intermediate lambda, partial learning
#preferences
betas= c(-1,1)
#attribute levels
rho_simp = c(2.25,2.75,5) #true simple
rho_comp = c( 1,1.5,3.5) #true complex
#rho_comp = c( 1,2.5,3.5) #true complex
rho_simp = c(2.25,2.75,5) #true simple
rho_comp = c( 2,2.5,6) #true complex
#possible realizations of complex attributes for both alternatives
states_correct=rbind(rho_comp,rep(0,length(rho_comp)))
states_wrong=rbind(rho_simp,rep(0,length(rho_simp)))
#################################################################################
#################### CORRECT SPECIFICATION, Price=simple ########################
#################################################################################
#simulate choices
# sim.list05 = RI_DCM2_Sim_new(betas, lambda05, rho_comp,
# simplelvls = rho_simp, N, 2, brand=FALSE)
nalt=3
states_correct = expand.grid(rep(list(rho_comp), nalt - 1 ))
states_correct = rbind(t(as.matrix(states_correct)),0)
states_wrong = expand.grid(rep(list(rho_simp), nalt - 1 ))
states_wrong = rbind(t(as.matrix(states_wrong)),0)
mu=rep(1/ncol(states_correct),ncol(states_correct))
for (i in 1:N){
Xsimp = rbind(cbind(sample(rho_simp,nalt-1,T)), 0)
outBA = BA_upd(Xsimp, betas[-2], lambda05, mu, states_correct, N, betas[2])
curr_state <- sample(1:length(mu), 1, prob = mu)
X = cbind(Xsimp, states_correct[,curr_state])
choice = sample(1:nrow(X), 1, replace = TRUE,
prob = outBA$condprobs[, curr_state])
# Output
choice.sim <- list(X, choice, curr_state, outBA$Omega,
outBA)
names(choice.sim) <- c("X", "choice", "curr_state", "Omega", "BAoutput")
# -----------------------------------------------------------------------
sim.list05[[i]] <- choice.sim
}
# create Data lists
Data05<- vector(mode = "list", N)
for(i in 1:N){
Data05[[i]]$X = sim.list05[[i]]$X
Data05[[i]]$y = sim.list05[[i]]$choice
Data05[[i]]$v = sim.list05[[i]]$curr_state
}
probcheck=rep(0,N)
for (i in 1:N){
if(sim.list05[[i]]$BAoutput$probs[1]<0.0001 |sim.list05[[i]]$BAoutput$probs[1]>0.999)
{probcheck[i]=1}
}
table(probcheck)
###############################################
############### Estimation ###################
###############################################
####################################################
########### INTERMEDIATE LAMBDA = 0.5 ##############
####################################################
#set.seed(66) #setting seed for replication
out05=RI_MH_singlecomplex_cpp(Data=Data05,R=5000,lambda=lambda05,
simp_idx=c(0),comp_idx=c(1), #cpp indexing
states=states_correct,
stepsizes=diag(0.001,2),beta_start=c(-2,2),
priomean=c(0,0),priovar=diag(100,2))
#results
windows()
par(mfrow=c(3,1))
plot(out05$betadraw[1,])
plot(out05$betadraw[2,])
plot(out05$lhdraw)
################################################################################
##################### WRONG SPEC ###############################################
################################################################################
Data_wrg_spec05 <- Data05
#row.names(states_wrong)=c("alt1","outside") #needed to find the state in wrong spec
nstates=ncol(states_wrong) #number of distinct states
comp_idx_wrg = 1 # index of the wrongly specified complex attribute
# finding the state number in the wrong specification
for(i in 1:N){
for(j in 1:nstates){
if(isTRUE(all.equal(Data_wrg_spec05[[i]]$X[,comp_idx_wrg],states_wrong[,j]))){Data_wrg_spec05[[i]]$v=j}
}
}
Data_cat=Data_wrg_spec05
for (i in 1:length(Data05)){
Data_cat[[i]]$X=matrix(rep(0,18),nrow=3,ncol=6)
if(Data05[[i]]$X[1,1]==rho_simp[1]){
Data_cat[[i]]$X[1,1]=1
}
if(Data05[[i]]$X[2,1]==rho_simp[1]){
Data_cat[[i]]$X[2,1]=1
}
if(Data05[[i]]$X[1,1]==rho_simp[2]){
Data_cat[[i]]$X[1,1]=1
}
if(Data05[[i]]$X[2,1]==rho_simp[2]){
Data_cat[[i]]$X[2,1]=1
}
if(Data05[[i]]$X[1,1]==rho_simp[3]){
Data_cat[[i]]$X[1,1]=1
}
if(Data05[[i]]$X[2,1]==rho_simp[3]){
Data_cat[[i]]$X[2,1]=1
}
#comp
if(Data05[[i]]$X[1,2]==rho_comp[2]){
Data_cat[[i]]$X[1,2]=1
}
if(Data05[[i]]$X[2,2]==rho_comp[2]){
Data_cat[[i]]$X[2,2]=1
}
if(Data05[[i]]$X[1,2]==rho_comp[3]){
Data_cat[[i]]$X[1,2]=1
}
if(Data05[[i]]$X[2,2]==rho_comp[3]){
Data_cat[[i]]$X[2,2]=1
}
}
likelihood_cat_wrong = function(betastar){
# Function that calculates the (log) likelihood of beta/lambda for the MH step
nobs = length(Data_cat)
beta = as.vector(betastar)
################
lambda = lambda05
states = cbind(c(betastar[1],betastar[1],0),c(betastar[2],betastar[1],0),c(betastar[3],betastar[1],0),
c(betastar[1],betastar[2],0),c(betastar[2],betastar[2],0),c(betastar[3],betastar[2],0),
c(betastar[1],betastar[3],0),c(betastar[2],betastar[3],0),c(betastar[3],betastar[3],0))
simp_idx=c(3,4,5)
comp_idx=c(0,1,2)
mu=rep((1/9),9)
lik = RI_lik_single_looped_cat(Data_cat, beta, lambda, mu, states, nobs, simp_idx, comp_idx,
optimal_action_weight)
# lpost = lik + dmvnorm(betastar,mean =rep(0,length(betastar)),sigma = diag(1,length(betastar)),log = T)
return(lik[[1]])
}
optimal_action_weight=matrix(rep(1/3,6000),nrow=3)
testbeta=c(-2,-4,-6,2,4,6)
testbeta=c(1,-0.2,-4,0.2,1,6)
testbeta=out_cat$betadraw[,5000]
wrong_cat=optim(testbeta,likelihood_cat_wrong,control=list(maxit=10000,fnscale=-1))
wrong_cat
states_cat=rbind(c(1,-0.5,-3,1,-0.5,-3,1,-0.5,-3),c(1,1,1,-0.5,-0.5,-0.5,-3,-3,-3),0)
states_cat=rbind(c(-2,-3,-4,-2,-3,-4,-2,-3,-4),c(-2,-2,-2,-3,-3,-3,-4,-4,-4),0)
out_cat=RI_MH_singlecomplex_cpp_cat(Data=Data_cat,R=5000,
lambda=lambda05,simp_idx=c(3,4,5),comp_idx=c(0,1,2),
states=states_cat,
stepsizes=diag(0.1,6),beta_start=c(-2,-3,-4,1,2,3),#wrong_cat$par,
priomean=rep(0,6),priovar=diag(100,6))
plot(out_cat$lhdraw[1001:5000])
###############################################
############### Estimation ###################
###############################################
####################################################
########## INTERMEDIATE LAMBDA = 0.5 ###############
####################################################
#set.seed(66) #setting seed for replication
out_wrg_spec05=RI_MH_singlecomplex_cpp(Data=Data_wrg_spec05,R=5000,
lambda=lambda05,simp_idx=c(1),comp_idx=c(0),
states=states_wrong,
stepsizes=diag(0.0001,2),beta_start=c(-2,2),
priomean=c(0,0),priovar=diag(100,2))
#results
windows()
par(mfrow=c(3,1))
plot(out_wrg_spec05$betadraw[1,])
plot(out_wrg_spec05$betadraw[2,])
plot(out_wrg_spec05$lhdraw)
likelihood_correct = function(betastar){
# Function that calculates the (log) likelihood of beta/lambda for the MH step
#nobs = length(Data_wrg_spec)
nobs = length(Data05)
beta = as.vector(betastar)
# beta[price_idx] = - exp(beta[price_idx])
################
lambda = lambda05
states = states_correct
simp_idx=c(0) #cpp indexing
comp_idx=c(1)
mu=rep((1/9),9)
lik = RI_lik_single_looped(Data05, beta, lambda, mu, states_correct, nobs, simp_idx, comp_idx,
optimal_action_weight)
return(lik[[1]])
}
optimal_action_weight=matrix(rep(1/3,6000),nrow=3)
correct=optim(c(-4,0.2),likelihood_correct,control=list(maxit=1000,fnscale=-1))
correct
likelihood_wrong = function(betastar){
# Function that calculates the (log) likelihood of beta/lambda for the MH step
#nobs = length(Data_wrg_spec)
nobs = length(Data_wrg_spec05)
beta = as.vector(betastar)
# beta[price_idx] = - exp(beta[price_idx])
################
lambda = lambda05
states = states_wrong
simp_idx=c(1) #cpp indexing
comp_idx=c(0)
mu=rep((1/3),3)
lik = RI_lik_single_looped(Data_wrg_spec05, beta, lambda, mu, states_correct, nobs, simp_idx, comp_idx,
optimal_action_weight)
return(lik[[1]])
}
optimal_action_weight=matrix(rep(0.5,4000),nrow=2)
wrong=optim(c(-3,10),likelihood_wrong,control=list(maxit=1000,fnscale=-1))
wrong
summary(out05$lhdraw[1001:5000])
summary(out_wrg_spec05$lhdraw[4801:5000])
round(c(correct$value,wrong$value),1)