---
title: "Simulation de chaînes de Markov"
author: S. Robin
date: septembre 2022
output: html_document
# output: pdf_document
# geometry: margin=1.5cm
# classoption: a4paper
---

# Chaînes de Markov à temps discret

On considère une chaîne de Markov $\{X_t\}_{t = 0, 1, \dots}$ à $K=3$ états et de matrice de transition
$$
\pi = \left[ \begin{array}{ccc} .6 & .3 & .1 \\ .1 & .8 & .1 \\ .4 & .4 & .2 \end{array} \right].
$$
```{r cmProbTrans, echo=TRUE}
pi <- matrix(c(.6, .3, .1, .1, .8, .1, .4, .4, .2), 3, 3, byrow=TRUE)
print(pi)
K <- ncol(pi)
```

1 - Déterminer sa distribution stationnaire $\mu$.

```{r cmDistStat, echo=TRUE}
eig <- eigen(t(pi))
muStat <- eig$vectors[, 1]; muStat <- muStat/sum(muStat)
print(muStat)
```

2 - Partant du premier état, déterminer la probabilié que la chaînes se trouve dans chacun des états après $t = 1, \dots 99$ intervalles de temps.

```{r cmProbEtat, echo=TRUE}
n <- 1e2; muInit <- c(1, 0, 0)
probEtat <- matrix(NA, n, K); probEtat[1, ] <- muInit
for(t in 2:n){probEtat[t, ] <- probEtat[t-1, ]%*%pi}
plot(0, 0, xlim=c(0, n), ylim=c(0, 1), col=0, xlab='', ylab='') # Crée un graphe vide
for(k in 1:K){lines(probEtat[, k], col=k, lwd=2)}
abline(h = muStat, col=(1:K), lty=2, lwd=2)
```

3 - Simuler une trajectoire de la chaine $\{X_t\}_{t = 0, 1, \dots}$.

```{r cmSimul, echo=TRUE}
n <- 1e2; etat <- rep(NA, n); etat[1] <- 1
for(t in 2:n){etat[t] <- sample(1:K, 1, prob=pi[etat[t-1], ])}
plot(etat, type='b')
```

4 - Comparer les fréquences des visites dans les différents états avec la distribution stationnaire.

```{r cmFreqVisites, echo=TRUE}
rbind(muStat, table(etat)/n)
```

# Chaînes de Markov à temps continu

On considère un processus de Markov à temps continu $\{X(t)\}_{t \geq 0}$ à $K=3$ états et de matrice de taux de transision
$$
\rho = \left[ \begin{array}{ccc} - & 1/3 & 1/6 \\ 1/2 & - & 1/2 \\ 1/5 & 1/5 & - \end{array} \right].
$$
```{r pmTauxTrans, echo=TRUE}
rho <- matrix(c(0, 1/3, 1/6, 1/2, 0, 1/2, 1/5, 1/5, 0), 3, 3, byrow=TRUE)
print(rho)
K <- ncol(rho)
```

1 - Déterminer son générateur infinitésimal $Q$ et sa distribution stationnaire $\mu$.

```{r pmGenerateur, echo=TRUE}
Q <- rho; diag(Q) <- -rowSums(Q)
cat('Q =')
print(Q)
eig <- eigen(t(Q))
muStat <- eig$vectors[, K]; muStat <- muStat/sum(muStat)
cat('mu =')
print(muStat)
```

2 - Simuler une suite d'états visités par ce processus.

```{r pmSimulEtat, echo=TRUE}
pi <- rho / rowSums(rho)
n <- 1e3; etat <- rep(NA, n); etat[1] <- 1
for(t in 2:n){etat[t] <- sample(1:K, 1, prob=pi[etat[t-1], ])}
```

3 - Simuler les temps de séjours dans chacun des états successivement visités.

```{r pmSimulTemps, echo=TRUE}
duree <- rexp(n-1)
duree <- duree/(-diag(Q)[etat[-n]])
temps <- c(0, cumsum(duree))
head(cbind(etat, temps))
plot(temps, etat, type='s')
```

4 - Comparer la proportion de temps passée dans chacun des états avec la distribution stationnaire $\mu$.

```{r pmTempsSejour, echo=TRUE}
sejour <- sapply(1:K, function(k){sum(duree[which(etat[-n]==k)])})
rbind(muStat, sejour/sum(sejour))
```