1\name{rbayesBLP} 2\alias{rbayesBLP} 3\concept{bayes} 4\concept{random coefficient logit} 5\concept{BLP} 6\concept{Metropolis Hasting} 7 8\title{Bayesian Analysis of Random Coefficient Logit Models Using Aggregate Data} 9 10\description{ 11\code{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. 12} 13 14\usage{rbayesBLP(Data, Prior, Mcmc)} 15 16\arguments{ 17 \item{Data }{list(X, share, J, Z)} 18 \item{Prior}{list(sigmasqR, theta_hat, A, deltabar, Ad, nu0, s0_sq, VOmega)} 19 \item{Mcmc }{list(R, keep, nprint, H, initial_theta_bar, initial_r, initial_tau_sq, initial_Omega, initial_delta, s, cand_cov, tol)} 20} 21 22\value{ 23 A list containing: 24 \item{thetabardraw }{\eqn{K x R/keep} matrix of random coefficient mean draws} 25 \item{Sigmadraw }{\eqn{K*K x R/keep} matrix of random coefficient variance draws} 26 \item{rdraw }{\eqn{K*K x R/keep} matrix of \eqn{r} draws (same information as in \code{Sigmadraw})} 27 \item{tausqdraw }{\eqn{R/keep x 1} vector of aggregate demand shock variance draws} 28 \item{Omegadraw }{\eqn{2*2 x R/keep} matrix of correlated endogenous shock variance draws} 29 \item{deltadraw }{\eqn{I x R/keep} matrix of endogenous structural equation coefficient draws} 30 \item{acceptrate }{scalor of acceptance rate of Metropolis-Hasting} 31 \item{s }{scale parameter used for Metropolis-Hasting} 32 \item{cand_cov }{var-cov matrix used for Metropolis-Hasting} 33} 34 35\section{Argument Details}{ 36 \emph{\code{Data = list(X, share, J, Z)}} [\code{Z} optional] 37 \tabular{ll}{ 38 \code{J: } \tab number of alternatives, excluding an outside option \cr 39 \code{X: } \tab \eqn{J*T x K} matrix (no outside option, which is normalized to 0). \cr 40 \code{ } \tab If IV is used, the last column of \code{X} is the endogeneous variable. \cr 41 \code{share: } \tab \eqn{J*T} vector (no outside option). \cr 42 \code{ } \tab Note that both the \code{share} vector and the \code{X} matrix are organized by the \eqn{jt} index. \cr 43 \code{ } \tab \eqn{j} varies faster than \eqn{t}, i.e. \eqn{(j=1,t=1), (j=2,t=1), ..., (j=J,T=1), ..., (j=J,t=T)} \cr 44 \code{Z: } \tab \eqn{J*T x I} matrix of instrumental variables (optional) 45 } 46 \emph{\code{Prior = list(sigmasqR, theta_hat, A, deltabar, Ad, nu0, s0_sq, VOmega)} [optional]} 47 \tabular{ll}{ 48 \code{sigmasqR: } \tab \eqn{K*(K+1)/2} vector for \eqn{r} prior variance (def: diffuse prior for \eqn{\Sigma}) \cr 49 \code{theta_hat: } \tab \eqn{K} vector for \eqn{\theta_bar} prior mean (def: 0 vector) \cr 50 \code{A: } \tab \eqn{K x K} matrix for \eqn{\theta_bar} prior precision (def: \code{0.01*diag(K)}) \cr 51 \code{deltabar: } \tab \eqn{I} vector for \eqn{\delta} prior mean (def: 0 vector) \cr 52 \code{Ad: } \tab \eqn{I x I} matrix for \eqn{\delta} prior precision (def: \code{0.01*diag(I)}) \cr 53 \code{nu0: } \tab d.f. parameter for \eqn{\tau_sq} and \eqn{\Omega} prior (def: K+1) \cr 54 \code{s0_sq: } \tab scale parameter for \eqn{\tau_sq} prior (def: 1) \cr 55 \code{VOmega: } \tab \eqn{2 x 2} matrix parameter for \eqn{\Omega} prior (def: \code{matrix(c(1,0.5,0.5,1),2,2)}) 56 } 57 \emph{\code{Mcmc = list(R, keep, nprint, H, initial_theta_bar, initial_r, initial_tau_sq, initial_Omega, initial_delta, s, cand_cov, tol)} [only \code{R} and \code{H} required]} 58 \tabular{ll}{ 59 \code{R: } \tab number of MCMC draws \cr 60 \code{keep: } \tab MCMC thinning parameter -- keep every \code{keep}th draw (def: 1) \cr 61 \code{nprint: } \tab print the estimated time remaining for every \code{nprint}'th draw (def: 100, set to 0 for no print) \cr 62 \code{H: } \tab number of random draws used for Monte-Carlo integration \cr 63 \code{initial_theta_bar: } \tab initial value of \eqn{\theta_bar} (def: 0 vector) \cr 64 \code{initial_r: } \tab initial value of \eqn{r} (def: 0 vector) \cr 65 \code{initial_tau_sq: } \tab initial value of \eqn{\tau_sq} (def: 0.1) \cr 66 \code{initial_Omega: } \tab initial value of \eqn{\Omega} (def: \code{diag(2)}) \cr 67 \code{initial_delta: } \tab initial value of \eqn{\delta} (def: 0 vector) \cr 68 \code{s: } \tab scale parameter of Metropolis-Hasting increment (def: automatically tuned) \cr 69 \code{cand_cov: } \tab var-cov matrix of Metropolis-Hasting increment (def: automatically tuned) \cr 70 \code{tol: } \tab convergence tolerance for the contraction mapping (def: 1e-6) 71 } 72} 73 74\section{Model Details}{ 75 \subsection{Model and Priors (without IV):}{ 76 \eqn{u_ijt = X_jt \theta_i + \eta_jt + e_ijt}\cr 77 \eqn{e_ijt} \eqn{\sim}{~} type I Extreme Value (logit)\cr 78 \eqn{\theta_i} \eqn{\sim}{~} \eqn{N(\theta_bar, \Sigma)}\cr 79 \eqn{\eta_jt} \eqn{\sim}{~} \eqn{N(0, \tau_sq)}\cr 80 81 This structure implies a logit model for each consumer (\eqn{\theta}). 82 Aggregate shares (\code{share}) are produced by integrating this consumer level 83 logit model over the assumed normal distribution of \eqn{\theta}. 84 85 \eqn{r} \eqn{\sim}{~} \eqn{N(0, diag(sigmasqR))}\cr 86 \eqn{\theta_bar} \eqn{\sim}{~} \eqn{N(\theta_hat, A^-1)}\cr 87 \eqn{\tau_sq} \eqn{\sim}{~} \eqn{nu0*s0_sq / \chi^2 (nu0)}\cr 88 89 Note: we observe the aggregate level market share, not individual level choices.\cr 90 91 Note: \eqn{r} is the vector of nonzero elements of cholesky root of \eqn{\Sigma}. 92 Instead of \eqn{\Sigma} we draw \eqn{r}, which is one-to-one correspondence with the positive-definite \eqn{\Sigma}. 93 } 94 \subsection{Model and Priors (with IV):}{ 95 \eqn{u_ijt = X_jt \theta_i + \eta_jt + e_ijt}\cr 96 \eqn{e_ijt} \eqn{\sim}{~} type I Extreme Value (logit)\cr 97 \eqn{\theta_i} \eqn{\sim}{~} \eqn{N(\theta_bar, \Sigma)}\cr 98 99 \eqn{X_jt = [X_exo_jt, X_endo_jt]}\cr 100 \eqn{X_endo_jt = Z_jt \delta_jt + \zeta_jt}\cr 101 \eqn{vec(\zeta_jt, \eta_jt)} \eqn{\sim}{~} \eqn{N(0, \Omega)}\cr 102 103 \eqn{r} \eqn{\sim}{~} \eqn{N(0, diag(sigmasqR))}\cr 104 \eqn{\theta_bar} \eqn{\sim}{~} \eqn{N(\theta_hat, A^-1)}\cr 105 \eqn{\delta} \eqn{\sim}{~} \eqn{N(deltabar, Ad^-1)}\cr 106 \eqn{\Omega} \eqn{\sim}{~} \eqn{IW(nu0, VOmega)}\cr 107 } 108} 109 110\section{MCMC and Tuning Details:}{ 111 \subsection{MCMC Algorithm:}{ 112 Step 1 (\eqn{\Sigma}):\cr 113 Given \eqn{\theta_bar} and \eqn{\tau_sq}, draw \eqn{r} via Metropolis-Hasting.\cr 114 Covert the drawn \eqn{r} to \eqn{\Sigma}.\cr 115 116 Note: if user does not specify the Metropolis-Hasting increment parameters 117 (\code{s} and \code{cand_cov}), \code{rbayesBLP} automatically tunes the parameters. 118 119 Step 2 without IV (\eqn{\theta_bar}, \eqn{\tau_sq}):\cr 120 Given \eqn{\Sigma}, draw \eqn{\theta_bar} and \eqn{\tau_sq} via Gibbs sampler.\cr 121 122 Step 2 with IV (\eqn{\theta_bar}, \eqn{\delta}, \eqn{\Omega}):\cr 123 Given \eqn{\Sigma}, draw \eqn{\theta_bar}, \eqn{\delta}, and \eqn{\Omega} via IV Gibbs sampler.\cr 124 } 125 126 \subsection{Tuning Metropolis-Hastings algorithm:}{ 127 r_cand = r_old + s*N(0,cand_cov)\cr 128 Fix the candidate covariance matrix as cand_cov0 = diag(rep(0.1, K), rep(1, K*(K-1)/2)).\cr 129 Start from s0 = 2.38/sqrt(dim(r))\cr 130 131 Repeat\{\cr 132 Run 500 MCMC chain.\cr 133 If acceptance rate < 30\% => update s1 = s0/5.\cr 134 If acceptance rate > 50\% => update s1 = s0*3.\cr 135 (Store r draws if acceptance rate is 20~80\%.)\cr 136 s0 = s1\cr 137 \} until acceptance rate is 30~50\% 138 139 Scale matrix C = s1*sqrt(cand_cov0)\cr 140 Correlation matrix R = Corr(r draws)\cr 141 Use C*R*C as s^2*cand_cov. 142 } 143} 144 145\references{For further discussion, see \emph{Bayesian Analysis of Random Coefficient Logit Models Using Aggregate Data} by Jiang, Manchanda, and Rossi, \emph{Journal of Econometrics}, 2009. \cr 146} 147 148\author{Keunwoo Kim, Anderson School, UCLA, \email{keunwoo.kim@gmail.com}.} 149 150\examples{ 151if(nchar(Sys.getenv("LONG_TEST")) != 0) { 152 153## Simulate aggregate level data 154simulData <- function(para, others, Hbatch) { 155 # Hbatch does the integration for computing market shares 156 # in batches of size Hbatch 157 158 ## parameters 159 theta_bar <- para$theta_bar 160 Sigma <- para$Sigma 161 tau_sq <- para$tau_sq 162 163 T <- others$T 164 J <- others$J 165 p <- others$p 166 H <- others$H 167 K <- J + p 168 169 ## build X 170 X <- matrix(runif(T*J*p), T*J, p) 171 inter <- NULL 172 for (t in 1:T) { inter <- rbind(inter, diag(J)) } 173 X <- cbind(inter, X) 174 175 ## draw eta ~ N(0, tau_sq) 176 eta <- rnorm(T*J)*sqrt(tau_sq) 177 X <- cbind(X, eta) 178 179 share <- rep(0, J*T) 180 for (HH in 1:(H/Hbatch)){ 181 ## draw theta ~ N(theta_bar, Sigma) 182 cho <- chol(Sigma) 183 theta <- matrix(rnorm(K*Hbatch), nrow=K, ncol=Hbatch) 184 theta <- t(cho)\%*\%theta + theta_bar 185 186 ## utility 187 V <- X\%*\%rbind(theta, 1) 188 expV <- exp(V) 189 expSum <- matrix(colSums(matrix(expV, J, T*Hbatch)), T, Hbatch) 190 expSum <- expSum \%x\% matrix(1, J, 1) 191 choiceProb <- expV / (1 + expSum) 192 share <- share + rowSums(choiceProb) / H 193 } 194 195 ## the last K+1'th column is eta, which is unobservable. 196 X <- X[,c(1:K)] 197 return (list(X=X, share=share)) 198} 199 200## true parameter 201theta_bar_true <- c(-2, -3, -4, -5) 202Sigma_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)) 203cho <- chol(Sigma_true) 204r_true <- c(log(diag(cho)), cho[1,2:4], cho[2,3:4], cho[3,4]) 205tau_sq_true <- 1 206 207## simulate data 208set.seed(66) 209T <- 300 210J <- 3 211p <- 1 212K <- 4 213H <- 1000000 214Hbatch <- 5000 215 216dat <- simulData(para=list(theta_bar=theta_bar_true, Sigma=Sigma_true, tau_sq=tau_sq_true), 217 others=list(T=T, J=J, p=p, H=H), Hbatch) 218X <- dat$X 219share <- dat$share 220 221## Mcmc run 222R <- 2000 223H <- 50 224Data1 <- list(X=X, share=share, J=J) 225Mcmc1 <- list(R=R, H=H, nprint=0) 226set.seed(66) 227out <- rbayesBLP(Data=Data1, Mcmc=Mcmc1) 228 229## acceptance rate 230out$acceptrate 231 232## summary of draws 233summary(out$thetabardraw) 234summary(out$Sigmadraw) 235summary(out$tausqdraw) 236 237### plotting draws 238plot(out$thetabardraw) 239plot(out$Sigmadraw) 240plot(out$tausqdraw) 241} 242} 243 244\keyword{models}