bayesm (version 3.1-6)

rbayesBLP: Bayesian Analysis of Random Coefficient Logit Models Using Aggregate Data

Description

rbayesBLP implements a hybrid MCMC algorithm for aggregate level sales data in a market with differentiated products. bayesm version 3.1-0 and prior verions contain an error when using instruments with this function; this will be fixed in a future version.

Usage

rbayesBLP(Data, Prior, Mcmc)

Value

A list containing:

thetabardraw

\(K x R/keep\) matrix of random coefficient mean draws

Sigmadraw

\(K*K x R/keep\) matrix of random coefficient variance draws

rdraw

\(K*K x R/keep\) matrix of \(r\) draws (same information as in Sigmadraw)

tausqdraw

\(R/keep x 1\) vector of aggregate demand shock variance draws

Omegadraw

\(2*2 x R/keep\) matrix of correlated endogenous shock variance draws

deltadraw

\(I x R/keep\) matrix of endogenous structural equation coefficient draws

acceptrate

scalor of acceptance rate of Metropolis-Hasting

s

scale parameter used for Metropolis-Hasting

cand_cov

var-cov matrix used for Metropolis-Hasting

Arguments

Data

list(X, share, J, Z)

Prior

list(sigmasqR, theta_hat, A, deltabar, Ad, nu0, s0_sq, VOmega)

Mcmc

list(R, keep, nprint, H, initial_theta_bar, initial_r, initial_tau_sq, initial_Omega, initial_delta, s, cand_cov, tol)

Argument Details

Data = list(X, share, J, Z) [Z optional]

J:	number of alternatives, excluding an outside option
X:	\(J*T x K\) matrix (no outside option, which is normalized to 0).
	If IV is used, the last column of X is the endogeneous variable.
share:	\(J*T\) vector (no outside option).
	Note that both the share vector and the X matrix are organized by the \(jt\) index.
	\(j\) varies faster than \(t\), i.e. \((j=1,t=1), (j=2,t=1), ..., (j=J,T=1), ..., (j=J,t=T)\)
Z:	\(J*T x I\) matrix of instrumental variables (optional)

Prior = list(sigmasqR, theta_hat, A, deltabar, Ad, nu0, s0_sq, VOmega) [optional]

sigmasqR:	\(K*(K+1)/2\) vector for \(r\) prior variance (def: diffuse prior for \(\Sigma\))
theta_hat:	\(K\) vector for \(\theta_bar\) prior mean (def: 0 vector)
A:	\(K x K\) matrix for \(\theta_bar\) prior precision (def: 0.01*diag(K))
deltabar:	\(I\) vector for \(\delta\) prior mean (def: 0 vector)
Ad:	\(I x I\) matrix for \(\delta\) prior precision (def: 0.01*diag(I))
nu0:	d.f. parameter for \(\tau_sq\) and \(\Omega\) prior (def: K+1)
s0_sq:	scale parameter for \(\tau_sq\) prior (def: 1)
VOmega:	\(2 x 2\) matrix parameter for \(\Omega\) prior (def: matrix(c(1,0.5,0.5,1),2,2))

Mcmc = list(R, keep, nprint, H, initial_theta_bar, initial_r, initial_tau_sq, initial_Omega, initial_delta, s, cand_cov, tol) [only R and H required]

R:	number of MCMC draws
keep:	MCMC thinning parameter -- keep every keepth draw (def: 1)
nprint:	print the estimated time remaining for every nprint'th draw (def: 100, set to 0 for no print)
H:	number of random draws used for Monte-Carlo integration
initial_theta_bar:	initial value of \(\theta_bar\) (def: 0 vector)
initial_r:	initial value of \(r\) (def: 0 vector)
initial_tau_sq:	initial value of \(\tau_sq\) (def: 0.1)
initial_Omega:	initial value of \(\Omega\) (def: diag(2))
initial_delta:	initial value of \(\delta\) (def: 0 vector)
s:	scale parameter of Metropolis-Hasting increment (def: automatically tuned)
cand_cov:	var-cov matrix of Metropolis-Hasting increment (def: automatically tuned)
tol:	convergence tolerance for the contraction mapping (def: 1e-6)

Model Details

Model and Priors (without IV):

\(u_ijt = X_jt \theta_i + \eta_jt + e_ijt\)
\(e_ijt\) \(\sim\) type I Extreme Value (logit)
\(\theta_i\) \(\sim\) \(N(\theta_bar, \Sigma)\)
\(\eta_jt\) \(\sim\) \(N(0, \tau_sq)\)

This structure implies a logit model for each consumer (\(\theta\)). Aggregate shares (share) are produced by integrating this consumer level logit model over the assumed normal distribution of \(\theta\).

\(r\) \(\sim\) \(N(0, diag(sigmasqR))\)
\(\theta_bar\) \(\sim\) \(N(\theta_hat, A^-1)\)
\(\tau_sq\) \(\sim\) \(nu0*s0_sq / \chi^2 (nu0)\)

Note: we observe the aggregate level market share, not individual level choices.

Note: \(r\) is the vector of nonzero elements of cholesky root of \(\Sigma\). Instead of \(\Sigma\) we draw \(r\), which is one-to-one correspondence with the positive-definite \(\Sigma\).

Model and Priors (with IV):

\(u_ijt = X_jt \theta_i + \eta_jt + e_ijt\)
\(e_ijt\) \(\sim\) type I Extreme Value (logit)
\(\theta_i\) \(\sim\) \(N(\theta_bar, \Sigma)\)

