# Decoy Effect LiF simulation in R

We ran the R simulation below to test our hypothesis that LiF should be able to optimize decoy effects (I will add comments further clarifying my train of thought as reflected by the script within the next few days). As the simulation seems to pan out quite well, our next step will be to test our hypothesis in a more realistic setting: online, with real live participants!

For those not familiar with the “decoy effect“: in the following short clip from the TV-series Numb3rs the concept is explained pretty succinctly!

The “decoy effect” is a longstanding marketing principle that was first demonstrated in 1982 by Joel Huber and others at Duke University and then expanded upon by Stanford GSB Professor Itamar Simonson. The premise is that a product is perceived as more valuable to a buyer if he or she can compare it to another less-desirable model on the shelf or a webpage.

``````library(ggplot2)

options( warn =-1)

############################## configuration ##############################

# set working directory

setwd("C:/Users/robin/Desktop/Monopolist")

# do spline? if false, we'll do line

is_spline <- TRUE

# if do_spline, which method?
#spline_method <- "monoH.FC"  # straight lines and bends through points
#spline_method <- "periodic"  # bends, but not necessarly through points, flattens toward sine
spline_method <- "fmm"        # default, alwyas bends, through points

# if not do_spline, which method?
#nspline_method <- "linear"   #lines
nspline_method <- "constant"  #blocks

# Number of customers %% 200

customers <- 20000

# LiF settings

amount_features_start <- 20
TT <- customers
T <- 8
w <- (2*pi)/T
A <- 20
gamma <- 0.9
window.length <- 1

############################## helper functions  #########################

push <- function(A, value){
if(any(is.na(A))){
i <- sum(!is.na(A))
A[i+1] <- value
} else {
A <- c(A,value)
A <- A[-1]
}
return(A)
}

slideFunct <- function(data, window, step){
total <- length(data)
spots <- seq(from=1, to=(total-window), by=step)
result <- vector(length = length(spots))
for(i in 1:length(spots)){
result[i] <- mean(data[spots[i]:(spots[i]+window)])
}
return(result)
}

############################## sketch some graphs  #####################

x_decoy_features <- c(30,60)  # A has 30 features B has 60 features
y_decoy_price <- c(30,60)     # A costs 30 B costs 60

plot(x_decoy_features, y_decoy_price, xlim = c(0, 90), ylim = c(90,0), col="red")
text(x_decoy_features+0.2, y_decoy_price-4, labels = c("A","B"), col="red")
abline(h=60,lty = 2)
text(25, 65, labels = "Decoy values of C")
text(25, 71, labels = "along dashed line")``````

``````x1   <- c(  0, 30, 60, 90   )
f1   <- c( .4, .4, .2, .2   )

x2   <- c(   0, 30, 60, 90  )
f2   <- c(   .4, .4, .7, .0 )

x3   <- c(   0, 30, 60, 90  )
f3   <- c(  .2, .2, .1, .8  )

if (is_spline) {

s1 <- splinefun(x1, f1,method=spline_method)
plot(x1, f1, ylim = c(0, 1))
curve(s1(x), add = TRUE, col = 2, n = 1001)

s2 <- splinefun(x2, f2,method=spline_method)
plot(x2, f2, ylim = c(0, 1))
curve(s2(x), add = TRUE, col = 2, n = 1001)

s3 <- function(x){
return(  1 - s1(x) - s2(x)  )
}
plot(x3, f3, ylim = c(0, 1))
curve(s3(x), add = TRUE, col = 2, n = 1001)

} else {

s1 <- approxfun(x1, f1,method = nspline_method)
plot(x1, f1, ylim = c(0, 1))
curve(s1(x), add = TRUE, col = 2, n = 1001)

s2 <- approxfun(x2, f2,method = nspline_method)
plot(x2, f2, ylim = c(0, 1))
curve(s2(x), add = TRUE, col = 2, n = 1001)

s3 <- function(x){
return(  1 - s1(x) - s2(x)  )
}
plot(x3, f3, ylim = c(0, 1))
curve(s3(x), add = TRUE, col = 2, n = 1001)

} ``````

