This is project was carried out with Vimel Meng, Hao Yu, Zijian Zhang, and Yuejing Lyu. The project involves using sales data, demograhic data and a pricing model to optimally position and price our product in a market with competiting products. This project also involved identifying customer segments through the purchase habits of certain demographics, to allow for efficient product positioning.

library("dummies")
## dummies-1.5.6 provided by Decision Patterns
library("AER")
## Loading required package: car
## Loading required package: carData
## Loading required package: lmtest
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## Loading required package: sandwich
## Loading required package: survival
library("mlogit")
## Loading required package: Formula
library("gmnl")
## Loading required package: maxLik
## Loading required package: miscTools
## 
## Please cite the 'maxLik' package as:
## Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1.
## 
## If you have questions, suggestions, or comments regarding the 'maxLik' package, please use a forum or 'tracker' at maxLik's R-Forge site:
## https://r-forge.r-project.org/projects/maxlik/
library("data.table")
library("sqldf")
## Loading required package: gsubfn
## Loading required package: proto
## Warning in doTryCatch(return(expr), name, parentenv, handler): unable to load shared object '/Library/Frameworks/R.framework/Resources/modules//R_X11.so':
##   dlopen(/Library/Frameworks/R.framework/Resources/modules//R_X11.so, 6): Library not loaded: /opt/X11/lib/libSM.6.dylib
##   Referenced from: /Library/Frameworks/R.framework/Versions/3.6/Resources/modules/R_X11.so
##   Reason: image not found
## Could not load tcltk.  Will use slower R code instead.
## Loading required package: RSQLite
library("dplyr")
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:data.table':
## 
##     between, first, last
## The following object is masked from 'package:car':
## 
##     recode
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

We begin by importing and cleaning the data

data=fread("kiwi_bubbles_P2.csv",stringsAsFactors = F)
demo = fread("demo_P2.csv",stringsAsFactors = F)

data=data[!(data$price.KB==99),]
data=data[!(data$price.KR==99),]
data=data[!(data$price.MB==99),]

Defining variables

id=data$id
purchaseKB=1*(data$choice=="KB") #Define an indicator of purchase
priceKB=data$price.KB
priceKR=data$price.KR
priceMB=data$price.MB

Converting data for MLE, running MLE model and estimating model

mlogitdata=mlogit.data(data,id="id",varying=4:7,choice="choice",shape="wide")
mle= gmnl(choice ~price, data = mlogitdata)
summary(mle)
## 
## Model estimated on: Sun Jul 26 15:14:46 2020 
## 
## Call:
## gmnl(formula = choice ~ price, data = mlogitdata, method = "nr")
## 
## Frequencies of categories:
## 
##       0      KB      KR      MB 
## 0.41564 0.18035 0.20039 0.20362 
## 
## The estimation took: 0h:0m:0s 
## 
## Coefficients:
##                Estimate Std. Error z-value  Pr(>|z|)    
## KB:(intercept)  4.25316    0.32821  12.959 < 2.2e-16 ***
## KR:(intercept)  4.36240    0.32945  13.241 < 2.2e-16 ***
## MB:(intercept)  4.20440    0.31331  13.419 < 2.2e-16 ***
## price          -3.73793    0.23671 -15.791 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Optimization of log-likelihood by Newton-Raphson maximisation
## Log Likelihood: -1909
## Number of observations: 1547
## Number of iterations: 4
## Exit of MLE: gradient close to zero
coef=mle$coefficients

Calculate own elasticity, evaluated at the average prices observed in the data - ProductKB

avgpKB = mean(priceKB)
avgpKR = mean(priceKR)
avgpMB = mean(priceMB)

probKB=exp(coef[1]+coef[4]*avgpKB)/(1+exp(coef[1]+coef[4]*avgpKB)+exp(coef[2]+coef[4]*avgpKR)+ exp(coef[3]+coef[4]*avgpMB))

