Replicating BLP (1995)

This assigment will attempt to replicate and extend the result of (Berry, Levinsohn, and Pakes 1995).

Problem 1: load and explore the data

The data from (Berry, Levinsohn, and Pakes 1995) is included in the “hdm” package.

if (!require(hdm)) install.packages("hdm")
library(hdm)
data(BLP)
## Note: despite what the hdm package's documentation says, price is not
## log price; it is price - mean(price). If you compare the quantiles
## of price with Table II from BLP (1995), you see that everything is
## shifted down by 11.761
tab2 <- c(3.393, 6.711, 8.728, 13.074, 68.597)
stopifnot(all(abs(quantile(BLP$price, c(0,0.25,0.5,0.75,1)) - tab2 +
                  11.761)<0.005))
BLP$price <- BLP$price + 11.761

Create some tables and figures to explore the data. You may want to reproduce table I and/or II from (Berry, Levinsohn, and Pakes 1995).

Problem 2: Logit Demand

Reproduce table III of (Berry, Levinsohn, and Pakes 1995). For the IV logit estimates, you can use hdm:::constructIV to create the instruments. See (Chernozhukov, Hansen, and Spindler 2016) for a description of hdm:::constructIV and an example of usage. (I could not obtain the same results as (Berry, Levinsohn, and Pakes 1995) in column 2.) Calculate the elasticities of demand implied by the logit model. Report how many own price elasticities have absolute value less than one (inelastic). Why are inelastic demand estimates undesirable?

Problem 3: shares and delta

The purpose of this problem is to write functions that calculate shares given delta and delta given observed shares. These functions will be used to calculate the objective function used in estimation. You can choose to organize your code however you want, but future code snippets will be more useful to you if use the same interface. The share of the \(j\)th car is \[ \int \frac{e^{\delta_j + \sum_k \sigma_k x_{jkt} v_{ik} + \alpha \log(y_i - p_j)}} { e^{\alpha \log(y_i)} + \sum_{\ell=1}^J e^{\delta_\ell + \sum_k \sigma_k x_{\ell kt} v_{ik} + \alpha \log(y_i - p_\ell)}} dF(y,v) \] We approximate this by taking \(S\) draws from the distribution of \(y\) and \(v\) and averaging. ((Berry, Levinsohn, and Pakes 1995) use a version of importance sampling, feel free to try their procedure instead if you want.) Each component of \(v\) is drawn from independent \(N(0,1)\). \(y\) is log-normal with mean and standard deviation set equal to the mean and standard deviation of log income in the CPS in each year. I downloaded the IPUMS-CPS and found that average log from 1971-1991 (income in thousands of 1983 $) = 3.082 with sd=0.840. There was little variation across years, but perhaps I miss-used the CPS. (There are changes in top-coding of income across years that I was not careful about). Anyway, I think it is adequate to let \(\log(y) \sim N(3.082, 0.840^2)\), but you could try to do better here if you want. Write a function to compute shares with the following interface

share.fn <- function(delta,  ## J vector
                     x,      ## J by K
                     log.y,  ## S vector of log(y_i)
                     log.yp, ## S by J of log(y_i - p_j) 
                     v,      ## S by K
                     alpha,  ## scalar
                     sigma)  ## K vector
{
  stop("## TODO: compute J vector of shares in this market")
  return(share)
}

Using this share function, you can compute delta given shares using the contraction mapping in equation 6.8 of (Berry, Levinsohn, and Pakes 1995).

delta.fn <- function(s,x,log.y,log.yp,v,alpha,sigma,
                     tol=1e-6, tol.s=1e-8,
                     max.iter=100)
{
  delta.new <- log(s) - log(1-sum(s)) ## initial guess
  dd <- 1 ## ||change in delta||
  ds <- 1 ## ||observed shares - share.fn(delta)||
  iter <- 0
  while ((dd>tol) | (ds>tol.s)) {
    delta.old <- delta.new
    delta.new <- stop("## TODO: update delta using contraction mapping")
    dd <- max(abs(delta.old - delta.new))
    ds <- max(abs(s - sm))
    iter = iter+1
    if (iter>max.iter) {
      warning(sprintf("Maximum iterations (%d) reached, returning with norm(delta.new - delta.old) = %.2g", 
                   max.iter, dd))
      break;
    }
  }
  return(delta.new)
}

Test your delta.fn by checking that delta.fn(share.fn(delta,…),..) = delta. Something like,

## Create instruments
## demand instruments
Z <- hdm:::constructIV(BLP$firm.id, BLP$cdid, BLP$id,
                       cbind(1,BLP[,c("hpwt","air","mpd","mpg","space","price")]))