\(X_jt = [X_exo_jt, X_endo_jt]\)
\(X_endo_jt = Z_jt \delta_jt + \zeta_jt\)
\(vec(\zeta_jt, \eta_jt)\) \(\sim\) \(N(0, \Omega)\)

\(r\) \(\sim\) \(N(0, diag(sigmasqR))\)
\(\theta_bar\) \(\sim\) \(N(\theta_hat, A^-1)\)
\(\delta\) \(\sim\) \(N(deltabar, Ad^-1)\)
\(\Omega\) \(\sim\) \(IW(nu0, VOmega)\)

MCMC and Tuning Details:

MCMC Algorithm:

Step 1 (\(\Sigma\)):
Given \(\theta_bar\) and \(\tau_sq\), draw \(r\) via Metropolis-Hasting.
Covert the drawn \(r\) to \(\Sigma\).

Note: if user does not specify the Metropolis-Hasting increment parameters (s and cand_cov), rbayesBLP automatically tunes the parameters.

Step 2 without IV (\(\theta_bar\), \(\tau_sq\)):
Given \(\Sigma\), draw \(\theta_bar\) and \(\tau_sq\) via Gibbs sampler.

Step 2 with IV (\(\theta_bar\), \(\delta\), \(\Omega\)):
Given \(\Sigma\), draw \(\theta_bar\), \(\delta\), and \(\Omega\) via IV Gibbs sampler.

Tuning Metropolis-Hastings algorithm:

r_cand = r_old + s*N(0,cand_cov)
Fix the candidate covariance matrix as cand_cov0 = diag(rep(0.1, K), rep(1, K*(K-1)/2)).
Start from s0 = 2.38/sqrt(dim(r))

Repeat{
Run 500 MCMC chain.
If acceptance rate < 30% => update s1 = s0/5.
If acceptance rate > 50% => update s1 = s0*3.
(Store r draws if acceptance rate is 20~80%.)
s0 = s1
} until acceptance rate is 30~50%

Scale matrix C = s1*sqrt(cand_cov0)
Correlation matrix R = Corr(r draws)
Use C*R*C as s^2*cand_cov.

Author

Peter Rossi and K. Kim, Anderson School, UCLA, .

References

For further discussion, see Bayesian Analysis of Random Coefficient Logit Models Using Aggregate Data by Jiang, Manchanda, and Rossi, Journal of Econometrics, 2009.

Examples

if(nchar(Sys.getenv("LONG_TEST")) != 0) {

## Simulate aggregate level data
simulData <- function(para, others, Hbatch) {
  # Hbatch does the integration for computing market shares
  #      in batches of size Hbatch

  ## parameters
  theta_bar <- para$theta_bar
  Sigma <- para$Sigma
  tau_sq <- para$tau_sq

  T <- others$T	
  J <- others$J	
  p <- others$p	
  H <- others$H	
  K <- J + p	

  ## build X	
  X <- matrix(runif(T*J*p), T*J, p)
  inter <- NULL
  for (t in 1:T) { inter <- rbind(inter, diag(J)) }
  X <- cbind(inter, X)

  ## draw eta ~ N(0, tau_sq)	
  eta <- rnorm(T*J)*sqrt(tau_sq)
  X <- cbind(X, eta)

  share <- rep(0, J*T)
  for (HH in 1:(H/Hbatch)){
    ## draw theta ~ N(theta_bar, Sigma)
    cho <- chol(Sigma)
    theta <- matrix(rnorm(K*Hbatch), nrow=K, ncol=Hbatch)
    theta <- t(cho)%*%theta + theta_bar

    ## utility
    V <- X%*%rbind(theta, 1)
    expV <- exp(V)
    expSum <- matrix(colSums(matrix(expV, J, T*Hbatch)), T, Hbatch)
    expSum <- expSum %x% matrix(1, J, 1)
    choiceProb <- expV / (1 + expSum)
    share <- share +  rowSums(choiceProb) / H
  }

  ## the last K+1'th column is eta, which is unobservable.
  X <- X[,c(1:K)]	
  return (list(X=X, share=share))
}

## true parameter
theta_bar_true <- c(-2, -3, -4, -5)
Sigma_true <- rbind(c(3,2,1.5,1), c(2,4,-1,1.5), c(1.5,-1,4,-0.5), c(1,1.5,-0.5,3))
cho <- chol(Sigma_true)
r_true <- c(log(diag(cho)), cho[1,2:4], cho[2,3:4], cho[3,4]) 
tau_sq_true <- 1

## simulate data
set.seed(66)
T <- 300
J <- 3
p <- 1
K <- 4
H <- 1000000
Hbatch <- 5000

dat <- simulData(para=list(theta_bar=theta_bar_true, Sigma=Sigma_true, tau_sq=tau_sq_true),
                 others=list(T=T, J=J, p=p, H=H), Hbatch)
X <- dat$X
share <- dat$share

## Mcmc run
R <- 2000
H <- 50
Data1 <- list(X=X, share=share, J=J)
Mcmc1 <- list(R=R, H=H, nprint=0)
set.seed(66)
out <- rbayesBLP(Data=Data1, Mcmc=Mcmc1)

## acceptance rate
out$acceptrate

## summary of draws
summary(out$thetabardraw)
summary(out$Sigmadraw)
summary(out$tausqdraw)

### plotting draws
plot(out$thetabardraw)
plot(out$Sigmadraw)
plot(out$tausqdraw)
}