elasKB=-coef[4]*avgpKB*(1-probKB)

Calculate own elasticity, evaluated at the average prices observed in the data - ProductKR

probKR=exp(coef[1]+coef[4]*avgpKR)/(1+exp(coef[1]+coef[4]*avgpKB)+exp(coef[2]+coef[4]*avgpKR)+ exp(coef[3]+coef[4]*avgpMB))
elasKR=-coef[4]*avgpKR*(1-probKR)  # 4.236821

Calculate cross-price elasticities between KB(our product) and MB (competitor’s product)

probMB=exp(coef[1]+coef[4]*avgpMB)/(1+exp(coef[1]+coef[4]*avgpKB)+exp(coef[2]+coef[4]*avgpKR)+exp(coef[3]+coef[4]*avgpMB))

elasMB_KB=-coef[4]*avgpMB*probMB

Calculate optimal prices for KB and KR (note that you have two products) to maximize the profit when MB price is 1.43.

pMB=1.43
#unit cost
uc=.5

demand=function(priceKB,priceKR,priceMB,para){
    probKB=exp(para[1]+para[4]*priceKB)/(1+exp(para[1]+para[4]*priceKB)+exp(para[2]+para[4]*priceKR)+
                                             exp(para[3]+para[4]*priceMB))
    probKR=exp(para[2]+para[4]*priceKR)/(1+exp(para[1]+para[4]*priceKB)+exp(para[2]+para[4]*priceKR)+
                                             exp(para[3]+para[4]*priceMB))
    probMB=exp(para[3]+para[4]*priceKR)/(1+exp(para[1]+para[4]*priceKB)+exp(para[2]+para[4]*priceKR)+
                                             exp(para[3]+para[4]*priceMB))
    #  probno=1/(1+exp(para[1]+para[4]*priceKB)+exp(para[2]+para[4]*priceKR)+
    #             exp(para[3]+para[4]*priceMB))
    return(cbind(probKB,probKR,probMB))
}


#Write profit as a function of prices we set and model parameters
profit=function(priceKB,priceKR,priceMB,para){
    profitKB=demand(priceKB,priceKR,priceMB,para)[,1]*(priceKB-uc)*1000
    profitKR=demand(priceKB,priceKR,priceMB,para)[,2]*(priceKR-uc)*1000
    profitMB=demand(priceKB,priceKR,priceMB,para)[,3]*(priceMB-uc)*1000
    return(cbind(profitKB,profitKR,profitMB))
}

# probablity for productKB, KR and MB
probdf=as.data.frame(demand(priceKB,priceKB,pMB,coef)) 

# price for Product KB, KR and MB
profdf=as.data.frame(profit(priceKB,priceKB,pMB,coef))

# add price to profdf
profdf$priceKB=priceKB
profdf$priceKR=priceKR
profdf$priceMB=priceMB

#sum the profit for our product 
profdf$own=profdf$profitKB+profdf$profitKR

#determine optimal price for our products based on profit
maxprof=max(profdf$own)
opKB=unique(profdf$priceKB[which(profdf$own==maxprof)]) #1.21
opKR=unique(profdf$priceKR[which(profdf$own==maxprof)]) #1.43

Determining customer segments, determining elasticity per segment and optimal price per segment per product (positioning)

data2=merge(data,demo,by="id")
mlogitdata=mlogit.data(data2,id="id",varying=4:7,choice="choice",shape="wide")


#best BIC. The important variables that affect choice. We could use these demo variables in segmentation
mle = mlogit(choice~price|
                 log(fam_size)+log(hh_edu)+log(fem_work)+log(male_educ)+male_occup+male_smoke,
             data=mlogitdata)