## supply instruments
W <- log(BLP[,c("hpwt","mpg","space","mpd")])
colnames(W) <- paste("log",colnames(W),sep=".")
Wiv <- hdm:::constructIV(BLP$firm.id, BLP$cdid, BLP$id, W)

## Draws for simluting integral
S <- 100 
T <- length(unique(BLP$cdid))
K <- 5
set.seed(41658)
y.s <- matrix(exp(rnorm(S*T,mean=3.082, sd=0.840)),nrow=T,ncol=S)
v.s <- array(rnorm(S*T*K, mean=0,sd=1), dim=c(T,S,K))

## Estimates from Table IV of BLP -- used for comparison and testing
est.blp <- list(alpha=43.501, sigma=c(3.612, 4.628, 1.818, 1.050,
                                      2.056),
                beta=c(-7.061, 2.883, 1.521, -0.122, 3.460),
                gamma=c(0.952, 0.477, 0.619, -.415, -.049, .019))

## Put data into more convenient structure for estimation
est.data <- list(x=as.matrix(cbind(1,BLP[,c("hpwt","air","mpd","space")])),
                 w=as.matrix(cbind(1,log(BLP$hpwt), BLP$air, log(BLP$mpg),
                                   log(BLP$space), BLP$trend)))
est.data$zd <- as.matrix(cbind(est.data$x, Z))
est.data$zs <- as.matrix(cbind(est.data$w, Wiv))
est.data$log.y <- log(y.s)
est.data$log.yp <- list()

## BLP uses log(y - p) in utility function, but some vehicles have
## p>y for some people. BLP do not say what they did in this cases.
## I will take a first order Taylor expansion of log to the left of
## some small number to get an almost log function that is defined
## everywhere. This is very arbitrary though ....
x0 <- 0.1 ## take linear taylor approx to log around x0 for y-p<x0 to
logx0 <- log(x0)
slope <- 1/x0
## avoid log(-)
my.log <- function(x)  ifelse(x>=x0, log(x), logx0 + slope*(x-x0))
dmy.log <- function(x) ifelse(x>=x0, 1/x, slope)
for (t in 1:T) {
  yp <- outer(drop(y.s[t,]), BLP$price[BLP$cdid==t],
              function(x,y) x-y)
  est.data$log.yp[[t]] <- my.log(yp)
  est.data$dlog.yp[[t]] <- dmy.log(yp)
}

## Testing of delta.fn 
t <- 1
inc <- BLP$cdid==t
delta <- rnorm(n=length(BLP$price[inc]))
s <- share.fn(delta, x=drop(est.data$x[inc,]),
              log.y=drop(est.data$log.y[t,]),
              log.yp=est.data$log.yp[[t]],
              v=drop(v.s[t,,]),
              alpha=est.blp$alpha,
              sigma=est.blp$sigma)
d.check <- delta.fn(s, x=drop(est.data$x[inc,]),
                    log.y=drop(est.data$log.y[t,]),
                    log.yp=est.data$log.yp[[t]],
                    v=drop(v.s[t,,]),
                    alpha=est.blp$alpha,
                    sigma=est.blp$sigma, max.iter=1000)
summary((abs(delta-d.check)))

The supply side of the model requires calculating the derivative of shares with respect to price. Create and test a function that does so.

dshare.dp <- function(delta,x,log.y, log.yp, dlog.yp, v,alpha,sigma)
{
  stop("TODO: compute dshare/dp")
  return(dshare)
}

## Testing dshare.dp
if (!require(numDeriv)) install.packages("numDeriv")
library(numDeriv)
t <- 1
dshare.num <- jacobian(function(p) {
  yp <- outer(drop(y.s[t,]), p,
              function(x,y) x-y)
  share.fn(delta, x=drop(est.data$x[inc,]),
           log.y=drop(est.data$log.y[t,]),
           log.yp=my.log(yp),
           v=drop(v.s[t,,]),
           alpha=est.blp$alpha,
           sigma=est.blp$sigma)
}, x=BLP$price[inc])
dshare <- dshare.dp(delta, x=drop(est.data$x[inc,]),
           log.y=drop(est.data$log.y[t,]),
           log.yp=est.data$log.yp[[t]], 
           dlog.yp=est.data$dlog.yp[[t]] ,
           v=drop(v.s[t,,]),
           alpha=est.blp$alpha,
           sigma=est.blp$sigma)
summary(as.vector(dshare-dshare.num))

