1% Generated by roxygen2: do not edit by hand
2% Please edit documentation in R/MCMCregressChange.R
3\name{MCMCregressChange}
4\alias{MCMCregressChange}
5\title{Markov Chain Monte Carlo for a linear Gaussian Multiple Changepoint Model}
6\usage{
7MCMCregressChange(
8  formula,
9  data = parent.frame(),
10  m = 1,
11  b0 = 0,
12  B0 = 0,
13  c0 = 0.001,
14  d0 = 0.001,
15  sigma.mu = NA,
16  sigma.var = NA,
17  a = NULL,
18  b = NULL,
19  mcmc = 1000,
20  burnin = 1000,
21  thin = 1,
22  verbose = 0,
23  seed = NA,
24  beta.start = NA,
25  P.start = NA,
26  random.perturb = FALSE,
27  WAIC = FALSE,
28  marginal.likelihood = c("none", "Chib95"),
29  ...
30)
31}
32\arguments{
33\item{formula}{Model formula.}
34
35\item{data}{Data frame.}
36
37\item{m}{The number of changepoints.}
38
39\item{b0}{The prior mean of \eqn{\beta}.  This can either be a scalar
40or a column vector with dimension equal to the number of betas. If this
41takes a scalar value, then that value will serve as the prior mean for all
42of the betas.}
43
44\item{B0}{The prior precision of \eqn{\beta}.  This can either be a
45scalar or a square matrix with dimensions equal to the number of betas.  If
46this takes a scalar value, then that value times an identity matrix serves
47as the prior precision of beta. Default value of 0 is equivalent to an
48improper uniform prior for beta.}
49
50\item{c0}{\eqn{c_0/2} is the shape parameter for the inverse Gamma
51prior on \eqn{\sigma^2} (the variance of the disturbances). The
52amount of information in the inverse Gamma prior is something like that from
53\eqn{c_0} pseudo-observations.}
54
55\item{d0}{\eqn{d_0/2} is the scale parameter for the inverse Gamma
56prior on \eqn{\sigma^2} (the variance of the disturbances). In
57constructing the inverse Gamma prior, \eqn{d_0} acts like the sum of
58squared errors from the \eqn{c_0} pseudo-observations.}
59
60\item{sigma.mu}{The mean of the inverse Gamma prior on
61\eqn{\sigma^2}.  \eqn{sigma.mu} and
62\eqn{sigma.var} allow users to choose the inverse Gamma prior by
63choosing its mean and variance.}
64
65\item{sigma.var}{The variacne of the inverse Gamma prior on
66\eqn{\sigma^2}.  \eqn{sigma.mu} and
67\eqn{sigma.var} allow users to choose the inverse Gamma prior by
68choosing its mean and variance.}
69
70\item{a}{\eqn{a} is the shape1 beta prior for transition probabilities.
71By default, the expected duration is computed and corresponding a and b
72values are assigned. The expected duration is the sample period divided by
73the number of states.}
74
75\item{b}{\eqn{b} is the shape2 beta prior for transition probabilities.
76By default, the expected duration is computed and corresponding a and b
77values are assigned. The expected duration is the sample period divided by
78the number of states.}
79
80\item{mcmc}{The number of MCMC iterations after burnin.}
81
82\item{burnin}{The number of burn-in iterations for the sampler.}
83
84\item{thin}{The thinning interval used in the simulation.  The number of
85MCMC iterations must be divisible by this value.}
86
87\item{verbose}{A switch which determines whether or not the progress of the
88sampler is printed to the screen.  If \code{verbose} is greater than 0 the
89iteration number, the \eqn{\beta} vector, and the error variance are
90printed to the screen every \code{verbose}th iteration.}
91
92\item{seed}{The seed for the random number generator.  If NA, the Mersenne
93Twister generator is used with default seed 12345; if an integer is passed
94it is used to seed the Mersenne twister.  The user can also pass a list of
95length two to use the L'Ecuyer random number generator, which is suitable
96for parallel computation.  The first element of the list is the L'Ecuyer
97seed, which is a vector of length six or NA (if NA a default seed of
98\code{rep(12345,6)} is used).  The second element of list is a positive
99substream number. See the MCMCpack specification for more details.}
100
101\item{beta.start}{The starting values for the \eqn{\beta} vector.
102This can either be a scalar or a column vector with dimension equal to the
103number of betas.  The default value of of NA will use the MLE estimate of
104\eqn{\beta} as the starting value.  If this is a scalar, that value
105will serve as the starting value mean for all of the betas.}
106
107\item{P.start}{The starting values for the transition matrix.  A user should
108provide a square matrix with dimension equal to the number of states.  By
109default, draws from the \code{Beta(0.9, 0.1)} are used to construct a proper
110transition matrix for each raw except the last raw.}
111
112\item{random.perturb}{If TRUE, randomly sample hidden states whenever
113regularly sampled hidden states have at least one single observation state
114(SOS). SOS is a sign of overfitting in non-ergodic hidden Markov models.}
115
116\item{WAIC}{Compute the Widely Applicable Information Criterion (Watanabe
1172010).}
118
119\item{marginal.likelihood}{How should the marginal likelihood be calculated?
120Options are: \code{none} in which case the marginal likelihood will not be
121calculated, and \code{Chib95} in which case the method of Chib (1995) is
122used.}
123
124\item{...}{further arguments to be passed}
125}
126\value{
127An mcmc object that contains the posterior sample.  This object can
128be summarized by functions provided by the coda package.  The object
129contains an attribute \code{prob.state} storage matrix that contains the
130probability of \eqn{state_i} for each period, the log-likelihood of
131the model (\code{loglike}), and the log-marginal likelihood of the model
132(\code{logmarglike}).
133}
134\description{
135This function generates a sample from the posterior distribution of a linear
136Gaussian model with multiple changepoints. The function uses the Markov
137chain Monte Carlo method of Chib (1998).  The user supplies data and priors,
138and a sample from the posterior distribution is returned as an mcmc object,
139which can be subsequently analyzed with functions provided in the coda
140package.
141}
142\details{
143\code{MCMCregressChange} simulates from the posterior distribution of the
144linear regression model with multiple changepoints.
145
146The model takes the following form:
147
148\deqn{y_t=x_t ' \beta_i + I(s_t=i)\varepsilon_{t},\;\; i=1, \ldots, k}
149
150Where \eqn{k} is the number of states and \eqn{I(s_t=i)} is an
151indicator function that becomes 1 when a state at \eqn{t} is
152\eqn{i} and otherwise 0.
153
154The errors are assumed to be Gaussian in each regime:
155
156\deqn{I(s_t=i)\varepsilon_{t} \sim \mathcal{N}(0, \sigma^2_i)}
157
158We assume standard, semi-conjugate priors:
159
160\deqn{\beta_i \sim \mathcal{N}(b_0,B_0^{-1}),\;\; i=1, \ldots, k}
161
162And:
163
164\deqn{\sigma^{-2}_i \sim \mathcal{G}amma(c_0/2, d_0/2),\;\; i=1, \ldots, k}
165
166Where \eqn{\beta_i} and \eqn{\sigma^{-2}_i} are assumed \emph{a
167priori} independent.
168
169The simulation proper is done in compiled C++ code to maximize efficiency.
170}
171\examples{
172
173\dontrun{
174set.seed(1119)
175n <- 100
176x1 <- runif(n)
177true.beta1 <- c(2, -2)
178true.beta2 <- c(0,  2)
179true.Sigma <- c(1, 2)
180true.s <- rep(1:2, each=n/2)
181
182mu1 <- cbind(1, x1[true.s==1])\%*\%true.beta1
183mu2 <- cbind(1, x1[true.s==2])\%*\%true.beta2
184
185y <- as.ts(c(rnorm(n/2, mu1, sd=sqrt(true.Sigma[1])), rnorm(n/2, mu2, sd=sqrt(true.Sigma[2]))))
186formula=y ~ x1
187
188ols1 <- lm(y[true.s==1] ~x1[true.s==1])
189ols2 <- lm(y[true.s==2] ~x1[true.s==2])
190
191## prior
192b0 <- 0
193B0 <- 0.1
194sigma.mu=sd(y)
195sigma.var=var(y)
196
197## models
198model0 <-  MCMCregressChange(formula, m=0, b0=b0, B0=B0, mcmc=100, burnin=100,
199           sigma.mu=sigma.mu, sigma.var=sigma.var, marginal.likelihood="Chib95")
200model1 <-  MCMCregressChange(formula, m=1, b0=b0, B0=B0, mcmc=100, burnin=100,
201           sigma.mu=sigma.mu, sigma.var=sigma.var, marginal.likelihood="Chib95")
202model2 <-  MCMCregressChange(formula, m=2, b0=b0, B0=B0, mcmc=100, burnin=100,
203           sigma.mu=sigma.mu, sigma.var=sigma.var, marginal.likelihood="Chib95")
204model3 <-  MCMCregressChange(formula, m=3, b0=b0, B0=B0, mcmc=100, burnin=100,
205           sigma.mu=sigma.mu, sigma.var=sigma.var, marginal.likelihood="Chib95")
206model4 <-  MCMCregressChange(formula, m=4, b0=b0, B0=B0, mcmc=100, burnin=100,
207           sigma.mu=sigma.mu, sigma.var=sigma.var, marginal.likelihood="Chib95")
208model5 <-  MCMCregressChange(formula, m=5, b0=b0, B0=B0, mcmc=100, burnin=100,
209           sigma.mu=sigma.mu, sigma.var=sigma.var, marginal.likelihood="Chib95")
210
211print(BayesFactor(model0, model1, model2, model3, model4, model5))
212plotState(model1)
213plotChangepoint(model1)
214
215}
216
217}
218\references{
219Jong Hee Park, 2012. ``Unified Method for Dynamic and
220  Cross-Sectional Heterogeneity: Introducing Hidden Markov Panel
221  Models.''  \emph{American Journal of Political Science}.56:
222  1040-1054. <doi: 10.1111/j.1540-5907.2012.00590.x>
223
224Sumio Watanabe. 2010. "Asymptotic equivalence of Bayes cross validation and
225widely applicable information criterion in singular learning theory"
226\emph{Journal of Machine Learning Research}. 11: 3571-3594.
227
228Siddhartha Chib. 1995. "Marginal Likelihood from the Gibbs Output."
229\emph{Journal of the American Statistical Association}. 90: 1313-1321.
230<doi: 10.1016/S0304-4076(97)00115-2>
231
232Siddhartha Chib. 1998. "Estimation and comparison of multiple change-point
233models."  \emph{Journal of Econometrics}. 86: 221-241.
234<doi: 10.1080/01621459.1995.10476635>
235
236Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011.
237``MCMCpack: Markov Chain Monte Carlo in R.'', \emph{Journal of
238Statistical Software}.  42(9): 1-21. \doi{10.18637/jss.v042.i09}.
239}
240\seealso{
241\code{\link{plotState}}, \code{\link{plotChangepoint}}
242}
243\keyword{models}
244