summary(mle)
## 
## Call:
## mlogit(formula = choice ~ price | log(fam_size) + log(hh_edu) + 
##     log(fem_work) + log(male_educ) + male_occup + male_smoke, 
##     data = mlogitdata, method = "nr")
## 
## Frequencies of alternatives:
##       0      KB      KR      MB 
## 0.41795 0.16087 0.21261 0.20857 
## 
## nr method
## 5 iterations, 0h:0m:0s 
## g'(-H)^-1g = 3.44E-05 
## successive function values within tolerance limits 
## 
## Coefficients :
##                    Estimate Std. Error  z-value  Pr(>|z|)    
## KB:(intercept)     6.638860   0.914984   7.2557 3.995e-13 ***
## KR:(intercept)     3.147239   0.887566   3.5459 0.0003912 ***
## MB:(intercept)     4.311252   0.874641   4.9292 8.258e-07 ***
## price             -3.933884   0.280913 -14.0039 < 2.2e-16 ***
## KB:log(fam_size)  -0.097136   0.212599  -0.4569 0.6477449    
## KR:log(fam_size)   0.029994   0.194842   0.1539 0.8776579    
## MB:log(fam_size)  -0.099016   0.192488  -0.5144 0.6069734    
## KB:log(hh_edu)    -1.278020   0.324754  -3.9353 8.308e-05 ***
## KR:log(hh_edu)    -0.085562   0.314238  -0.2723 0.7854048    
## MB:log(hh_edu)     1.039705   0.330445   3.1464 0.0016531 ** 
## KB:log(fem_work)   0.063910   0.235500   0.2714 0.7860993    
## KR:log(fem_work)   1.239377   0.241725   5.1272 2.941e-07 ***
## MB:log(fem_work)   0.693538   0.229459   3.0225 0.0025070 ** 
## KB:log(male_educ) -0.071736   0.319063  -0.2248 0.8221094    
## KR:log(male_educ)  0.539310   0.309023   1.7452 0.0809480 .  
## MB:log(male_educ) -1.303521   0.286296  -4.5531 5.287e-06 ***
## KB:male_occup     -0.028542   0.026085  -1.0942 0.2738657    
## KR:male_occup     -0.094846   0.023696  -4.0025 6.267e-05 ***
## MB:male_occup      0.011278   0.022749   0.4958 0.6200547    
## KB:male_smoke      0.776350   0.217609   3.5676 0.0003602 ***
## KR:male_smoke     -0.804728   0.272122  -2.9572 0.0031041 ** 
## MB:male_smoke     -0.293862   0.214105  -1.3725 0.1699045    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log-Likelihood: -1443.3
## McFadden R^2:  0.11248 
## Likelihood ratio test : chisq = 365.85 (p.value = < 2.22e-16)
BIC=log(length(data2$id))*length(mle$coefficients)-2*mle$logLik[1]
BIC
## [1] 3043.285
#3043.285



###let's carry out a cluster analysis to give us an ideal of the number of clusters in our demographic data

require("cluster")
## Loading required package: cluster
require("fpc")
## Loading required package: fpc
require("factoextra")
## Loading required package: factoextra
## Loading required package: ggplot2
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
require("gridExtra")
## Loading required package: gridExtra
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
library(cluster)
library(fpc)
library(factoextra)
library(gridExtra)

clustTest = function(toClust,print=TRUE,scale=TRUE,maxClusts=15,seed=12345,nstart=20,iter.max=100){
    if(scale){ toClust = scale(toClust);}
    set.seed(seed);   # set random number seed before doing cluster analysis
    wss <- (nrow(toClust)-1)*sum(apply(toClust,2,var))
    for (i in 2:maxClusts) wss[i] <- sum(kmeans(toClust,centers=i,nstart=nstart,iter.max=iter.max)$withinss)
    
    gpw = fviz_nbclust(toClust,kmeans,method="wss",iter.max=iter.max,nstart=nstart,k.max=maxClusts) #alternative way to get wss elbow chart.
    pm1 = pamk(toClust,scaling=TRUE)
    
    gps = fviz_nbclust(toClust,kmeans,method="silhouette",iter.max=iter.max,nstart=nstart,k.max=maxClusts) 
    if(print){
        grid.arrange(gpw,gps, nrow = 1)
    }
    list(wss=wss,pm1=pm1$nc,gpw=gpw,gps=gps)
}

