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}