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