runClusts = function(toClust,nClusts,print=TRUE,maxClusts=15,seed=12345,nstart=20,iter.max=100){
    if(length(nClusts)>4){
        warning("Using only first 4 elements of nClusts.")
    }
    kms=list(); ps=list();
    for(i in 1:4){
        kms[[i]] = kmeans(toClust,nClusts[i],iter.max = iter.max, nstart=nstart)
        ps[[i]] = fviz_cluster(kms[[i]], geom = "point", data = toClust) + ggtitle(paste("k =",nClusts[i]))
        
    }
    library(gridExtra)
    if(print){
        tmp = marrangeGrob(ps, nrow = 2,ncol=2)
        print(tmp)
    }
    list(kms=kms,ps=ps)
}


plotClust = function(km,toClust,discPlot=FALSE){
    nc = length(km$size)
    if(discPlot){par(mfrow=c(2,2))}
    else {par(mfrow=c(3,1))}
    percsize = paste(1:nc," = ",format(km$size/sum(km$size)*100,digits=2),"%",sep="")
    pie(km$size,labels=percsize,col=1:nc)
    
    clusplot(toClust, km$cluster, color=TRUE, shade=TRUE,
             labels=2, lines=0,col.clus=1:nc); #plot clusters against principal components
    
    if(discPlot){
        plotcluster(toClust, km$cluster,col=km$cluster); #plot against discriminant functions ()
    }
    rng = range(km$centers)
    dist = rng[2]-rng[1]
    locs = km$centers+.05*dist*ifelse(km$centers>0,1,-1)
    bm = barplot(km$centers,beside=TRUE,col=1:nc,main="Cluster Means",ylim=rng+dist*c(-.1,.1))
    text(bm,locs,formatC(km$centers,format="f",digits=1))
}

demo_sub = demo[, c(1,2,4,9,12,13,15)]
#calculate optinal number of clusters, with plot
checks = clustTest(demo_sub,print=TRUE,scale=TRUE,maxClusts=15,seed=12345,nstart=20,iter.max=100) #3 clusters

clusts = runClusts(demo_sub,c(2,3,4,5),print=TRUE,maxClusts=15,seed=12345,nstart=20,iter.max=100)

x = sqldf("select distinct(id) from data")#select unique ids 
count(x)#count the number of people 
N = 329 

#4 clusters is best using some demographic variables
demo_cluster = kmeans(x=demo_sub[,2:7], centers = 6, nstart = 1000)

cluster_id = data.frame(id =demo$id)
cluster_id$cluster = demo_cluster$cluster
data3 = merge(data,cluster_id, by = 'id', all.x = T)

data3$cluster[is.na(data3$cluster)] = 4

seg.share = c(table(demo_cluster$cluster),N - sum(table(demo_cluster$cluster))) / N

seg.share
##          1          2          3          4          5          6 
## 0.10638298 0.18541033 0.14893617 0.17021277 0.07902736 0.17021277 
##            
## 0.13981763
#0.3161094 0.2340426 0.3100304 0.1398176 

coef.est = data.frame(segment = 1:4, intercept.KB = NA, intercept.KR = NA, 
                      intercept.MB = NA, price.coef = NA) 
#Write a for-loop
for (seg in 1:4) {
    # During each loop, pick subset of data of consumers from each segment.
    data.sub = subset(data3, cluster == seg)
    
    #Using that data, the rest remains the same.
    mlogitdata=mlogit.data(data.sub,id="id",varying=4:7,choice="choice",shape="wide")
    
    #Run MLE.
    mle= gmnl(choice ~  price, data = mlogitdata)
    mle
    #Store the outcome in the coef.est matrix.
    coef.est[seg, 2:5] = mle$coefficients
}

