The time-series model in the paper is
- multivariate with $N$ time series of duration $T$
- regime-switching, meaning that the model parameters are changing at random times (with indicators $s_t$), with an unknown number $S$ of switches
- clustered, meaning that during one regime (time interval) $s$, the components of the time series are driven by one of $K$ models, with all components in the same cluster $k$ following the same model (sharing the same parameter $\theta_{sk}$
When considering a Bayesian modelling of this setting, simulating the parameters and the latent variables (regimes and clusters) can be done by Gibbs sampling. For the cluster allocation, this can be done conditional on the regimes $\mathbf s^T$ and the parameters, since the clustering is independent between regimes. Each component $n$ of the time series vector can then be allocated to one of the $K$ clusters, conditional on the other components, with a probability
$$\mathbb P(c_{n,s}=\kappa)\propto \eta_{s\kappa} \prod_k \int f(\mathbf y_{ns}|\theta_{s\kappa})\pi(\theta s_k|\mathbf y_{-n,s})\,\text d\theta_{sk}\tag{1}$$
that has integrated out the first level parameters, thanks to the conjugacy of the priors (cf. (A-19) in the Appendix A-3). Since the rhs is available in closed form, the simulation of $c_{n,s}$ is straightforward as a discrete distribution on $\{1,2,\ldots,K\}$, as indicated in John Madden's comment.
calculate the expressions for c = 1,2,3 and store those in a vector.
Then, divide that vector by its sum. Then, sample (1,2,3) with
probability given by that vector.
This is one step of a partly integrated Gibbs sampler that is perfectly valid. Note that, due to this integration of the first level group- and cluster-dependent parameters, the components $\mathbf y_{ms}$ are no longer (conditionally) independent, which is why they all appear in the right hand side of (1).
