matrix_question.R

# Make a 3D transition matrix T where z-dimension is time
T_states <- 3
T_terms <- 4
T <- array(runif(T_states * T_states * Tz)/1.5, dim = c(T_states, T_states, Tz))

# Normalise by rows 
for (z in seq(Tz)) {
  T[, , z] <- T[, , z] / apply(T[, , z], 1, sum)
}

# Starting state
start <- matrix(0, nrow = 1, ncol = T_states)
start[1] <- 1 

# Process across terms
end <- start
for (term in seq(T_terms)) {
  end <- end %*% T[, , term]
}