The function below counts the initial uniformly chosen state as the first of the visited states, and then generates transitions using the appropriate row of :
simchain <- function(n, P) {
  if (!is.numeric(n) || length(n) != 1L || !is.finite(n) ||
      n < 1 || n != floor(n)) stop("n must be a positive integer")
  if (!is.matrix(P) || !is.numeric(P) || nrow(P) == 0L ||
      nrow(P) != ncol(P) || any(!is.finite(P)) ||
      any(P < 0) || any(abs(rowSums(P) - 1) > 1e-12)) {
    stop("P must be a finite square stochastic matrix")
  }
  k <- nrow(P)
  state <- sim1(rep(1/k, k))
  counts <- numeric(k)
  counts[state] <- 1
  if (n > 1) {
    for (step in seq_len(n - 1)) {
      state <- sim1(P[state, ])
      counts[state] <- counts[state] + 1
    }
  }
  counts / n
}
Conditional on the current state , sim1(P[i, ]) selects with probability , proving the Markov property with the desired transition matrix. Initialization is uniform by the same sampler. Dividing the visit counts by gives the requested Empirical occupation proportions of a finite Markov chain, including a valid result when .