The field of exact analysis of independent discrete outcomes data with covariates has come of age in the recent years. Development of efficient numerical algorithms has been the prime factor fueling the growth. Exact inference on correlated or time-varying discrete data, on the other hand, has garnered little attention. In this paper, we lay the foundation for exact analysis for one class of longitudinal data problems – Markov chain models with covariates. In particular, we focus on a two-state first-order stationary Markov chain model with a multi-level stratification factor and a binary covariate of interest. We show that exact distributions for such models derive from conditional convolutions of elemental Bose–Einstein distributions. We develop a recursive polynomial multiplication algorithm with checks for infeasibility to compute such distributions, and empirically demonstrate its efficiency. We utilize the algorithm to compare exact inference on Markov chain data with the traditional large sample method. Given the transition counts, the latter method is insensitive as to whether the data are generated from a few long chains or many short chains. The exact method takes this feature into account. Our findings point towards a potential problem with large sample inference, namely the use of observed instead of expected information in finite samples. This issue needs further study. It is also seen that with Markov chain data, inferences drawn from the exact method may differ from that from the asymptotic method at transition count values that may not be deemed small.