coef.est
###own-price elasticities
##KB elasticity
demand=function(priceKB,priceKR,priceMB,para){
    prob=exp(para[1]+para[4]*priceKB)/(1+exp(para[1]+para[4]*priceKB)+exp(para[2]+para[4]*priceKR)+exp(para[3]+para[4]*priceMB))
    return(prob)
}
agg_choice=function(priceKB,priceKR,priceMB) {
    
    agg_choice=seg.share[1]*demand(priceKB,priceKR,priceMB,as.numeric(coef.est[1,2:5]))+
        seg.share[2]*demand(priceKB,priceKR,priceMB,as.numeric(coef.est[2,2:5]))+
        seg.share[3]*demand(priceKB,priceKR,priceMB,as.numeric(coef.est[3,2:5]))+
        seg.share[4]*demand(priceKB,priceKR,priceMB,as.numeric(coef.est[4,2:5]))
    
    
    return(agg_choice)
}



KB_elas = -(avgpKB/agg_choice(avgpKB,avgpKR,avgpMB))*
    (seg.share[[1]]*coef.est[1,5]*demand(avgpKB,avgpKR,avgpMB,as.numeric(coef.est[1,2:5]))*(1-demand(avgpKB,avgpKR,avgpMB,as.numeric(coef.est[1,2:5])))+
         seg.share[[2]]*coef.est[2,5]*demand(avgpKB,avgpKR,avgpMB,as.numeric(coef.est[2,2:5]))*(1-demand(avgpKB,avgpKR,avgpMB,as.numeric(coef.est[2,2:5])))+
         seg.share[[3]]*coef.est[3,5]*demand(avgpKB,avgpKR,avgpMB,as.numeric(coef.est[3,2:5]))*(1-demand(avgpKB,avgpKR,avgpMB,as.numeric(coef.est[3,2:5])))+
         seg.share[[4]]*coef.est[4,5]*demand(avgpKB,avgpKR,avgpMB,as.numeric(coef.est[4,2:5]))*(1-demand(avgpKB,avgpKR,avgpMB,as.numeric(coef.est[4,2:5]))))

KB_elas
##        1 
## 3.776154
##-4.33573 

##KR_elas
demand=function(priceKB,priceKR,priceMB,para){
    prob=exp(para[1]+para[4]*priceKR)/(1+exp(para[1]+para[4]*priceKB)+exp(para[2]+para[4]*priceKR)+exp(para[3]+para[4]*priceMB))
    return(prob)
    
}


KR_elas = -(avgpKR/agg_choice(avgpKB,avgpKR,avgpMB))*
    (seg.share[[1]]*coef.est[1,5]*demand(avgpKB,avgpKR,avgpMB,as.numeric(coef.est[1,2:5]))*(1-demand(avgpKB,avgpKR,avgpMB,as.numeric(coef.est[1,2:5])))+
         seg.share[[2]]*coef.est[2,5]*demand(avgpKB,avgpKR,avgpMB,as.numeric(coef.est[2,2:5]))*(1-demand(avgpKB,avgpKR,avgpMB,as.numeric(coef.est[2,2:5])))+
         seg.share[[3]]*coef.est[3,5]*demand(avgpKB,avgpKR,avgpMB,as.numeric(coef.est[3,2:5]))*(1-demand(avgpKB,avgpKR,avgpMB,as.numeric(coef.est[3,2:5])))+
         seg.share[[4]]*coef.est[4,5]*demand(avgpKB,avgpKR,avgpMB,as.numeric(coef.est[4,2:5]))*(1-demand(avgpKB,avgpKR,avgpMB,as.numeric(coef.est[4,2:5]))))


KR_elas
##        1 
## 3.758245
##-4.314258 
pricespace=seq(0.5,1.8,0.01)

