# generate a matrix that contains sample paths and classify it as a "dsdtmc" object simulate_MC_sample_paths <- function(transition_probability_matrix, number_sample_paths, initial_states, number_iterations) { sample_paths <- matrix(data = NA, ncol = number_sample_paths, nrow = number_iterations) for (j in 1:number_sample_paths) { sample_paths[1, j] <- initial_states[j] for (i in 2:number_iterations) { ##### Your Code Starts Here #####  p <- transition_probability_matrix[sample_paths[, j], ] ##### Your Code Ends Here #####  sample_paths[i, j] <- which(rmultinom(n = 1, size = 1, prob = p) == 1) } } # return structure(sample_paths, class = "dsdtmc") # assign sample_paths to the "dsdtmc" class } # add a method to the generic function plot() for the "dsdtmc" class plot.dsdtmc <- function(sample_paths) { pch <- 16 cex <- 0.5 plot(sample_paths[, 1], type = "l", col = 1, xlab = "Time", ylab = "State") points(sample_paths[, 1], pch = pch, cex = cex, col = 1) if (ncol(sample_paths) > 1) { for (j in 2:ncol(sample_paths)) { lines(sample_paths[, j], col = j) points(sample_paths[, j], pch = pch, cex = cex, col = j) } } }