# Introduction to Continuous Time Markov Chains

``````from symbulate import *

from matplotlib import pyplot as plt``````

## Seagull

``````Q = [[-1, 0.4, 0.6],
[0.8 / 3,  -1 / 3, 0.2 / 3],
[0.9 / 5, 0.1 / 5,  -1 / 5]]

pi0 = [1, 0, 0]

states = [1, 2, 3]

X = ContinuousTimeMarkovChain(Q, pi0, states)``````
``````plt.figure();
X.sim(1).plot()
plt.show();``````

``````plt.figure();
X.sim(3).plot()
plt.show();``````

``````plt.figure();
X[1 / 60].sim(10000).plot()
plt.show();``````

``````plt.figure();
X[0.5].sim(10000).plot()
plt.show();``````

``````plt.figure();
X[1].sim(10000).plot()
plt.show();``````

``````plt.figure();
X[5].sim(10000).plot()
plt.show();``````

``````plt.figure();
X[60].sim(10000).plot()
plt.show();``````

Times between jumps

``````P = ContinuousTimeMarkovChainProbabilitySpace(Q, pi0, states)
W = RV(P, interarrival_times)``````
``````plt.figure();
W[0].sim(10000).plot()
plt.show();``````

``````plt.figure();
W[1].sim(10000).plot()
plt.show();``````

``````plt.figure();
(W[0] & W[1]).sim(10000).plot()
plt.show();``````

A few paths from scratch

``````states = 1:3 # state space

pi0 = c(1, 0, 0) # initial distribution

P = rbind(
c(0, 0.6, 0.4),
c(0.8, 0, 0.2),
c(0.9, 0.1, 0)
) # transition matrix of embedded DTMC

lambda = c(1, 1 / 3, 1 / 5)  # transition rates for each state

Nruns = 3 # number of paths to simulate
Njump = 10  # number of transitions to sim for each path

Xt = matrix(rep(NA, Nruns * (Njump + 1)), nrow = Njump + 1) # +1 to include time 0
Wn = Xt

for (r in 1:Nruns){
Xt[1, r] = sample(states, 1, prob = pi0) # initial state
Wn[1, r] = 0 # time 0
for (n in 2:(Njump+1)){
Tn = rexp(1) / lambda[Xt[n - 1, r]] # time between jumps, exp(lambda_i)
Xt[n, r] = sample(states, 1, prob = P[Xt[n - 1, r], ]) # new state
Wn[n,r] = Wn[n - 1, r] + Tn # time of nth jump
}
}
# Plot the paths
Npaths = Nruns
plot(range(Wn),
type = 'n',
ylab = "X(t)", xlab = "t",
main = paste(Npaths," paths of a CTMC"))
colors = c("black","orange","skyblue")

for (r in 1:Npaths){
segments(y0 = Xt[-(Njump + 1), r] + plotadj[r],
x0 = Wn[-(Njump + 1), r],
y1 = Xt[-(Njump + 1), r] + plotadj[r],
x1 = Wn[-1, r], col = colors[r], lty = r, lwd = 2)
points(Wn[, r], Xt[,r] + plotadj[r], pch = 16, col = colors[r])
points(Wn[-1, r], Xt[-(Njump + 1), r] + plotadj[r], col = colors[r])
}``````

## M/M/1 Queue: Sample Path

``````lambda = 0.2

mu = 0.5

n_jumps = 20

X_t = rep(NA, n_jumps + 1)
T_n = rep(NA, n_jumps + 1)

X_t[1] = 0
T_n[1] = 0

for (n in 2:n_jumps){
if (X_t[n - 1] == 0){
T_n[n] = T_n[n - 1] + rexp(1) / lambda
X_t[n] = X_t[n - 1] + 1
} else {
T_n[n] = T_n[n - 1] + rexp(1) / (lambda + mu)
X_t[n] = X_t[n - 1] + sample(c(-1, 1), 1, prob = c(mu, lambda))
}
}

plot(T_n, X_t,
type = "s",
xlab = "t", ylab = "X(t)",
main = "Sample path of M/M/1 queue")``````

## M/M/s Queue: Sample path

``````lambda = 2.3

mu = 0.5

s = 5

n_jumps = 20

X_t = rep(NA, n_jumps + 1)
W_n = rep(NA, n_jumps + 1)

X_t[1] = 0
W_n[1] = 0

for (n in 2:n_jumps){
W_n[n] = W_n[n - 1] + rexp(1) / (lambda + mu * min(s, X_t[n - 1]))
X_t[n] = X_t[n - 1] + sample(c(-1, 1), 1,
prob = c(mu * min(s, X_t[n - 1]), lambda))
}

plot(W_n, X_t,
type = "s",
xlab = "t", ylab = "X(t)",
main = paste("Sample path of M/M/s queue with s = ", s))``````