##kB choice prob per seg
demand=function(priceKB,priceKR,priceMB,para){
    prob=exp(para[1]+para[4]*priceKB)/(1+exp(para[1]+para[4]*priceKB)+exp(para[2]+para[4]*priceKR)+exp(para[3]+para[4]*priceMB))
    return(prob)
}


plot(pricespace,demand(pricespace,avgpKR,avgpMB,as.numeric(coef.est[1,2:5])),type='l',xlab='Prices',
     ylab='Probability of purchase for KB',col="blue",lwd=20*seg.share[1],
     cex=2,cex.lab=1.5, cex.axis=1.5, cex.main=1.5, cex.sub=1.5)
lines(pricespace,demand(pricespace,avgpKR,avgpMB,as.numeric(coef.est[3,2:5])),col="red",lwd=20*seg.share[3])
lines(pricespace,demand(pricespace,avgpKR,avgpMB,as.numeric(coef.est[4,2:5])),col="orange",lwd=20*seg.share[4])
lines(pricespace,demand(pricespace,avgpKR,avgpMB,as.numeric(coef.est[2,2:5])),col="green",lwd=20*seg.share[2])
legend(1.3, 0.9, legend=c("Segment 1", "Segment 2", "Segment 3", "Segment 4"),
       col=c( "blue","green", "red", "orange"), lty=1, cex=1.1)

##KR choice prob per seg
demand=function(priceKB,priceKR,priceMB,para){
    prob=exp(para[2]+para[4]*priceKR)/(1+exp(para[1]+para[4]*priceKB)+exp(para[2]+para[4]*priceKR)+exp(para[3]+para[4]*priceMB))
    return(prob)
}


plot(pricespace,demand(avgpKB,pricespace,avgpMB,as.numeric(coef.est[1,2:5])),type='l',xlab='Prices',
     ylab='Probability of purchase for KR',col="blue",lwd=20*seg.share[1],
     cex=2,cex.lab=1.5, cex.axis=1.5, cex.main=1.5, cex.sub=1.5)
lines(pricespace,demand(avgpKB,pricespace,avgpMB,as.numeric(coef.est[3,2:5])),col="red",lwd=20*seg.share[3])
lines(pricespace,demand(avgpKB,pricespace,avgpMB,as.numeric(coef.est[4,2:5])),col="orange",lwd=20*seg.share[4])
lines(pricespace,demand(avgpKB,pricespace,avgpMB,as.numeric(coef.est[2,2:5])),col="green",lwd=20*seg.share[2])
legend(1.3, 0.7, legend=c("Segment 1", "Segment 2", "Segment 3", "Segment 4"),
       col=c( "blue","green", "red", "orange"), lty=1, cex=1.1)

##MB choice prob per seg
demand=function(priceKB,priceKR,priceMB,para){
    prob=exp(para[3]+para[4]*priceMB)/(1+exp(para[1]+para[4]*priceKB)+exp(para[2]+para[4]*priceKR)+exp(para[3]+para[4]*priceMB))
    return(prob)
}


plot(pricespace,demand(avgpKB,avgpKR,pricespace,as.numeric(coef.est[1,2:5])),type='l',xlab='Prices',
     ylab='Probability of purchase for MB',col="blue",lwd=20*seg.share[1],
     cex=2,cex.lab=1.5, cex.axis=1.5, cex.main=1.5, cex.sub=1.5)
lines(pricespace,demand(avgpKB,avgpKR,pricespace,as.numeric(coef.est[3,2:5])),col="red",lwd=20*seg.share[3])
lines(pricespace,demand(avgpKB,avgpKR,pricespace,as.numeric(coef.est[4,2:5])),col="orange",lwd=20*seg.share[4])
lines(pricespace,demand(avgpKB,avgpKR,pricespace,as.numeric(coef.est[2,2:5])),col="green",lwd=20*seg.share[2])
legend(1.3, 0.9, legend=c("Segment 1", "Segment 2", "Segment 3", "Segment 4"),
       col=c( "blue","green", "red", "orange"), lty=1, cex=1.1)