``````# Plot functions, see if conform sketch

p <- ggplot(data.frame(x=c(0,90)), aes(x))
p <- p + scale_y_continuous(limits = c(0, 1))
p <- p + stat_function(fun=function(x){s1(x)}, geom="line", colour="red")
p <- p + stat_function(fun=function(x){s1(x)+s2(x)}, geom="line", colour="blue")
p <- p + stat_function(fun=function(x){s1(x)+s2(x)+s3(x)}, geom="line", colour="green")
print (p)``````

``````# So, maximize sales of *B* by changing how many features *C* (both at same price)

p <- ggplot(data.frame(x=c(0,90)), aes(x))
p <- p + scale_y_continuous(limits = c(0, 1))
p <- p + stat_function(fun=function(x){s2(x)}, geom="line", colour="blue")
print (p)``````

``````############################## run a test of probs  #####################

# set amount of features to 60

feat = 60

choices_sampling <- sample( c("A","B","C"), customers, replace=TRUE, prob=c(s1(feat), s2(feat), s3(feat)))

# plot barchart of prob

choicetable = prop.table(table(choices_sampling))

barplot(choicetable, col = heat.colors(3), ylim=c(0,1))``````

``````################################ Do lif! ################################

# now see if LIF can find the optimum C
# which optimizes B, where B is most profitable

yws <- rep(NA, T)

ytm <- .3

length_features <- 90

profit_per_sale <- rep(NA,length_features)

arr.features <- rep(NA, TT)
arr.amount_features_start <- rep(NA,TT)
arr.sampleChoices <- rep(NA, TT)

############################## Lif sim loop ################################

for(i in 1:customers){

# Oscillate around features
features <- amount_features_start + A*cos(w*i)

# This is only necessary for non-spline constant
if(features > 90) features <- 90
if(features < 0) features <- 0

sampleChoice <- sample( c("A","B","C"), 1, replace=TRUE, prob=c(s1(features), s2(features), s3(features)))

if (sampleChoice == "A") {
yt <- 30*0.05     # 5% profit on \$30 of A
} else if (sampleChoice == "B"){
yt <- 60*0.06     # 6% profit of \$60 for B
} else {
yt <- 60*0.04     # 4% profit on \$60 for C
}

# rolling average of yt
ytm <- ytm + ((yt - ytm) / window.length)

yws <- push(yws, ytm * cos(w*i))

# reset amount_features_start using lif
if(i > T){
yw <- sum(yws) / T
amount_features_start <- amount_features_start + (gamma / T) * yw
}

# Log me!
profit_per_sale[i] <- yt
arr.sampleChoices[i] <- sampleChoice
arr.features[i] <- features
arr.amount_features_start[i] <- amount_features_start

}

############################## Model profit ###################################

profit_plot <- rep(NA,length_features)
for (j in 1:length_features) profit_plot[j] <- s1(j)*30*0.05 + s2(j)*60*0.06  +s3(j)*60*0.04

############################## Result #########################################

max_profit_model = which(profit_plot==max(profit_plot))

print(paste0("Maximum profit at feature X of profit function: ", max_profit_model))``````
``## [1] "Maximum profit at feature X of profit function: 63"``
``print(paste0("Maximum profit at feature X found by LiF ", amount_features_start))``
``## [1] "Maximum profit at feature X found by LiF 62.7493552775659"``
``````############################## do some plotting  ##############################

# Do some simple plots

plot(arr.features, type="l",ylim=c(0,90))``````

``````plot(profit_plot, type="l")
lines(c(max_profit_model[1],max_profit_model[1]),c(-5,max(profit_plot)), col="red")``````

``plot(slideFunct(profit_per_sale,100,50), type="l")``

``````plot(arr.amount_features_start, type="l", xlab="Time in stream", ylab="features",ylim=c(0,90))
lines(arr.amount_features_start, col="red",ylim=c(0,90))``````