I am following the Hidden Markov Models guide text for Time Series An Introduction Using R (Walter Zucchini). Chapter 7. Bayesian inference for Poisson-hidden Markov models. Specifically in section 7.2 talk about how to select the right model from the Bayes factor between two models.
Then use a code in rstan to estimate the Poisson-Hidden Markov model parameters. For the base data Earthquakes: Number of major earthquakes (magnitude 7 or greater) in the World, 1900-2006.
library(rstan)
dat <- read.table("http://www.hmms-for-time-series.de/second/data/earthquakes.txt")
stan.data <- list(n=dim(dat)[1], m=2, x=dat$V2)
PHMM <- '
data {
int<lower=0> n; // length of the time series
int<lower=0> x[n]; // data
int<lower=1> m; // number of states
}
parameters{
simplex[m] Gamma[m]; // tpm
positive_ordered[m] lambda; // mean of poisson - ordered
}
model{
vector[m] log_Gamma_tr[m]; // log, transposed tpm
vector[m] lp; // for forward variables
vector[m] lp_p1; // for forward variables
lambda ~ gamma(0.1, 0.01); // assigning exchangeable priors
//(lambdas´s are ordered for sampling purposes)
// transposing tpm and taking the log of each entry
for(i in 1:m)
for(j in 1:m)
log_Gamma_tr[j, i] = log(Gamma[i, j]);
lp = rep_vector(-log(m), m); //
for(i in 1:n) {
for(j in 1:m)
lp_p1[j] = log_sum_exp(log_Gamma_tr[j] + lp) + poisson_lpmf(x[i] | lambda[j]);
lp = lp_p1;
}
target += log_sum_exp(lp);
}'
The output is as follows:
modelo <- stan(model_code = PHMM, data = stan.data, iter = 1000, chains = 1)
print(modelo,digits_summary = 3)
Inference for Stan model: f4fa370704db17d55928f5a217f01f5d.
1 chains, each with iter=1000; warmup=500; thin=1;
post-warmup draws per chain=500, total post-warmup draws=500.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
Gamma[1,1] 0.905 0.003 0.049 0.786 0.881 0.912 0.939 0.973 280 1.004
Gamma[1,2] 0.095 0.003 0.049 0.027 0.061 0.088 0.119 0.214 280 1.004
Gamma[2,1] 0.142 0.003 0.068 0.032 0.093 0.135 0.178 0.308 501 0.998
Gamma[2,2] 0.858 0.003 0.068 0.692 0.822 0.865 0.907 0.968 501 0.998
lambda[1] 15.068 0.053 0.953 13.121 14.497 15.184 15.736 16.641 324 0.998
lambda[2] 25.611 0.082 1.450 22.948 24.675 25.657 26.504 28.830 310 0.998
lp__ -350.196 0.108 1.666 -354.707 -350.908 -349.797 -348.959 -348.187 237 0.998
Samples were drawn using NUTS(diag_e) at Thu Feb 07 14:24:17 2019.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).
7.2 Bayesian estimation of the number of states.
In the Bayesian approach to model selection, the number of states, $m$, is a parameter whose value is assessed from its posterior distribution, $p(m | x^{(T)})$. Computing this posterior distribution is, however, not an easy problem; indeed, it has been described as `notoriously difficult to calculate' (Scott, James and Sugar, 2005). Using $p$ as a general symbol for probability mass or density functions, one has
$$p(m | x^{(T)}) = p(m) p(x^{(T)} | m) / p(x^{(T)}) \propto p(m) p(x^{(T)} | m) \ \ \ \ \ (7.3)$$
where $p(x^{(T)} | m)$ is called the integrated likelihood. If only two models are being compared, the posterior odds are equal to the product of the `Bayes factor' and the prior odds:
$$\frac{p(m_2 | x^{(t)})}{p(m_1 | x^{(t)})} = \frac{p(x^{(T)} | m_2)}{p(x^{(T)} | m_1)} \times \frac{p(m_2)}{p(m_1)} \ \ \ \ \ (7.4)$$
7.2.1 Use of the integrated likelihood
In order to use (7.3) or (7.4) we need to estimate the integrated likelihood
$$p(x^{(T)} | m) = \int p(\theta_m, x^{(T)} | m) d \theta_m \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ = \int p(x^{(T)} | m, \theta_m) p(\theta_m | m) d \theta_m$$
One way of doing so would be to simulate from $p(\theta_m | m)$, the prior distribution of the parameters $\theta_m$ of the m-state model. But it is convenient and especially if the prior is diffuse and more efficient to use a method that requires instead a sample from the posterior distribution, $p(\theta_m | x(T), m)$. Such a method is as follows.
Write the integrated likelihood as
$$\int p(x^{(T)} | m, \theta_m) \frac{p(\theta_m | m)}{p^*(\theta_m)} p*(\theta_m) d \theta_m$$
that is, write it in a form suitable for the use of a sample from some convenient density $p^*(\theta_m)$ for the parameters $\theta_m$. Since we have available a sample $\theta_m^{(j)} \ (j=1,2,...,B)$ from the posterior distribution, we can use that sample; that is, we can take $p^*(\theta_m) = p(\theta_m | x^{(T)}, m)$. Newton and Raftery (1994) therefore suggest inter alia that the integrated likelihood can be estimated by
$$\hat{I} = \sum_{j=1}^B w_j \ p(x^{(T)} | m, \theta_m^{(j)}) \Big/ \sum_{j=1}^B w_j, \ \ \text{where} \ \ w_j = \frac{p(\theta_m^{(j)} | m)}{p(\theta_m^{(j)} | x^{(T)}, m)} \ \ \ \ \ \ (7.5)$$
After some manipulation this simplies to the harmonic mean of the likelihood values of a sample from the posterior,
$$\hat{I} = \left(B^{-1} \sum_{j=1}^B \left( p(x^{(T)} | m, \theta_m^{(j)}) \right)^{-1} \right); \ \ \ \ \ \ (7.6)$$
In section 7.2.2 Zucchini, indicates that another and simpler way to select the best model by parallel sampling.
The truth is that I have no idea how to elaborate the code in R, which allows me to calculate the factor of bays either by method 1. using the integrated likelihood or by method 2 using parallel sampling and thus choosing the most suitable model. Any suggestion is well-received.
New information
A hidden Markov model $\lbrace X_t : t \in \mathbb{N} \rbrace$ is a particular kind of dependent mixture. With $X^{(t)}$ and $C^{(t)}$ representing the histories from time $1$ to time $t$, one can summarize the simplest model of this kind by:
$$P_r (C_t | C^{(t-1)}) = P_r (C_t | C_{t-1}), \ \ \ t = 2,3,... \ \ \ \ \ (2.1)$$
$$P_r(X_t | X^{(t-1)}, C^{(t)}) = P_r(X_t | C_t) \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (2.2)$$
The model consists of two parts: firstly, an unobserved 'parameter process' $\lbrace C_t: t = 1,2,... \rbrace$ satisfying the Markov property; and secondly, the 'state-dependent process' $\brace X_t: t = 1,2,... \rbrace$, in which the distribution of $X_t$ depends only on the current state $C_t$ and not on previous states or observations. This structure is represented by the directed graph in Figure 2.2.
rstancode. $\endgroup$