Before we start, let’s first define two helper functions. MatrixShift is used in LiF. It enqueues a row at the end of a matrix M while dequeuing a row from the start of the matrix. It then returns the new matrix. PlotAnimation returns the animated GIF at the end of this page (if doPlot is set to TRUE). This function makes use of R libraries “animation” and “cairo”.

doPlot <- FALSE

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

plotAnimation <- function() {

    Cairo(file = paste0("S", i/200, ".png"), bg = "white", width = 1300,height = 900)  
    par(mfrow=c(3,2))
    plot(qd, op, type="l", xlab="Quantity", ylab="Dollars", main="Demand", lwd=2)
    abline(h=0, lwd=2)
    lines(qd, mc, col="red", lwd=2)
    plot(minmax(op),minmax(tr,tp), type="n", ylab="", xlab="Price", main="Optimal pricing")
    grid()
    abline(h=0, lwd=2)
    abline(v=optprice, col="red", lwd=2)
    lines(op,tr, col="blue", lwd=2)
    lines(op,tp, col="red", lwd=2)
    plot(arr.wtp, ylim=c(35,60))
    plot(arr.P, type="l",ylim=c(35,60))
    plot(arr.doesbuy, type="l")
    plot(arr.p0, type="l", xlab="Time in stream", ylab="selected price",ylim=c(35,60))
    lines(arr.p0, col="red",ylim=c(35,60))
    dev.off() 
}
minmax <- function(...) c(min(...),max(...))

Now on to the meat of our simulation! To keep our model relatively simple, we will simulate a single price monopolist model of pricing, marginal cost, demand and revenue.

The general assumptions of a monopoly market simulation are:

• Only one firm is in the market, with a market share of 100%.
• A unique product with no close substitutes
• Entry of firms is not permitted (completely blocked).

In other words, a single price monopolist is the only supplier of a particular good. The monopolist therefore has the power to set a price to sell the product at. Those who have a willingness to pay which is greater than the price will buy the good. Those who have a willingness to pay for the good which is less than the chosen price will not buy it.

Let’s first generate our consumers:

npeep <- 2000                           # Number of potential consumers
wtp <- 45 + rnorm(npeep,0)*4            # Each person has a different willingness to pay
wtp[wtp<0] <- 0                         # if less than zero, they really don't want to buy, set to zero
wtp_sliced <- c()

To figure out the demand curve we count the number of people willing to pay at least as much as the offering price.

maxop <- 120                            # Maximum offering price
op <- 0:maxop                           # Offering Prices (ranges from 0 to maxop)
qd <- rep(NA,maxop+1)                   # Quantity Demanded

Set LiF parameters:

TT <- npeep                             # Length of stream                    
T <- 10                                 # Integration time
w <- (2*pi)/T                           # Omega
A <- 3                                  # Amlitude LIF
gamma <- 0.05                           # Learnrate
p0 <- 55                                # Startvalue of price

Some variables to help us log our sales:

arr.P <- rep(NA, TT)                    # Price log
arr.doesbuy <- rep(NA, TT)              # Does or doesn't buy log
arr.p0 <- rep(NA,TT)                    # Startprice log
arr.wtp <- rep(NA,TT)                   # Willingness To Pay log

yws <- rep(NA, T)  
doesbuy <- 0

We will also need Total Revenue, Total Profit, and Total Cost vectors:

tr <- tp <- tc <- rep(NA,maxop+1) 

Now iterate through all of our customers, where each iteration selling prices is calculated by LiF while logging our sales:

for(i in 1:npeep){
  
  # Log willingness to pay
  wtp_sliced <- c(wtp_sliced, wtp[i])
  arr.wtp[i] <- wtp[i]
  
  for (j in 1:maxop+1) qd[j] <- sum(wtp_sliced>=op[j])
  
  # Marginal cost is increasing in our model - let's see if LiF can keep up!
  mc <- qd*0.005
  
  # Calculate Quantity Demanded
  qd.gain <- qd[-length(qd)]-qd[-1]
  qd.gain[length(qd.gain)+1] <- qd.gain[length(qd.gain)]
  
  # Calculate Total Cost
  for (j in 1:length(op)) tc[j] <- sum((mc*qd.gain)[length(qd):j])
  
  # Calculate Total Revenue
  tr <- qd*op
  
  # Calculate Total Price
  tp <- tr-tc
  
  # Log indexmax and Optprice for our plots
  if (doPlot) { 
    indexmax <- which.max(tp)
    optprice <- op[indexmax]
  }

  # LIF: Oscillate around P
  P <- p0 + A*cos(w*i)
  arr.P[i] <- P 
  
  # If willingness to pay is smaller than the price, don't buy
  # If willingness to pay is equal to, or greater than the price, buy
  if (wtp[i]<P) {
    doesbuy <- 0
  } else {
    doesbuy <- 1
  }
  
  # Log doesbuy
  arr.doesbuy[i] <- doesbuy

  # Compute Profit:
  yt <- doesbuy*(P-i*0.005)
  yws <- matrixShift(yws, yt * cos(w*i))
  
  # LIF: Reset p0
  if(i > T){
    yw <- sum(yws) / T
    p0 <- p0 + (gamma / T) * yw
  }
  
  # Log p0
  arr.p0[i] <- p0
  
  # plot a snapshot of the logs every 200 iterations
  if (i %% 200 == 0 && doPlot) { 
    plotAnimation
  }
}

Now let’s see if LiF has indeed found, and kept track of, our selling price:

print(p0)

## [1] 54.6034

Seems a reasonable enough a price, but to really know how the price as determined by LiF corresponds to the optimal value we can infer from our monopolist assumptions, we need to take a closer look at our logs.

Or maybe we can do better – let’s plot the log files as animated GIFs, to get an even better feel for how LiF, and our simulated model, behaved:

if (doPlot) {
    ani.options(convert = shQuote("convert.exe"),ani.width = 1300, ani.height = 900, interval = 0.05, ani.dev = "png",ani.type = "png")
    bm.files = sprintf("S%i.png", 1:(npeep/200) )
    im.convert(bm.files, output = paste0("ani",format(Sys.time(), "%d%b%H%M%S"),".gif"), extra.opts = "-dispose Background", clean = T)  
}

Link to full-size plot.

Let’s assume my plots are numbered from 1 to 6, left to right, then top to bottom, like so:

… then the numbered plots would describe the following processes: