This notebook shows how to perform price optimization when willingness to pay (wtp) data are available from different individuals. We show how to fit a normal distribution and a logistic distribution to the data and compute the implied demand curve. We also show how to compute the optimal profits and optimal prices, given variable costs and estimates of the market size.
Reading Data
priceData <- read.csv("/Users/gss_wacom_mm/Desktop/Asim/Fwd Files For 11:10:2017/wtpData.csv")
summary(priceData)
Wtp
Min. :39.00
1st Qu.:47.00
Median :49.00
Mean :49.52
3rd Qu.:52.00
Max. :61.00
We compute the standard deviation of the sample.
sd(priceData$Wtp)
[1] 4.443589
WTP Distributions
We make a density histogram of the WTP data below.

Normal Distribution
We can superimpose the normal density curve on the histogram to visually assess the fit of the distribution. We use the method of moments (MOM), which involves equating the sample moments (the sample mean and the sample standard deviaion) to the theoretical moments of the normal distribution.

We can also use th fitdistrplus package to fit the normal density to the data, using the method of maximum likelihood. These give the maximum likelihood estimates of the normal mean and standard deviations, along with the uncertainty estimates given by the standard errors for these parameters. We see that the ML estimates are very similar to the MOM estimates.
r
r fitn<-fitdist(priceData$Wtp, ) summary(fitn)
Fitting of the distribution ' norm ' by maximum likelihood
Parameters :
estimate Std. Error
mean 49.516667 0.5688642
sd 4.406403 0.4022476
Loglikelihood: -174.1198 AIC: 352.2397 BIC: 356.4284
Correlation matrix:
mean sd
mean 1 0
sd 0 1
Having fit a normal, we will now explore how well a logistic distribution fits the data. As above, we will first use the method of moments, and then the MLE.
Logistic Distribution
We now fit a logistic distribution to the above histogram by equating the theoretical moments (i.e, mean and variance) of the distribution with those calculated from the sample.

We now fit the logistic distribution using maximum likelihood methods.
NaNs produced
r
r summary(fitw)
Fitting of the distribution ' logis ' by maximum likelihood
Parameters :
estimate Std. Error
location 49.397316 0.5397594
scale 2.427447 0.2647926
Loglikelihood: -173.5216 AIC: 351.0431 BIC: 355.2318
Correlation matrix:
location scale
location 1.00000000 0.02085409
scale 0.02085409 1.00000000
As the normal and the logistic are both two parameter distributions, we can directly compare the likelihood. Alternatively, we can compute the BIC statistics for the two fits. We can see that the logistic results in a slightly better fit than the normal distribution, so we will use the logistic distribution for the rest of the document.
Demand Curve for Logistic WTP
We now write a function that computes the demand function based on the logistic distribution. Note that the demand at a given price p is given by (1 - cdf(p)) multiplied by the marketSize.
r
r demand<-function(x,m,s, mktSize){ mktSize*(1-plogis(x, location = m, scale = s)) }
We can plot the demand curve as follows.