##segment 1 choice prob
demandKB=function(priceKB,priceKR,priceMB,para){
    prob=exp(para[1]+para[4]*priceKB)/(1+exp(para[1]+para[4]*priceKB)+exp(para[2]+para[4]*priceKR)+exp(para[3]+para[4]*priceMB))
    return(prob)
}

demandKR=function(priceKB,priceKR,priceMB,para){
    prob=exp(para[2]+para[4]*priceKR)/(1+exp(para[1]+para[4]*priceKB)+exp(para[2]+para[4]*priceKR)+exp(para[3]+para[4]*priceMB))
    return(prob)
}

demandMB=function(priceKB,priceKR,priceMB,para){
    prob=exp(para[3]+para[4]*priceMB)/(1+exp(para[1]+para[4]*priceKB)+exp(para[2]+para[4]*priceKR)+exp(para[3]+para[4]*priceMB))
    return(prob)
}

plot(pricespace,demandKB(pricespace,avgpKR,avgpMB,as.numeric(coef.est[1,2:5])),type='l',xlab='Prices',
     ylab='Probability of purchase for Segment 1',col="blue",lwd=20*seg.share[1],
     cex=2,cex.lab=1.5, cex.axis=1.5, cex.main=1.5, cex.sub=1.5)
lines(pricespace,demandKR(avgpKB,pricespace,avgpMB,as.numeric(coef.est[1,2:5])),col="red",lwd=20*seg.share[1])
lines(pricespace,demandMB(avgpKB,avgpKR,pricespace,as.numeric(coef.est[1,2:5])),col="orange",lwd=20*seg.share[1])
legend(1.3, 0.7, legend=c("KB", "KR", "MB"),
       col=c( "blue","red", "orange"), lty=1, cex=1.1)

#kiwi regular least price sensitive

##segment 2 choice prob
plot(pricespace,demandKB(pricespace,avgpKR,avgpMB,as.numeric(coef.est[2,2:5])),type='l',xlab='Prices',
     ylab='Probability of purchase for Segment 2',col="blue",lwd=20*seg.share[2],
     cex=2,cex.lab=1.5, cex.axis=1.5, cex.main=1.5, cex.sub=1.5)
lines(pricespace,demandKR(avgpKB,pricespace,avgpMB,as.numeric(coef.est[2,2:5])),col="red",lwd=20*seg.share[2])
lines(pricespace,demandMB(avgpKB,avgpKR,pricespace,as.numeric(coef.est[2,2:5])),col="orange",lwd=20*seg.share[2])
legend(1.3, 0.7, legend=c("KB", "KR", "MB"),
       col=c( "blue","red", "orange"), lty=1, cex=1.1)

#Mango bubbles least price sensitive

##segment 3 choice prob
plot(pricespace,demandKB(pricespace,avgpKR,avgpMB,as.numeric(coef.est[3,2:5])),type='l',xlab='Prices',
     ylab='Probability of purchase for Segment 3',col="blue",lwd=20*seg.share[3],
     cex=2,cex.lab=1.5, cex.axis=1.5, cex.main=1.5, cex.sub=1.5)
lines(pricespace,demandKR(avgpKB,pricespace,avgpMB,as.numeric(coef.est[3,2:5])),col="red",lwd=20*seg.share[3])
lines(pricespace,demandMB(avgpKB,avgpKR,pricespace,as.numeric(coef.est[3,2:5])),col="orange",lwd=20*seg.share[3])
legend(1.3, 0.7, legend=c("KB", "KR", "MB"),
       col=c( "blue","red", "orange"), lty=1, cex=1.1)

#kiwi regular least price sensitive

##segment 4 choice prob
plot(pricespace,demandKB(pricespace,avgpKR,avgpMB,as.numeric(coef.est[4,2:5])),type='l',xlab='Prices',
     ylab='Probability of purchase for Segment 4',col="blue",lwd=20*seg.share[4],
     cex=2,cex.lab=1.5, cex.axis=1.5, cex.main=1.5, cex.sub=1.5)
