1% Generated by roxygen2: do not edit by hand
2% Please edit documentation in R/MCMCoprobitChange.R
3\name{MCMCoprobitChange}
4\alias{MCMCoprobitChange}
5\title{Markov Chain Monte Carlo for Ordered Probit Changepoint Regression Model}
6\usage{
7MCMCoprobitChange(
8  formula,
9  data = parent.frame(),
10  m = 1,
11  burnin = 1000,
12  mcmc = 1000,
13  thin = 1,
14  tune = NA,
15  verbose = 0,
16  seed = NA,
17  beta.start = NA,
18  gamma.start = NA,
19  P.start = NA,
20  b0 = NULL,
21  B0 = NULL,
22  a = NULL,
23  b = NULL,
24  marginal.likelihood = c("none", "Chib95"),
25  gamma.fixed = 0,
26  ...
27)
28}
29\arguments{
30\item{formula}{Model formula.}
31
32\item{data}{Data frame.}
33
34\item{m}{The number of changepoints.}
35
36\item{burnin}{The number of burn-in iterations for the sampler.}
37
38\item{mcmc}{The number of MCMC iterations after burnin.}
39
40\item{thin}{The thinning interval used in the simulation.  The number of
41MCMC iterations must be divisible by this value.}
42
43\item{tune}{The tuning parameter for the Metropolis-Hastings step. Default
44of NA corresponds to a choice of 0.05 divided by the number of categories in
45the response variable.}
46
47\item{verbose}{A switch which determines whether or not the progress of the
48sampler is printed to the screen.  If \code{verbose} is greater than 0 the
49iteration number, the \eqn{\beta} vector, and the error variance are
50printed to the screen every \code{verbose}th iteration.}
51
52\item{seed}{The seed for the random number generator.  If NA, the Mersenne
53Twister generator is used with default seed 12345; if an integer is passed
54it is used to seed the Mersenne twister.  The user can also pass a list of
55length two to use the L'Ecuyer random number generator, which is suitable
56for parallel computation.  The first element of the list is the L'Ecuyer
57seed, which is a vector of length six or NA (if NA a default seed of
58\code{rep(12345,6)} is used).  The second element of list is a positive
59substream number. See the MCMCpack specification for more details.}
60
61\item{beta.start}{The starting values for the \eqn{\beta} vector.
62This can either be a scalar or a column vector with dimension equal to the
63number of betas.  The default value of of NA will use the MLE estimate of
64\eqn{\beta} as the starting value.  If this is a scalar, that value
65will serve as the starting value mean for all of the betas.}
66
67\item{gamma.start}{The starting values for the \eqn{\gamma} vector.
68This can either be a scalar or a column vector with dimension equal to the
69number of gammas.  The default value of of NA will use the MLE estimate of
70\eqn{\gamma} as the starting value.  If this is a scalar, that value
71will serve as the starting value mean for all of the gammas.}
72
73\item{P.start}{The starting values for the transition matrix.  A user should
74provide a square matrix with dimension equal to the number of states.  By
75default, draws from the \code{Beta(0.9, 0.1)} are used to construct a proper
76transition matrix for each raw except the last raw.}
77
78\item{b0}{The prior mean of \eqn{\beta}.  This can either be a scalar
79or a column vector with dimension equal to the number of betas. If this
80takes a scalar value, then that value will serve as the prior mean for all
81of the betas.}
82
83\item{B0}{The prior precision of \eqn{\beta}.  This can either be a
84scalar or a square matrix with dimensions equal to the number of betas.  If
85this takes a scalar value, then that value times an identity matrix serves
86as the prior precision of beta. Default value of 0 is equivalent to an
87improper uniform prior for beta.}
88
89\item{a}{\eqn{a} is the shape1 beta prior for transition probabilities.
90By default, the expected duration is computed and corresponding a and b
91values are assigned. The expected duration is the sample period divided by
92the number of states.}
93
94\item{b}{\eqn{b} is the shape2 beta prior for transition probabilities.
95By default, the expected duration is computed and corresponding a and b
96values are assigned. The expected duration is the sample period divided by
97the number of states.}
98
99\item{marginal.likelihood}{How should the marginal likelihood be calculated?
100Options are: \code{none} in which case the marginal likelihood will not be
101calculated, and \code{Chib95} in which case the method of Chib (1995) is
102used.}
103
104\item{gamma.fixed}{1 if users want to constrain \eqn{\gamma} values
105to be constant. By default, \eqn{\gamma} values are allowed to vary
106across regimes.}
107
108\item{...}{further arguments to be passed}
109}
110\value{
111An mcmc object that contains the posterior sample.  This object can
112be summarized by functions provided by the coda package.  The object
113contains an attribute \code{prob.state} storage matrix that contains the
114probability of \eqn{state_i} for each period, the log-likelihood of
115the model (\code{loglike}), and the log-marginal likelihood of the model
116(\code{logmarglike}).
117}
118\description{
119This function generates a sample from the posterior distribution of an
120ordered probit regression model with multiple parameter breaks. The function
121uses the Markov chain Monte Carlo method of Chib (1998).  The user supplies
122data and priors, and a sample from the posterior distribution is returned as
123an mcmc object, which can be subsequently analyzed with functions provided
124in the coda package.
125}
126\details{
127\code{MCMCoprobitChange} simulates from the posterior distribution of an
128ordinal probit regression model with multiple parameter breaks. The
129simulation of latent states is based on the linear approximation method
130discussed in Park (2011).
131
132The model takes the following form:
133
134\deqn{\Pr(y_t = 1) = \Phi(\gamma_{c, m} - x_i'\beta_m) - \Phi(\gamma_{c-1, m} - x_i'\beta_m)\;\; m = 1, \ldots, M}
135
136Where \eqn{M} is the number of states, and \eqn{\gamma_{c, m}} and
137\eqn{\beta_m} are paramters when a state is \eqn{m} at \eqn{t}.
138
139We assume Gaussian distribution for prior of \eqn{\beta}:
140
141\deqn{\beta_m \sim \mathcal{N}(b_0,B_0^{-1}),\;\; m = 1, \ldots, M}
142
143And:
144
145\deqn{p_{mm} \sim \mathcal{B}eta(a, b),\;\; m = 1, \ldots, M}
146
147Where \eqn{M} is the number of states.
148
149Note that when the fitted changepoint model has very few observations in any
150of states, the marginal likelihood outcome can be ``nan," which indicates
151that too many breaks are assumed given the model and data.
152}
153\examples{
154
155set.seed(1909)
156N <- 200
157x1 <- rnorm(N, 1, .5);
158
159## set a true break at 100
160z1 <- 1 + x1[1:100] + rnorm(100);
161z2 <- 1 -0.2*x1[101:200] + rnorm(100);
162z <- c(z1,  z2);
163y <- z
164
165## generate y
166y[z < 1] <- 1;
167y[z >= 1 & z < 2] <- 2;
168y[z >= 2] <- 3;
169
170## inputs
171formula <- y ~ x1
172
173## fit multiple models with a varying number of breaks
174out1 <- MCMCoprobitChange(formula, m=1,
175      	mcmc=100, burnin=100, thin=1, tune=c(.5, .5), verbose=100,
176     	b0=0, B0=0.1, marginal.likelihood = "Chib95")
177out2 <- MCMCoprobitChange(formula, m=2,
178      	mcmc=100, burnin=100, thin=1, tune=c(.5, .5, .5), verbose=100,
179     	b0=0, B0=0.1, marginal.likelihood = "Chib95")
180
181## Do model comparison
182## NOTE: the chain should be run longer than this example!
183BayesFactor(out1, out2)
184
185## draw plots using the "right" model
186plotState(out1)
187plotChangepoint(out1)
188
189}
190\references{
191Jong Hee Park. 2011. ``Changepoint Analysis of Binary and
192Ordinal Probit Models: An Application to Bank Rate Policy Under the Interwar
193Gold Standard."  \emph{Political Analysis}. 19: 188-204. <doi:10.1093/pan/mpr007>
194
195Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011.  ``MCMCpack:
196Markov Chain Monte Carlo in R.'', \emph{Journal of Statistical Software}.
19742(9): 1-21.  \doi{10.18637/jss.v042.i09}.
198
199Siddhartha Chib. 1998. ``Estimation and comparison of multiple change-point
200models.'' \emph{Journal of Econometrics}. 86: 221-241.
201}
202\seealso{
203\code{\link{plotState}}, \code{\link{plotChangepoint}}
204}
205\keyword{models}
206