T <- 30 t <- 0 l <- 1 m <- 1/2 Q <- matrix(c(-l, m, 0, 0, 0, l, -l-m, m, 0, 0, 0, l, -l-m, m, 0, 0, 0, l, -l-m, m, 0, 0, 0, l, -m), 5, 5) alpha <- c(1/2, 0, 0, 0, 1/2) # Computing the entries of the embedded Markov chain J <- matrix(0, 5, 5) for(i in seq(1,5,1)){ for(j in seq(1,5,1)){ if(i == j){ J[i,j] = 0; } else{ J[i,j] = -Q[i,j] / Q[i,i] } } } # Simulation of X_{0} X <- min(which(runif(1) < cumsum(alpha))) jump_times <- c(0) jump_chain <- c(X) while(t < T){ # Simulation of the time until the next jump occurs t <- t + log(runif(1)) / Q[X,X] print(t) if(t > T){ break } else{ jump_times <- c(jump_times, t) # Simulation of the next state X <- min(which(runif(1) < cumsum(J[X,]))) jump_chain <- c(jump_chain, X) } } plot(jump_times, jump_chain, type="s")