As an estimator gets more involved and the code for it becomes longer, it gets more and more important to test each part of the code individually. If you can think of anything else to test, do so.

Problem 4: code optimization

With the above, delta.fn, we can compute \(\delta\) for any \(\alpha\) and \(\sigma\). Given \(\delta\) the value of \(\beta\) and \(\gamma\) that minimizes the objective function is given a weighted regression. Therefore, \(\beta\) and \(\gamma\) can be “concentrated out” so that we only have to numerically minimize over \(\alpha\) and \(\sigma\). The following code computes the objective function.

moments <- function(alpha,sigma, s,p,x,log.yp, log.y,dlog.yp,
                    v,w,zd,zs, W,
                    market.id,firm.id, delta.tol=1e-10,
                    max.iter=100)
{
  ## Find delta and omega for each market
  delta <- rep(NA, length(s))
  omega <- rep(NA, length(s))
  for (t in unique(market.id)) {
    inc <- market.id==t
    delta[inc] <- delta.fn(s[inc], x[inc,],drop(log.y[t,]),
                           log.yp[[t]] ,
                           v=drop(v[t,,]),
                           alpha,sigma, tol=delta.tol,
                           tol.s=1e-6, max.iter=max.iter)
    dShare <- dshare.dp(delta[inc], x[inc,], drop(log.y[t,]),
                        log.yp[[t]], dlog.yp[[t]] , drop(v[t,,]), alpha,sigma)
    dShare <- dShare* outer(firm.id[inc],firm.id[inc],
                            function(x,y) x==y)
    b <- solve(dShare) %*% s[inc]
    omega[inc] <- log(p[inc]-b)
  }

  ## Solve for beta and gamma
  X <- as.matrix(rbind(cbind(x,0*w), cbind(0*x,w)))
  Y <- c(delta, omega)
  Z <- as.matrix(rbind(cbind(zd,0*zs), cbind(0*zd, zs)))
  B <- solve(t(X) %*% Z %*% W %*% t(Z) %*% X) %*%
    (t(X) %*% Z %*% W %*% t(Z) %*% Y)
  beta <- B[1:ncol(x)]
  gamma <- B[(1+ncol(x)):(ncol(x)+ncol(w))]

  ## Compute GMM objective
  E <- Y - X %*% B
  g <- drop(E)*Z
  G <- colMeans(g)
  obj <- nrow(g)*t(G) %*% W %*% G
  return(list(obj=obj, beta=beta, gamma=gamma))
}

Minimization will require evaluating the objective function hundreds or thousands of times, so it is important for the code to run quickly. Use the R profiler to identify which parts of your code are slowest.

Rprof("blp.prof", line.profiling=TRUE) # start the profiler
q <- moments(est.blp$alpha,est.blp$sigma,s=BLP$share,p=BLP$price,
             x=est.data$x, log.y=est.data$log.y,
             log.yp=est.data$log.yp,
             dlog.yp=est.data$dlog.yp,
             v=v.s,w=est.data$w,zd=est.data$zd,zs=est.data$zs,
             W=diag(1,nrow=(ncol(est.data$zd)+ncol(est.data$zs))),
             market.id=BLP$cdid, firm.id=BLP$firm.id)
Rprof(NULL) # stop the profiler

summaryRprof("blp.prof", lines="both") # show the results

Attempt to speed up your code. Begin by modifying the parts that take the longest. Describe what changes you made that successfully sped up your code.

Problem 5: estimation

Estimate the parameters of the model. Compare your estimates with those in Table IV of (Berry, Levinsohn, and Pakes 1995). Compute price elasticities as in Table VI.

Problem 6: inference

Calculate standard errors and/or confidence regions for your estimates and elasticities. Refer to the section on inference in the solutions to assignment 1.

(I expect that these last two questions are rather difficult. It is sufficient to complete only one of them.)

Problem 7: bayesian estimation

Estimate the model using Bayesian methods. You may use the rbayesBLP function from the bayesm package, or attempt to use Rstan to estimate the model.

(I expect that these last two questions are rather difficult. It is sufficient to complete only one of them.)

References

Berry, Steven, James Levinsohn, and Ariel Pakes. 1995. “Automobile Prices in Market Equilibrium.” Econometrica 63 (4). The Econometric Society: pp. 841–90. http://www.jstor.org/stable/2171802.

Chernozhukov, Victor, Chris Hansen, and Martin Spindler. 2016. “Hdm: High-Dimensional Metrics.” ArXiv Preprint ArXiv:1608.00354. https://arxiv.org/abs/1608.00354.