The problem is that even though you may obtain a unbiased estimate $\hat{S}$ of the stochastic matrix $S$, $\hat{S}^k$ is not an unbiased estimate of $S^k$, the k-steps transition matrix.

To account for the convexity, you need to put a prior on your transition, for instance an independent Dirichlet prior on each column. Observing the Markov chain will give you a Bayesian update to that distribution (which is conveniently a conjugate prior). 

You can then sample whole chains, or even marginals within the chain, but every single chain will unfortunately represent a point hypothesis over $S$.