rbayesBLP {bayesm} R Documentation

## 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)`

### 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)

### 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

### 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 Σ) `theta_hat: ` K vector for θ_bar prior mean (def: 0 vector) `A: ` K x K matrix for θ_bar prior precision (def: `0.01*diag(K)`) `deltabar: ` I vector for δ prior mean (def: 0 vector) `Ad: ` I x I matrix for δ prior precision (def: `0.01*diag(I)`) `nu0: ` d.f. parameter for τ_sq and Ω prior (def: K+1) `s0_sq: ` scale parameter for τ_sq prior (def: 1) `VOmega: ` 2 x 2 matrix parameter for Ω 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 `keep`th 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 θ_bar (def: 0 vector) `initial_r: ` initial value of r (def: 0 vector) `initial_tau_sq: ` initial value of τ_sq (def: 0.1) `initial_Omega: ` initial value of Ω (def: `diag(2)`) `initial_delta: ` initial value of δ (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 θ_i + η_jt + e_ijt
e_ijt ~ type I Extreme Value (logit)
θ_i ~ N(θ_bar, Σ)
η_jt ~ N(0, τ_sq)

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

r ~ N(0, diag(sigmasqR))
θ_bar ~ N(θ_hat, A^-1)
τ_sq ~ nu0*s0_sq / χ^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 Σ. Instead of Σ we draw r, which is one-to-one correspondence with the positive-definite Σ.

#### Model and Priors (with IV):

u_ijt = X_jt θ_i + η_jt + e_ijt
e_ijt ~ type I Extreme Value (logit)
θ_i ~ N(θ_bar, Σ)

X_jt = [X_exo_jt, X_endo_jt]
X_endo_jt = Z_jt δ_jt + ζ_jt
vec(ζ_jt, η_jt) ~ N(0, Ω)

r ~ N(0, diag(sigmasqR))
θ_bar ~ N(θ_hat, A^-1)
Ω ~ IW(nu0, VOmega)

### MCMC and Tuning Details:

#### MCMC Algorithm:

Step 1 (Σ):
Given θ_bar and τ_sq, draw r via Metropolis-Hasting.
Covert the drawn r to Σ.

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

Step 2 without IV (θ_bar, τ_sq):
Given Σ, draw θ_bar and τ_sq via Gibbs sampler.

Step 2 with IV (θ_bar, δ, Ω):
Given Σ, draw θ_bar, δ, and Ω 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(s)

Keunwoo Kim, Anderson School, UCLA, keunwoo.kim@gmail.com.

### 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)