Model selection with this model of a large number of components

I have a discrete time Markov Chain $$\{X_n: n \in \mathbb{N}_0\}$$ with unknown transition matrix $$P \in \mathbb{R}^{M \times M}$$ on the state space $$\mathcal{S}_X = \{1,2, \dots, M\}$$, with $$M \geq 2$$. At each time step, the random variable $$Y_n$$ is defined by $$Y_n = g(X_n)$$, where $$g: \mathcal{S}_X \to \{0,1\}$$ is a deterministic function.

I can work out the emission probabilities $$\mathbb{P}(Y_n = 0|X_n)$$ and $$\mathbb{P}(Y_n=1|X_n)$$, given $$P$$, i.e. these probabilities are known and computable functions of the entries of $$P$$.

Let's say I now observe $$Z_n = \sum_{i =1}^k Y_n^i = \sum_{i =1}^k g(X_n^i)$$, where $$k \geq 1$$ is unknown, where the sum is taken over $$k$$ iid Markov Chains $$X_n^i$$. Note that while I can write down the likelihood for $$Z_1, Z_2, \dots$$, the function is intractable for even moderately small values of $$k$$.

I am interested in estimating both $$P$$ and $$k$$ from my data $$Z_n$$. I know that I can use some sort of simulation based particle filter to estimate $$P$$ from the data in $$Z$$, however, I am wondering how I can estimate $$k$$. The value for $$k$$ is typically large (>1000). I am not sure that this problem requires advanced/transdimensional MCMC methods (e.g. reversible jump), since $$P$$ is independent of $$k$$, but I am not sure how to proceed as $$k$$ is large. It seems inefficient to fix $$k$$ and fit $$P$$ to each corresponding model, and then performing Bayesian model selection.

Does anyone know of any quick but accurate methods to perform model selection?

Thanks!