R/dmc_period.R

In paulnorthrop/stat2003: stat2003: Introduction to Applied Probability at UCL

#' Finding period of each state for a discrete-time Markov chain
#'
#' \code{dmc_period} returns the period of each state for a given discrete-time Markov chain
#'
#'
#' @param MC is an object in class 'stat2003.d'
#'
#' @return \code{dmc_period } returns the period of each state in a number or
#' in 'aperiodic' if period is 1
#'
#' @seealso \code{\link{stat2003.d-class}}  A class type of discrete-time finite
#' state space Markov chain in stat2003 package.
#' @seealso \code{\link{dmc_simu}} simulates a discrete-time Markov chain by returning
#' one possible sequence and a states against steps plot.
#' @seealso \code{\link{dmc_tp}} can calculate transient probabilities at a specific step,
#' and also can give a graph about transient probabilities from step zero to that specific step.
#' @seealso \code{\link{dmc_equi}} returns the equilibrium
#' distribution for a discrete-time Markov chain.
#' @seealso \code{\link{dmc_inv}} returns the invariant
#' distribution for a discrete-time Markov chain.
#' @seealso \code{\link{dmc_irreclass}} focuses on irreducible classes for
#' a given discrete-time Markov chain
#'
#' @examples
#'
#' m <- matrix(c(0, 1, 0, 0,
#'               1, 0, 0, 0,
#'               1/4, 1/4, 1/4,1/4,
#'               0, 0, 1, 0), nr = 4, nc=4, byrow = TRUE)
#' A <- new("stat2003.d", p_start = c(1/2, 0, 0, 1/2), p = m,
#'          statespace = letters[1:4] )
#' dmc_period(A)
#'
#'@export
dmc_period <- function(MC) {
  if (class(MC) != "stat2003.d") {
    stop("MC has to be an object in class 'stat2003.d' !")
  }

  # created a function to find the greatest common divisor
  gcd <- function(x,y) {
    r <- x%%y;
    return(ifelse(r, gcd(y, r), y))
  }
  # created a variable to record period
  period_recorder <- matrix(0, nrow = 1, ncol = nrow(MC@p),
                            dimnames = list("periodicity", MC@statespace))

  for (i in 1 : nrow(MC@p)) {
    n <- NULL
    p_n <- MC@p

    # the next loop is used to find time j when pii > 0,
    # i.e. when the process comes back to the state.
    for (j in 1 : (300 * nrow(MC@p))) {
      if (p_n[i, i] > 0) {
        n <- append(n, j)  # recording j if pii > 0 at time j
      }
      p_n <- p_n %*% MC@p
    }

    # If there are no elements in n
    if (is.null(n) == TRUE) {
      # if pii = 0 for all j>=1, we also say the periodicity of
      # this state is 1 (aperiodic) by convention
      period_recorder[1, i] <- 1
    } else if (is.null(n) == FALSE & length(n) == 1) {
      # if there is only one number in n, the number is the periodicity of this state
      period_recorder[1, i] <- n[1]
    } else {
      # if there are more than one numbers in n,
      # we need to find the greatest common divisor of these numbers
      common_divisor <- gcd(n[1], n[2])
      for (k in 1 : length(n)) {
        common_divisor <- gcd(common_divisor, n[k])
      }
      period_recorder[1, i] <- common_divisor  # recording periodicity of each state
    }
  }

  list(replace(period_recorder, period_recorder == 1, "aperiodic"))
}