## Section 6.9 notes

## Approximating Transition probabilities

########################################################
# function for calculating matrix powers, by Alberto Monteiro
# Copy and paste this into R first:
#
matpowfast <- function(mat, n)
{
  if (n == 1) return(mat)
  result <- diag(1, ncol(mat))
  while (n > 0) {
    if (n %% 2 != 0) {
      result <- result %*% mat
      n <- n - 1
    }
    mat <- mat %*% mat
    n <- n / 2
  }
  return(result)
}
#
########################################################

# Example 1 (Two-window restaurant)

my.R <- matrix(c(-1/2, 1/2, 0, 0, -1, 1, 1/3, 0, -1/3), nrow=3, ncol=3, byrow=T)

I.3 <- diag(3)

# Method 1:

my.n <- 128

P.5 <- matpowfast(mat = I.3 + my.R*(5/my.n), n = my.n)

print(P.5)

# Method 2:

my.n <- 128

P.5 <- matpowfast(mat = solve(I.3 - my.R*(5/my.n)), n = my.n)

print(P.5)

# The methods give similar answers.