ppmSuite (version 0.3.4)

icp_ppm: Function that fits the multivariate independent product partition change point model

Description

icp_ppm is a function that fits a Bayesian product partition change point model. Each series is treated independently.

Usage

icp_ppm(ydata,
         a0, b0,
         mltypes,
         thetas,
         nburn, nskip, nsave,
         verbose = FALSE)

Value

The function returns a list containing arrays filled with MCMC iterates corresponding to model parameters. In order to provide more detail, in what follows let \(M\) be the number of MCMC iterates collected. The output list contains the following:

• C. An \(M \times \{L(n - 1)\}\) matrix containing MCMC iterates associated with each series indicators of a change point. The \(m\)th row in C is divided into \(L\) blocks; the first \((n - 1)\) change point indicators for time series 1, the next \((n - 1)\) change point indicators for time series 2, and so on.

• P. An \(M \times \{L(n - 1)\}\) matrix containing MCMC iterates associated with each series probability of a change point. The \(m\)th row in P is divided into \(L\) blocks; the first \((n - 1)\) change point probabilities for time series 1, the next \((n - 1)\) change point probabilities for time series 2, and so on.

Arguments

ydata

An \(L \times n\) data matrix, where \(L\) is the number of time series and \(n\), the number of time points.

a0

Vector of dimension \(L\) with shape 1 Beta parameters (see Details).

b0

Vector of dimension \(L\) with shape 2 Beta parameters (see Details).

mltypes

Type of marginal likelihood. Currently only available is:

• mltypes = 1. Observations within a block are conditionally independent \(Normal(\mu, \sigma^2)\) variates with mean \(\mu\) and variance \(\sigma^2\). The desired marginal likelihood is obtained after integrating \((\mu, \sigma^2)\) with respect to a \(Normal-Inverse-Gamma(\mu_0, \kappa_0, \alpha_0, \beta_0)\) prior.

thetas

An \(L \times q\) matrix containing hyperparameters associated with the marginal likelihood. The number of rows \((L)\) corresponds to the number of series. The number of columns \((q)\) depend on the marginal likelihood:

• If mltypes = 1, then \(q = 4\) and thetas equals the hyperparameter \((\mu_{0}, \kappa_{0}, \alpha_{0}, \beta_{0})\) of the Normal-Inverse-Gamma prior.

nburn

The number of initial MCMC iterates to be discarded as burn-in.

nskip

The amount to thinning that should be applied to the MCMC chain.

nsave

Then number of MCMC iterates to be stored.

verbose

Logical indicating whether to print to screen the MCMC progression. The default value is verbose = FALSE.

Details

As described in Barry and Hartigan (1992) and Loschi and Cruz (2002), for each time series \(\boldsymbol{y}_{i} = (y_{i,1}, \ldots , y_{i,n})'\):

$$\boldsymbol{y}_{i} \mid \rho_{i} \sim \prod_{j = 1}^{b_{i}}\mathcal{F}(\boldsymbol{y}_{i,j} \mid \boldsymbol{\theta}_{i})$$

$$\rho_{i} \mid p_{i} \sim p_{i}^{b_{i} - 1}(1 - p_{i})^{n - b_{i}}$$

$$p_{i} \sim Beta(a_{i,0}, b_{i,0}).$$

Here, \(\rho_{i} = \{S_{i,1}, \ldots , S_{i,b_{i}}\}\) is a partition of the set \(\{1, \ldots , n\}\) into \(b_{i}\) contiguous blocks, and \(\boldsymbol{y}_{i,j} = (y_{i,t} : t \in S_{i,j})'\). Also, \(\mathcal{F}( \cdot \mid \boldsymbol{\theta}_{i})\) is a marginal likelihood function which depends on the nature of \(\boldsymbol{y}_{i}\), indexed by a hyperparameter \(\boldsymbol{\theta}_{i}\). Notice that \(p_{i}\) is the probability of observing a change point in series \(i\), at each time \(t \in \{2, \ldots , n\}\).

Examples

# \donttest{
# Generate data that has two series, each with 100 observations
y1 <- replicate(25, rnorm(4, c(-1, 0, 1, 2), c(0.1, 0.25, 0.5, 0.75)))
y2 <- replicate(25, rnorm(4, c(2, 1, 0, -2), c(0.1, 0.25, 0.5, 0.75)))
y <- rbind(c(t(y1)), c(t(y2)))
n <- ncol(y)
# Marginal likelihood parameters
thetas <- matrix(1, nrow = 2, ncol = 4)
thetas[1,] <- c(0, 1, 2, 1)
thetas[2,] <- c(0, 1, 2, 1)

# Fit the Bayesian ppm change point model
fit <- icp_ppm(ydata = y,
               a0 = c(1, 1),
               b0 = c(1, 1),
               mltypes = c(1, 1),
               thetas = thetas,
               nburn = 1000, nskip = 1, nsave = 1000)

cpprobsL <- matrix(apply(fit$C,2,mean), nrow=n-1, byrow=FALSE)


# }