1. Definiciones básicas

Sea \(\{X_t\}_{t\ge 0}\) una CMTC en espacio de estados finito \(S=\{1,\dots,m\}\). La matriz generadora \(Q=(q_{ij})\) cumple:

La matriz de transición \(P(t)=(p_{ij}(t))\) está dada por \[ p_{ij}(t)=\Pr\{X_t=j\mid X_0=i\},\qquad P(0)=I. \]

2. Ecuaciones de Chapman–Kolmogórov (CK)

Para todo \(s,t\ge 0\): \[ P(t+s)=P(t)\,P(s)=P(s)\,P(t). \]

Demostración (borrador): Para \(i,j\in S\), \[ p_{ij}(t+s)=\Pr\{X_{t+s}=j\mid X_0=i\} = \sum_{k\in S}\Pr\{X_{t+s}=j,X_t=k\mid X_0=i\} = \sum_k p_{ik}(t)\,p_{kj}(s), \] usando la propiedad de Markov y homogeneidad temporal. En forma matricial, esto es \(P(t+s)=P(t)P(s)\) y, por simetría de los roles de \(t,s\), también \(=P(s)P(t)\).

3. Ecuaciones diferenciales de Kolmogórov

Derivando CK en \(s\) y evaluando en \(s=0\) se obtienen:

Para CMTC homogéneas en espacio finito, ambas se satisfacen y la solución común es \[ P(t)=e^{Qt} \;=\; \sum_{n=0}^\infty \frac{(Qt)^n}{n!}. \]

4. Distribución en el tiempo, estacionaria y límite

Sea \(\pi(t)\) el vector fila de distribución en tiempo \(t\). Si \(\pi(0)=\mu\), entonces \[ \pi(t)=\mu\,P(t) = \mu\,e^{Qt}. \]

Una distribución estacionaria \(\pi\) satisface \[ \pi Q = 0,\qquad \pi\mathbf{1}=1,\quad \pi\ge 0. \]

Si la CMTC es irreducible y positivamente recurrente (ergódica), entonces existe única \(\pi\) y además ocurre convergencia al equilibrio: \[ \lim_{t\to\infty} P(t) = \mathbf{1}\,\pi \quad\Rightarrow\quad \lim_{t\to\infty}\pi(0)P(t)=\pi. \]

5. Ejemplo: proceso nacimiento–muerte finito

Considere \(S=\{0,1,2\}\) con nacimientos tasa \(\lambda\) y muertes tasa \(\mu\): \[ Q=\begin{pmatrix} -\lambda & \lambda & 0\\ \mu & -(\lambda+\mu) & \lambda\\ 0 & \mu & -\mu \end{pmatrix}. \]

6. Cálculo de \(P(t)=e^{Qt}\) y simulación en R

# Paquetes
if (!require(expm)) install.packages("expm")
## Cargando paquete requerido: expm
## Warning: package 'expm' was built under R version 4.4.3
## Cargando paquete requerido: Matrix
## 
## Adjuntando el paquete: 'expm'
## The following object is masked from 'package:Matrix':
## 
##     expm
library(expm)

# Parámetros del ejemplo
lambda <- 2; mu <- 1
Q <- rbind(
  c(-lambda, lambda, 0),
  c(mu, -(lambda+mu), lambda),
  c(0, mu, -mu)
)

# Matriz de transición a t
t <- 1.5
P_t <- expm(Q * t)
round(P_t, 4)
##        [,1]   [,2]   [,3]
## [1,] 0.2016 0.3092 0.4892
## [2,] 0.1546 0.2916 0.5538
## [3,] 0.1223 0.2769 0.6008
# Distribución inicial mu (por ejemplo, empezamos en 0 con prob 1)
mu0 <- c(1,0,0)           # vector fila
pi_t <- mu0 %*% P_t       # distribución en t
round(pi_t, 4)
##        [,1]   [,2]   [,3]
## [1,] 0.2016 0.3092 0.4892
set.seed(123)

simulate_ctmc <- function(Q, x0=1, T_end=10){
  Q <- as.matrix(Q)
  S <- nrow(Q)
  q_out <- -diag(Q)                  # tasas de salida
  t <- 0; x <- x0
  times <- c(0); states <- c(x0)
  while (t < T_end){
    rate <- q_out[x]
    if (rate <= 0) break
    dt <- rexp(1, rate = rate)
    if (t + dt > T_end) {
      t <- T_end; times <- c(times, t); states <- c(states, x); break
    } else {
      t <- t + dt
      probs <- ifelse(1:S != x, Q[x,]/rate, 0)
      probs[x] <- 0
      probs <- probs / sum(probs)
      x <- sample(1:S, 1, prob = probs)
      times <- c(times, t); states <- c(states, x)
    }
  }
  data.frame(time = times, state = states)
}

path <- simulate_ctmc(Q, x0 = 1, T_end = 50)
head(path, 10)
##         time state
## 1  0.0000000     1
## 2  0.4217286     2
## 3  0.8647469     1
## 4  0.8805356     2
## 5  0.9148256     3
## 6  1.2290529     2
## 7  2.1377984     1
## 8  2.8455223     2
## 9  2.9740579     3
## 10 3.3511757     2
if (!require(deSolve)) install.packages("deSolve")
## Cargando paquete requerido: deSolve
## Warning: package 'deSolve' was built under R version 4.4.3
library(deSolve)

# Vectorizamos por filas P(t); dimension m x m
forward_ode <- function(t, y, parms){
  m <- parms$m
  Q <- parms$Q
  P <- matrix(y, nrow = m, byrow = TRUE)
  dP <- P %*% Q
  list(as.vector(t(dP)))  # devolver vector fila-por-fila
}

m <- nrow(Q)
y0 <- as.vector(t(diag(m)))  # I vectorizado por filas
sol <- ode(y = y0, times = seq(0, 2, by = 0.1), func = forward_ode, parms = list(m=m, Q=Q))
# Extraer P(t*) y comparar con expm
t_star <- 1.5
idx <- which.min(abs(sol[,"time"] - t_star))
P_num <- matrix(sol[idx, -1], nrow = m, byrow = TRUE)
cbind(num = round(P_num,4), expm = round(expm(Q*t_star),4))
##        [,1]   [,2]   [,3]   [,4]   [,5]   [,6]
## [1,] 0.2016 0.3092 0.4892 0.2016 0.3092 0.4892
## [2,] 0.1546 0.2916 0.5538 0.1546 0.2916 0.5538
## [3,] 0.1223 0.2769 0.6008 0.1223 0.2769 0.6008