lines(pricespace,demandKR(avgpKB,pricespace,avgpMB,as.numeric(coef.est[4,2:5])),col="red",lwd=20*seg.share[4])
lines(pricespace,demandMB(avgpKB,avgpKR,pricespace,as.numeric(coef.est[4,2:5])),col="orange",lwd=20*seg.share[4])
legend(1.3, 0.7, legend=c("KB", "KR", "MB"),
       col=c( "blue","red", "orange"), lty=1, cex=1.1)

#kiwi bubbles less price sensitive


#KB Vs KR
plot(coef.est[1,2]-coef.est[1,3],coef.est[1,5],cex=20*seg.share[1],xlim=c(-3,3),ylim=c(-9,-1.5),
     col = "blue",pch=16,cex.lab=1.5, cex.axis=1.5, cex.main=1.5, cex.sub=1.5,
     xlab="beta_0^KB-beta_0^KR",ylab=("beta_1"))
points(coef.est[2,2]-coef.est[2,3],coef.est[2,5],cex=20*seg.share[2],col = "green",pch=16)
points(coef.est[3,2]-coef.est[3,3],coef.est[3,5],cex=20*seg.share[3],col = "chocolate",pch=16)
points(coef.est[4,2]-coef.est[4,3],coef.est[4,5],cex=20*seg.share[4],col = "red",pch=16)
legend(1.3, 0.001, legend=c("Segment 1", "Segment 2", "Segment 3", "Segment 4"),
       col=c( "blue","green", "chocolate", "red"), lty=1, cex=1.1)

#segment 4 like kiwibubbles best, followed by segnent 2

#KB Vs MB
plot(coef.est[1,2]-coef.est[1,4],coef.est[1,5],cex=20*seg.share[1],xlim=c(-3,3),ylim=c(-9,-1.5),
     col = "blue",pch=16,cex.lab=1.5, cex.axis=1.5, cex.main=1.5, cex.sub=1.5,
     xlab="beta_0^KB-beta_0^MB",ylab=("beta_1"))
points(coef.est[2,2]-coef.est[2,4],coef.est[2,5],cex=20*seg.share[2],col = "green",pch=16)
points(coef.est[3,2]-coef.est[3,4],coef.est[3,5],cex=20*seg.share[3],col = "chocolate",pch=16)
points(coef.est[4,2]-coef.est[4,4],coef.est[4,5],cex=20*seg.share[4],col = "red",pch=16)
legend(1.3, 0.001, legend=c("Segment 1", "Segment 2", "Segment 3", "Segment 4"),
       col=c( "blue","green", "chocolate", "red"), lty=1, cex=1.1)

#segment 4 like kiwibubbles best, followed by segment 1

#KR Vs MB
plot(coef.est[1,3]-coef.est[1,4],coef.est[1,5],cex=20*seg.share[1],xlim=c(-3,3),ylim=c(-9,-1.5),
     col = "blue",pch=16,cex.lab=1.5, cex.axis=1.5, cex.main=1.5, cex.sub=1.5,
     xlab="beta_0^KR-beta_0^MB",ylab=("beta_1"))
points(coef.est[2,3]-coef.est[2,4],coef.est[2,5],cex=20*seg.share[2],col = "green",pch=16)
points(coef.est[3,3]-coef.est[3,4],coef.est[3,5],cex=20*seg.share[3],col = "chocolate",pch=16)
points(coef.est[4,3]-coef.est[4,4],coef.est[4,5],cex=20*seg.share[4],col = "red",pch=16)
legend(1.3, 0.001, legend=c("Segment 1", "Segment 2", "Segment 3", "Segment 4"),
       col=c( "blue","green", "chocolate", "red"), lty=1, cex=1.1)