Profit function
We can similarly write a function for computing the profits at any given price. Notice that we write below a function that returns the negative of the profit.
r
r profit<-function(price, cost, m, s, mktSize){ (price-cost)*demand(price, m, s, mktSize) }
We can compute the profit at two different prices as shown below.
r
r profit(35, 15, m, s, 10000)
[1] 199467.3
r
r profit(20, 15, m, s, 10000)
[1] 49999.71
Optimal Price and Profit
We can now compute the optimal profit by using the optim function. The optim command minimizes the function, so we use the negative profit, as minimizing -profit is the same as maximizing profit. We see that the optimal price is 43.70 and the optimal profit is 262560.
Bootstrapping Demand Curves
We can use bootstapping to construct 1000 bootstrap samples, compute the logistic demand parameters based on these samples, compute the profit function based on these parameter values and finally compute the profit that we obtain we set the price to our optimal price for each of the parameter values.
r
r n<-60 bootProfit<-rep(0.0, 1000) cost<-15 for(i in 1:1000){ bsample<-sample(priceData$Wtp, n, replace=TRUE) mb<-mean(bsample) sb<-sqrt(var(bsample)*3/pi^2) bootProfit[i]<-(profit(optPrice, cost, mb, sb, 10000)) }
Finally we can compute a bootstrap interval for the range of profits that we could realize, if we were to set the price to the optimal price calculated based on our original sample.
r
r quantile(bootProfit, probs=c(0.025, 0.975))
2.5% 97.5%
248058.7 274438.7
LS0tCnRpdGxlOiAiUHJpY2luZyBPcHRpbWl6YXRpb24iCm91dHB1dDogaHRtbF9ub3RlYm9vawotLS0KClRoaXMgbm90ZWJvb2sgc2hvd3MgaG93IHRvIHBlcmZvcm0gcHJpY2Ugb3B0aW1pemF0aW9uIHdoZW4gd2lsbGluZ25lc3MgdG8gcGF5ICh3dHApIGRhdGEgYXJlIGF2YWlsYWJsZSBmcm9tIGRpZmZlcmVudCBpbmRpdmlkdWFscy4gV2Ugc2hvdyBob3cgdG8gZml0IGEgbm9ybWFsIGRpc3RyaWJ1dGlvbiBhbmQgYSBsb2dpc3RpYyBkaXN0cmlidXRpb24gdG8gdGhlIGRhdGEgYW5kIGNvbXB1dGUgdGhlIGltcGxpZWQgZGVtYW5kIGN1cnZlLiBXZSBhbHNvIHNob3cgaG93IHRvIGNvbXB1dGUgdGhlIG9wdGltYWwgcHJvZml0cyBhbmQgb3B0aW1hbCBwcmljZXMsIGdpdmVuIHZhcmlhYmxlIGNvc3RzIGFuZCBlc3RpbWF0ZXMgb2YgdGhlIG1hcmtldCBzaXplLiAKCiMjIFJlYWRpbmcgRGF0YQoKYGBge3J9CnByaWNlRGF0YSA8LSByZWFkLmNzdigiL1VzZXJzL2dzc193YWNvbV9tbS9EZXNrdG9wL0FzaW0vRndkIEZpbGVzIEZvciAxMToxMDoyMDE3L3d0cERhdGEuY3N2IikKYGBgCgpgYGB7cn0Kc3VtbWFyeShwcmljZURhdGEpCgpgYGAKV2UgY29tcHV0ZSB0aGUgc3RhbmRhcmQgZGV2aWF0aW9uIG9mIHRoZSBzYW1wbGUuIAoKYGBge3J9CnNkKHByaWNlRGF0YSRXdHApCmBgYAoKIyMgV1RQIERpc3RyaWJ1dGlvbnMKCldlIG1ha2UgYSBkZW5zaXR5IGhpc3RvZ3JhbSBvZiB0aGUgV1RQIGRhdGEgYmVsb3cuCgpgYGB7ciwgZmlnLndpZHRoPTR9CgpsaWJyYXJ5KGdncGxvdDIpCgpiYXNlPC1nZ3Bsb3QocHJpY2VEYXRhLCBhZXMoeD1XdHApKSArIGdlb21faGlzdG9ncmFtKGFlcyh5PS4uZGVuc2l0eS4uKSwgYmlud2lkdGg9MywgY29sb3I9ImJsYWNrIixmaWxsPSIjRkNDRjdBIiwgc2l6ZT0wLjEpK3RoZW1lX21pbmltYWwoKSt5bGFiKCJEZW5zaXR5IikKCmJhc2UKYGBgCgojIyBOb3JtYWwgRGlzdHJpYnV0aW9uCgpXZSBjYW4gc3VwZXJpbXBvc2UgdGhlIG5vcm1hbCBkZW5zaXR5IGN1cnZlIG9uIHRoZSBoaXN0b2dyYW0gdG8gdmlzdWFsbHkgYXNzZXNzIHRoZSBmaXQgb2YgdGhlIGRpc3RyaWJ1dGlvbi4gV2UgdXNlIHRoZSBtZXRob2Qgb2YgbW9tZW50cyAoTU9NKSwgd2hpY2ggaW52b2x2ZXMgZXF1YXRpbmcgdGhlIHNhbXBsZSBtb21lbnRzICh0aGUgc2FtcGxlIG1lYW4gYW5kIHRoZSBzYW1wbGUgc3RhbmRhcmQgZGV2aWFpb24pIHRvIHRoZSB0aGVvcmV0aWNhbCBtb21lbnRzIG9mIHRoZSBub3JtYWwgZGlzdHJpYnV0aW9uLiAKCmBgYHtyLCBmaWcud2lkdGg9NH0KCm08LW1lYW4ocHJpY2VEYXRhJFd0cCkKczwtc2QocHJpY2VEYXRhJFd0cCkKCmJhc2Urc3RhdF9mdW5jdGlvbihmdW49ZG5vcm0sIGFyZ3M9bGlzdChtZWFuPW0sIHNkPXMpLCBjb2xvcj0iYmxhY2siKQpgYGAKCldlIGNhbiBhbHNvIHVzZSB0aCBmaXRkaXN0cnBsdXMgcGFja2FnZSB0byBmaXQgdGhlIG5vcm1hbCBkZW5zaXR5IHRvIHRoZSBkYXRhLCB1c2luZyB0aGUgbWV0aG9kIG9mIG1heGltdW0gbGlrZWxpaG9vZC4gVGhlc2UgZ2l2ZSB0aGUgbWF4aW11bSBsaWtlbGlob29kIGVzdGltYXRlcyBvZiB0aGUgbm9ybWFsIG1lYW4gYW5kIHN0YW5kYXJkIGRldmlhdGlvbnMsIGFsb25nIHdpdGggdGhlIHVuY2VydGFpbnR5IGVzdGltYXRlcyBnaXZlbiBieSB0aGUgc3RhbmRhcmQgZXJyb3JzIGZvciB0aGVzZSBwYXJhbWV0ZXJzLiAgV2Ugc2VlIHRoYXQgdGhlIE1MIGVzdGltYXRlcyBhcmUgdmVyeSBzaW1pbGFyIHRvIHRoZSBNT00gZXN0aW1hdGVzLgoKYGBge3J9CmxpYnJhcnkoZml0ZGlzdHJwbHVzKQpmaXRuPC1maXRkaXN0KHByaWNlRGF0YSRXdHAsICJub3JtIikKc3VtbWFyeShmaXRuKQpgYGAKCkhhdmluZyBmaXQgYSBub3JtYWwsIHdlIHdpbGwgbm93IGV4cGxvcmUgaG93IHdlbGwgYSBsb2dpc3RpYyBkaXN0cmlidXRpb24gZml0cyB0aGUgZGF0YS4gQXMgYWJvdmUsIHdlIHdpbGwgZmlyc3QgdXNlIHRoZSBtZXRob2Qgb2YgbW9tZW50cywgYW5kIHRoZW4gdGhlIE1MRS4gCgojIyBMb2dpc3RpYyBEaXN0cmlidXRpb24KCldlIG5vdyBmaXQgYSBsb2dpc3RpYyBkaXN0cmlidXRpb24gdG8gdGhlIGFib3ZlIGhpc3RvZ3JhbSBieSBlcXVhdGluZyB0aGUgdGhlb3JldGljYWwgbW9tZW50cyAoaS5lLCBtZWFuIGFuZCB2YXJpYW5jZSkgb2YgdGhlIGRpc3RyaWJ1dGlvbiB3aXRoIHRob3NlIGNhbGN1bGF0ZWQgZnJvbSB0aGUgc2FtcGxlLiAKCmBgYHtyLCBmaWcud2lkdGg9NH0KCm08LW1lYW4ocHJpY2VEYXRhJFd0cCkKczwtc3FydCh2YXIocHJpY2VEYXRhJFd0cCkqMy9waV4yKQoKYmFzZStzdGF0X2Z1bmN0aW9uKGZ1bj1kbG9naXMsIGFyZ3M9bGlzdChsb2NhdGlvbj1tLCBzY2FsZT1zKSwgY29sb3I9ImJsYWNrIikKYGBgCldlIG5vdyBmaXQgdGhlIGxvZ2lzdGljIGRpc3RyaWJ1dGlvbiB1c2luZyBtYXhpbXVtIGxpa2VsaWhvb2QgbWV0aG9kcy4KCmBgYHtyfQpsaWJyYXJ5KGZpdGRpc3RycGx1cykKZml0dzwtZml0ZGlzdChwcmljZURhdGEkV3RwLCAibG9naXMiKQpzdW1tYXJ5KGZpdHcpCmBgYApBcyB0aGUgbm9ybWFsIGFuZCB0aGUgbG9naXN0aWMgYXJlIGJvdGggdHdvIHBhcmFtZXRlciBkaXN0cmlidXRpb25zLCB3ZSBjYW4gZGlyZWN0bHkgY29tcGFyZSB0aGUgbGlrZWxpaG9vZC4gQWx0ZXJuYXRpdmVseSwgd2UgY2FuIGNvbXB1dGUgdGhlIEJJQyBzdGF0aXN0aWNzIGZvciB0aGUgdHdvIGZpdHMuIFdlIGNhbiBzZWUgdGhhdCB0aGUgbG9naXN0aWMgcmVzdWx0cyBpbiBhIHNsaWdodGx5IGJldHRlciBmaXQgdGhhbiB0aGUgbm9ybWFsIGRpc3RyaWJ1dGlvbiwgc28gd2Ugd2lsbCB1c2UgdGhlIGxvZ2lzdGljIGRpc3RyaWJ1dGlvbiBmb3IgdGhlIHJlc3Qgb2YgdGhlIGRvY3VtZW50LiAKCiMjIERlbWFuZCBDdXJ2ZSBmb3IgTG9naXN0aWMgV1RQCgpXZSBub3cgd3JpdGUgYSBmdW5jdGlvbiB0aGF0IGNvbXB1dGVzIHRoZSBkZW1hbmQgZnVuY3Rpb24gYmFzZWQgb24gdGhlIGxvZ2lzdGljIGRpc3RyaWJ1dGlvbi4gTm90ZSB0aGF0IHRoZSBkZW1hbmQgYXQgYSBnaXZlbiBwcmljZSBwIGlzIGdpdmVuIGJ5ICgxIC0gY2RmKHApKSBtdWx0aXBsaWVkIGJ5IHRoZSBtYXJrZXRTaXplLiAKCmBgYHtyfQpkZW1hbmQ8LWZ1bmN0aW9uKHgsbSxzLCBta3RTaXplKXsKICBta3RTaXplKigxLXBsb2dpcyh4LCBsb2NhdGlvbiA9IG0sIHNjYWxlID0gcykpCn0KYGBgCgpXZSBjYW4gcGxvdCB0aGUgZGVtYW5kIGN1cnZlIGFzIGZvbGxvd3MuIAoKYGBge3IsIGZpZy53aWR0aD00fQpta3RTaXplPC0xMDAwMAoKZ2dwbG90KHByaWNlRGF0YSwgYWVzKHg9V3RwKSkrc3RhdF9mdW5jdGlvbihmdW49ZGVtYW5kLCBhcmdzPWxpc3QobSxzLCBta3RTaXplKSwgY29sb3I9IkJsYWNrIikrdGhlbWVfbWluaW1hbCgpK3lsYWIoIkRlbWFuZCIpK3hsYWIoInByaWNlIikKYGBgCgoKIyMgUHJvZml0IGZ1bmN0aW9uCgpXZSBjYW4gc2ltaWxhcmx5IHdyaXRlIGEgZnVuY3Rpb24gZm9yIGNvbXB1dGluZyB0aGUgcHJvZml0cyBhdCBhbnkgZ2l2ZW4gcHJpY2UuIE5vdGljZSB0aGF0IHdlIHdyaXRlIGJlbG93IGEgZnVuY3Rpb24gdGhhdCByZXR1cm5zIHRoZSBuZWdhdGl2ZSBvZiB0aGUgcHJvZml0LiAKCmBgYHtyfQpwcm9maXQ8LWZ1bmN0aW9uKHByaWNlLCBjb3N0LCBtLCBzLCBta3RTaXplKXsKICAocHJpY2UtY29zdCkqZGVtYW5kKHByaWNlLCBtLCBzLCBta3RTaXplKQp9CmBgYAoKV2UgY2FuIGNvbXB1dGUgdGhlIHByb2ZpdCBhdCB0d28gZGlmZmVyZW50IHByaWNlcyBhcyBzaG93biBiZWxvdy4KCmBgYHtyfQpwcm9maXQoMzUsIDE1LCBtLCBzLCAxMDAwMCkKcHJvZml0KDIwLCAxNSwgbSwgcywgMTAwMDApCmBgYAoKIyMgT3B0aW1hbCBQcmljZSBhbmQgUHJvZml0CgpXZSBjYW4gbm93IGNvbXB1dGUgdGhlIG9wdGltYWwgcHJvZml0IGJ5IHVzaW5nIHRoZSBvcHRpbSBmdW5jdGlvbi4gVGhlIG9wdGltIGNvbW1hbmQgbWluaW1pemVzIHRoZSBmdW5jdGlvbiwgc28gd2UgdXNlIHRoZSBuZWdhdGl2ZSBwcm9maXQsIGFzIG1pbmltaXppbmcgLXByb2ZpdCBpcyB0aGUgc2FtZSBhcyBtYXhpbWl6aW5nIHByb2ZpdC4gV2Ugc2VlIHRoYXQgdGhlIG9wdGltYWwgcHJpY2UgaXMgNDMuNzAgYW5kIHRoZSBvcHRpbWFsIHByb2ZpdCBpcyAyNjI1NjAuIAoKYGBge3J9CnJlczwtb3B0aW0oMjIscHJvZml0LCBjb3N0PTE1LCBtPW0sIHM9cywgbWt0U2l6ZT0xMDAwMCxtZXRob2Q9IkJGR1MiLCBjb250cm9sPWxpc3QoZm5zY2FsZT0tMSkpCgpvcHRQcmljZTwtcmVzJHBhcgpvcHRQcm9maXQ8LShyZXMkdmFsdWUpCgpvcHRQcmljZQoKb3B0UHJvZml0CgpgYGAKCiMjIEJvb3RzdHJhcHBpbmcgRGVtYW5kIEN1cnZlcwoKV2UgY2FuIHVzZSBib290c3RhcHBpbmcgdG8gY29uc3RydWN0IDEwMDAgYm9vdHN0cmFwIHNhbXBsZXMsIGNvbXB1dGUgdGhlIGxvZ2lzdGljIGRlbWFuZCBwYXJhbWV0ZXJzIGJhc2VkIG9uIHRoZXNlIHNhbXBsZXMsIGNvbXB1dGUgdGhlIHByb2ZpdCBmdW5jdGlvbiBiYXNlZCBvbiB0aGVzZSBwYXJhbWV0ZXIgdmFsdWVzICBhbmQgZmluYWxseSBjb21wdXRlIHRoZSBwcm9maXQgdGhhdCB3ZSBvYnRhaW4gd2Ugc2V0IHRoZSBwcmljZSB0byBvdXIgb3B0aW1hbCBwcmljZSBmb3IgZWFjaCBvZiB0aGUgcGFyYW1ldGVyIHZhbHVlcy4gCgpgYGB7cn0KCm48LTYwCgpib290UHJvZml0PC1yZXAoMC4wLCAxMDAwKQpjb3N0PC0xNQoKZm9yKGkgaW4gMToxMDAwKXsKICBic2FtcGxlPC1zYW1wbGUocHJpY2VEYXRhJFd0cCwgbiwgcmVwbGFjZT1UUlVFKQogIG1iPC1tZWFuKGJzYW1wbGUpCiAgc2I8LXNxcnQodmFyKGJzYW1wbGUpKjMvcGleMikKICBib290UHJvZml0W2ldPC0ocHJvZml0KG9wdFByaWNlLCBjb3N0LCBtYiwgc2IsIDEwMDAwKSkKfQoKYGBgCgpGaW5hbGx5IHdlIGNhbiBjb21wdXRlIGEgYm9vdHN0cmFwIGludGVydmFsIGZvciB0aGUgcmFuZ2Ugb2YgcHJvZml0cyB0aGF0IHdlIGNvdWxkIHJlYWxpemUsIGlmIHdlIHdlcmUgdG8gc2V0IHRoZSBwcmljZSB0byB0aGUgb3B0aW1hbCBwcmljZSBjYWxjdWxhdGVkIGJhc2VkIG9uIG91ciBvcmlnaW5hbCBzYW1wbGUuIAoKYGBge3J9CnF1YW50aWxlKGJvb3RQcm9maXQsIHByb2JzPWMoMC4wMjUsIDAuOTc1KSkKYGBgCgoK