Chapter 3: Continuous-Time Markov Processes

Haiyue
24min

Chapter 3: Continuous-Time Markov Processes

3.1 Definition of Continuous-Time Markov Processes

3.1.1 Basic Definition

Let {X(t),t≥0}\{X(t), t \geq 0\} be a stochastic process defined on a continuous time set with state space SS. If for any 0≤t0<t1<⋯<tn<t0 \leq t_0 < t_1 < \cdots < t_n < t and any states i0,i1,…,in,j∈Si_0, i_1, \ldots, i_n, j \in S, we have:

P(X(t)=j∣X(t0)=i0,X(t1)=i1,…,X(tn)=in)=P(X(t)=j∣X(tn)=in)P(X(t) = j | X(t_0) = i_0, X(t_1) = i_1, \ldots, X(t_n) = i_n) = P(X(t) = j | X(t_n) = i_n)

then {X(t)}\{X(t)\} is called a continuous-time Markov process.

3.1.2 Homogeneity

If the transition probability depends only on the time interval and not on the starting time, i.e.: P(X(t+s)=j∣X(s)=i)=P(X(t)=j∣X(0)=i)P(X(t+s) = j | X(s) = i) = P(X(t) = j | X(0) = i)

then it is called a time-homogeneous continuous-time Markov process.

3.1.3 Right Continuity

Continuous-time Markov processes typically assume paths are right-continuous, i.e.: lim⁡s↓tX(s)=X(t)\lim_{s \downarrow t} X(s) = X(t)

This ensures mathematical completeness in the treatment of the process.

3.2 Transition Functions and Generator Matrix

3.2.1 Transition Function

Define the transition function: Pij(t)=P(X(t)=j∣X(0)=i),t≥0P_{ij}(t) = P(X(t) = j | X(0) = i), \quad t \geq 0

Properties:

  1. Initial Condition: Pij(0)=δijP_{ij}(0) = \delta_{ij} (Kronecker delta)
  2. Chapman-Kolmogorov Equation: Pij(s+t)=∑k∈SPik(s)Pkj(t)P_{ij}(s+t) = \sum_{k \in S} P_{ik}(s) P_{kj}(t)
  3. Matrix Form: P(s+t)=P(s)P(t)P(s+t) = P(s)P(t)

3.2.2 Generator Matrix (Q Matrix)

Definition: The generator matrix Q=(qij)Q = (q_{ij}) is defined as: qij=lim⁡h↓0Pij(h)−δijhq_{ij} = \lim_{h \downarrow 0} \frac{P_{ij}(h) - \delta_{ij}}{h}

Physical Meaning:

  • qii=−∑j≠iqijq_{ii} = -\sum_{j \neq i} q_{ij}: Total departure rate from state ii
  • qijq_{ij} (when i≠ji \neq j): Instantaneous transition rate from state ii to state jj

3.2.3 Kolmogorov Differential Equations

Forward Equation: ddtPij(t)=∑k∈SPik(t)qkj\frac{d}{dt} P_{ij}(t) = \sum_{k \in S} P_{ik}(t) q_{kj}

Backward Equation: ddtPij(t)=∑k∈SqikPkj(t)\frac{d}{dt} P_{ij}(t) = \sum_{k \in S} q_{ik} P_{kj}(t)

Matrix Form: ddtP(t)=P(t)Q=QP(t)\frac{d}{dt} P(t) = P(t) Q = Q P(t)

Solution: P(t)=eQt=∑n=0∞(Qt)nn!P(t) = e^{Qt} = \sum_{n=0}^{\infty} \frac{(Qt)^n}{n!}

3.3 Poisson Process

3.3.1 Definition and Properties

A Poisson process {N(t),t≥0}\{N(t), t \geq 0\} is a counting process with the following properties:

  1. Independent Increments: Increments over disjoint time intervals are independent
  2. Stationary Increments: Distribution of increments depends only on interval length
  3. Finiteness: Finite number of events in finite time
🔄 正在渲染 Mermaid 图表...

3.3.2 Mathematical Characteristics

For a Poisson process with intensity λ\lambda:

Probability Mass Function: P(N(t)=n)=(λt)ne−λtn!,n=0,1,2,…P(N(t) = n) = \frac{(\lambda t)^n e^{-\lambda t}}{n!}, \quad n = 0, 1, 2, \ldots

Mean and Variance: E[N(t)]=λt,Var[N(t)]=λtE[N(t)] = \lambda t, \quad \text{Var}[N(t)] = \lambda t

Inter-arrival Times: Let TnT_n be the time of the nn-th event, then the inter-arrival time Sn=Tn−Tn−1S_n = T_n - T_{n-1} follows an exponential distribution with parameter λ\lambda: fSn(s)=λe−λs,s>0f_{S_n}(s) = \lambda e^{-\lambda s}, \quad s > 0

3.3.3 Compound Poisson Process

Definition: Let {N(t)}\{N(t)\} be a Poisson process with intensity λ\lambda, and {Yi}\{Y_i\} be i.i.d. random variables, then: X(t)=∑i=1N(t)YiX(t) = \sum_{i=1}^{N(t)} Y_i

is called a compound Poisson process.

Characteristic Function: ϕX(t)(u)=exp⁡(λt(ϕY(u)−1))\phi_{X(t)}(u) = \exp\left(\lambda t (\phi_Y(u) - 1)\right)

where ϕY(u)\phi_Y(u) is the characteristic function of YiY_i.

3.4 Birth-Death Process

3.4.1 Definition

A birth-death process is a continuous-time Markov process with state space of non-negative integers {0,1,2,…}\{0, 1, 2, \ldots\}, where transitions can only occur between adjacent states:

  • From state nn to state n+1n+1 (birth) with transition rate λn\lambda_n
  • From state nn to state n−1n-1 (death) with transition rate μn\mu_n

3.4.2 Generator Matrix

The generator matrix of a birth-death process has a tridiagonal structure:

-\lambda_0 & \lambda_0 & 0 & 0 & \cdots \\ \mu_1 & -(\lambda_1 + \mu_1) & \lambda_1 & 0 & \cdots \\ 0 & \mu_2 & -(\lambda_2 + \mu_2) & \lambda_2 & \cdots \\ \vdots & \vdots & \vdots & \vdots & \ddots \end{pmatrix}$$ ### 3.4.3 Stationary Distribution For an irreducible birth-death process, the stationary distribution $\pi = (\pi_0, \pi_1, \ldots)$ satisfies: $$\pi Q = 0$$ **Detailed Balance Condition**: $$\pi_n \lambda_n = \pi_{n+1} \mu_{n+1}, \quad n \geq 0$$ **Solution**: $$\pi_n = \pi_0 \prod_{k=1}^n \frac{\lambda_{k-1}}{\mu_k}, \quad n \geq 1$$ where $\pi_0$ is determined by the normalization condition: $$\pi_0 = \left(1 + \sum_{n=1}^{\infty} \prod_{k=1}^n \frac{\lambda_{k-1}}{\mu_k}\right)^{-1}$$ ## 3.5 Markov Jump Processes ### 3.5.1 Construction Method A Markov jump process can be constructed as follows: 1. **Waiting Time**: The waiting time $T_i$ in state $i$ follows an exponential distribution with parameter $\nu_i = -q_{ii}$ 2. **Jump Probability**: The probability of jumping from state $i$ to state $j$ is: $$p_{ij} = \frac{q_{ij}}{\nu_i}, \quad j \neq i$$ ### 3.5.2 Jump Times and Jump Chain Let $\{T_n\}$ be the jump time sequence and $\{Y_n\}$ be the jump chain, then: $$X(t) = Y_n, \quad \text{when } T_n \leq t < T_{n+1}$$ **Transition Matrix of Jump Chain**: $$P = (p_{ij})_{i,j \in S}$$ ```mermaid graph TB A[State i] -->|Waiting time ~ Exp(νᵢ)| B[Ready to Jump] B -->|Probability pᵢⱼ| C[State j] B -->|Probability pᵢₖ| D[State k] B -->|Probability pᵢₗ| E[State l] ``` ## 3.6 Python Implementation and Applications ```python import numpy as np import matplotlib.pyplot as plt from scipy.linalg import expm from scipy.stats import expon, poisson import seaborn as sns class ContinuousTimeMarkovChain: """Continuous-Time Markov Chain Class""" def __init__(self, generator_matrix, state_names=None): """ Initialize continuous-time Markov chain Parameters: ----------- generator_matrix : array-like Generator matrix Q state_names : list, optional State names """ self.Q = np.array(generator_matrix, dtype=float) self.n_states = self.Q.shape[0] if state_names is None: self.state_names = [f"State_{i}" for i in range(self.n_states)] else: self.state_names = state_names self._validate_generator_matrix() def _validate_generator_matrix(self): """Validate generator matrix""" if self.Q.shape[0] != self.Q.shape[1]: raise ValueError("Generator matrix must be square") # Check diagonal elements are negative, off-diagonal elements are non-negative for i in range(self.n_states): if self.Q[i, i] > 0: raise ValueError(f"Diagonal element Q[{i},{i}] must be non-positive") for j in range(self.n_states): if i != j and self.Q[i, j] < 0: raise ValueError(f"Off-diagonal element Q[{i},{j}] must be non-negative") # Check row sums are 0 row_sums = np.sum(self.Q, axis=1) if not np.allclose(row_sums, 0.0, atol=1e-10): raise ValueError("Each row of generator matrix must sum to 0") def transition_matrix(self, t): """Calculate transition matrix P(t) = exp(Qt) at time t""" return expm(self.Q * t) def stationary_distribution(self): """Calculate stationary distribution""" # Solve πQ = 0, i.e., Q^T π = 0 # Add normalization constraint A = np.vstack([self.Q.T, np.ones(self.n_states)]) b = np.append(np.zeros(self.n_states), 1) # Least squares solution pi, _, _, _ = np.linalg.lstsq(A, b, rcond=None) # Ensure non-negativity pi = np.abs(pi) pi = pi / np.sum(pi) return pi def simulate_jump_chain(self, initial_state, max_time, seed=None): """Simulate Markov jump process""" if seed is not None: np.random.seed(seed) # Calculate jump rates and jump probabilities rates = -np.diag(self.Q) P_jump = np.zeros_like(self.Q) for i in range(self.n_states): if rates[i] > 0: for j in range(self.n_states): if i != j: P_jump[i, j] = self.Q[i, j] / rates[i] # Simulate path times = [0.0] states = [initial_state] current_state = initial_state current_time = 0.0 while current_time < max_time: # Waiting time if rates[current_state] > 0: waiting_time = np.random.exponential(1.0 / rates[current_state]) else: # Absorbing state break current_time += waiting_time if current_time >= max_time: break # Choose next state if np.sum(P_jump[current_state, :]) > 0: next_state = np.random.choice( self.n_states, p=P_jump[current_state, :] ) else: break times.append(current_time) states.append(next_state) current_state = next_state return times, states def plot_sample_path(self, times, states, figsize=(12, 6)): """Plot sample path""" plt.figure(figsize=figsize) # Draw step function for i in range(len(times) - 1): plt.hlines(states[i], times[i], times[i+1], colors='blue', linewidth=2) if i < len(times) - 2: plt.vlines(times[i+1], states[i], states[i+1], colors='red', linestyles='dashed', alpha=0.7) # Last segment if len(times) > 1: plt.hlines(states[-1], times[-1], times[-1] + (times[-1] - times[-2]) * 0.1, colors='blue', linewidth=2) plt.xlabel('Time') plt.ylabel('State') plt.title('Continuous-Time Markov Chain Sample Path') plt.grid(True, alpha=0.3) # Set y-axis ticks to state names plt.yticks(range(self.n_states), self.state_names) plt.tight_layout() plt.show() class PoissonProcess: """Poisson Process Class""" def __init__(self, rate): """ Initialize Poisson process Parameters: ----------- rate : float Poisson process intensity parameter λ """ self.rate = rate def simulate(self, max_time, seed=None): """Simulate Poisson process""" if seed is not None: np.random.seed(seed) times = [0.0] current_time = 0.0 while current_time < max_time: # Waiting time for next event inter_arrival = np.random.exponential(1.0 / self.rate) current_time += inter_arrival if current_time < max_time: times.append(current_time) return times def count_process(self, times, max_time): """Generate counting process from event times""" # Create fine-grained time grid t_grid = np.linspace(0, max_time, int(max_time * 100) + 1) counts = np.zeros_like(t_grid) for i, t in enumerate(t_grid): counts[i] = np.sum(np.array(times) <= t) return t_grid, counts def plot_sample_path(self, times, max_time, figsize=(12, 6)): """Plot Poisson process sample path""" t_grid, counts = self.count_process(times, max_time) plt.figure(figsize=figsize) plt.step(t_grid, counts, where='post', linewidth=2, label=f'λ = {self.rate}') # Mark event times for t in times[1:]: # Skip time 0 plt.axvline(x=t, color='red', linestyle='--', alpha=0.5) plt.xlabel('Time') plt.ylabel('Event Count N(t)') plt.title(f'Poisson Process Sample Path (λ = {self.rate})') plt.legend() plt.grid(True, alpha=0.3) plt.tight_layout() plt.show() return t_grid, counts class BirthDeathProcess: """Birth-Death Process Class""" def __init__(self, birth_rates, death_rates, state_names=None): """ Initialize birth-death process Parameters: ----------- birth_rates : array-like Birth rate sequence [λ₀, λ₁, λ₂, ...] death_rates : array-like Death rate sequence [μ₀, μ₁, μ₂, ...], μ₀ is usually 0 """ self.birth_rates = np.array(birth_rates) self.death_rates = np.array(death_rates) self.n_states = len(birth_rates) if state_names is None: self.state_names = [f"State_{i}" for i in range(self.n_states)] else: self.state_names = state_names # Construct generator matrix self.Q = self._build_generator_matrix() def _build_generator_matrix(self): """Construct generator matrix""" Q = np.zeros((self.n_states, self.n_states)) for i in range(self.n_states): # Birth transition if i < self.n_states - 1: Q[i, i + 1] = self.birth_rates[i] # Death transition if i > 0: Q[i, i - 1] = self.death_rates[i] # Diagonal element Q[i, i] = -(Q[i, :].sum() - Q[i, i]) return Q def stationary_distribution(self): """Calculate stationary distribution""" pi = np.zeros(self.n_states) pi[0] = 1.0 # Recursive calculation for n in range(1, self.n_states): if self.death_rates[n] > 0: product = 1.0 for k in range(1, n + 1): product *= self.birth_rates[k - 1] / self.death_rates[k] pi[n] = product else: pi[n] = 0.0 # Normalization total = np.sum(pi) if total > 0: pi = pi / total else: # If infinite series diverges, return uniform distribution pi = np.ones(self.n_states) / self.n_states return pi # Application Example 1: M/M/1 Queue Model def mm1_queue_example(): """M/M/1 Queueing System Example""" print("M/M/1 Queueing System Analysis") print("=" * 30) # Parameters arrival_rate = 2.0 # Arrival rate λ service_rate = 3.0 # Service rate μ max_capacity = 10 # Maximum capacity # Construct birth-death process birth_rates = [arrival_rate] * max_capacity # Last state has no birth birth_rates[-1] = 0 # Capacity constraint death_rates = [0] + [service_rate] * (max_capacity - 1) queue = BirthDeathProcess(birth_rates, death_rates, [f"Queue_{i}" for i in range(max_capacity)]) # Calculate stationary distribution pi = queue.stationary_distribution() print("\nStationary distribution (probability of each queue length):") for i, prob in enumerate(pi): if prob > 1e-6: print(f"Queue length {i}: {prob:.4f}") # Performance metrics print("\nPerformance metrics:") mean_customers = np.sum(np.arange(max_capacity) * pi) prob_empty = pi[0] prob_full = pi[-1] print(f"Average number of customers in system: {mean_customers:.3f}") print(f"System idle probability: {prob_empty:.3f}") print(f"System full probability: {prob_full:.3f}") # Theoretical value (infinite capacity M/M/1) rho = arrival_rate / service_rate if rho < 1: theoretical_mean = rho / (1 - rho) print(f"Theoretical average customers (infinite capacity): {theoretical_mean:.3f}") return queue # Application Example 2: Telephone Exchange Model def telephone_exchange_example(): """Telephone exchange birth-death process model""" print("\nTelephone Exchange Model") print("=" * 25) # Parameters call_rate = 5.0 # Call arrival rate line_capacity = 8 # Number of lines service_rate = 1.0 # Service rate per line # Birth-death process parameters birth_rates = [call_rate] * (line_capacity + 1) birth_rates[-1] = 0 # No new calls when full death_rates = [0] + [i * service_rate for i in range(1, line_capacity + 1)] exchange = BirthDeathProcess(birth_rates, death_rates, [f"Busy_{i}" for i in range(line_capacity + 1)]) # Calculate stationary distribution pi = exchange.stationary_distribution() print("\nStationary distribution:") for i, prob in enumerate(pi): print(f"Busy lines {i}: {prob:.4f}") # Blocking probability (numerical result of Erlang B formula) blocking_prob = pi[-1] print(f"\nCall blocking probability: {blocking_prob:.4f}") # System utilization utilization = np.sum(np.arange(1, line_capacity + 1) * pi[1:]) / line_capacity print(f"System utilization: {utilization:.3f}") return exchange # Application Example 3: Population Model def population_model_example(): """Simple population growth model""" print("\nPopulation Growth Model") print("=" * 20) # Parameters max_population = 15 birth_rate_per_individual = 0.5 death_rate_per_individual = 0.3 # Birth-death process parameters (dependent on population size) birth_rates = [i * birth_rate_per_individual for i in range(max_population)] birth_rates.append(0) # Reached maximum capacity death_rates = [0] + [i * death_rate_per_individual for i in range(1, max_population + 1)] population = BirthDeathProcess(birth_rates, death_rates, [f"Population_{i}" for i in range(max_population + 1)]) # Calculate stationary distribution pi = population.stationary_distribution() print("\nPopulation distribution:") for i, prob in enumerate(pi): if prob > 1e-6: print(f"Population {i}: {prob:.4f}") # Expected population expected_population = np.sum(np.arange(max_population + 1) * pi) print(f"\nExpected population size: {expected_population:.2f}") # Extinction probability extinction_prob = pi[0] print(f"Population extinction probability: {extinction_prob:.4f}") return population # Poisson process example def poisson_process_example(): """Poisson process application example""" print("\nPoisson Process Applications") print("=" * 20) # Create Poisson processes with different intensities rates = [1.0, 2.0, 3.0] max_time = 10.0 plt.figure(figsize=(15, 5)) for i, rate in enumerate(rates): poisson_proc = PoissonProcess(rate) times = poisson_proc.simulate(max_time, seed=42 + i) plt.subplot(1, 3, i + 1) t_grid, counts = poisson_proc.count_process(times, max_time) plt.step(t_grid, counts, where='post', linewidth=2) plt.axhline(y=rate * max_time, color='red', linestyle='--', label=f'Expected value λt = {rate * max_time}') # Mark event times for t in times[1:]: plt.axvline(x=t, color='gray', linestyle=':', alpha=0.5) plt.xlabel('Time') plt.ylabel('Event Count N(t)') plt.title(f'Poisson Process (λ = {rate})') plt.legend() plt.grid(True, alpha=0.3) plt.tight_layout() plt.show() # Verify Poisson distribution print("\nVerify Poisson distribution (λ=2, t=5):") rate = 2.0 time_interval = 5.0 expected_count = rate * time_interval # Simulate multiple realizations n_simulations = 1000 counts = [] poisson_proc = PoissonProcess(rate) for i in range(n_simulations): times = poisson_proc.simulate(time_interval, seed=i) counts.append(len(times) - 1) # Subtract initial time 0 # Compare empirical and theoretical distributions empirical_mean = np.mean(counts) empirical_var = np.var(counts) print(f"Empirical mean: {empirical_mean:.3f}") print(f"Theoretical mean: {expected_count:.3f}") print(f"Empirical variance: {empirical_var:.3f}") print(f"Theoretical variance: {expected_count:.3f}") return counts # Run all examples if __name__ == "__main__": print("Chapter 3: Continuous-Time Markov Process Application Examples\n") # M/M/1 queue system mm1_queue = mm1_queue_example() # Telephone exchange exchange = telephone_exchange_example() # Population model population = population_model_example() # Poisson process poisson_counts = poisson_process_example() # Visualize generator matrices of birth-death processes try: plt.figure(figsize=(12, 4)) models = [ (mm1_queue, "M/M/1 Queue"), (exchange, "Telephone Exchange"), (population, "Population Model") ] for i, (model, title) in enumerate(models): plt.subplot(1, 3, i + 1) sns.heatmap(model.Q, annot=True, fmt='.2f', cmap='RdBu_r', center=0, square=True, cbar_kws={'shrink': 0.8}) plt.title(f"{title}\nGenerator Matrix") plt.xlabel("State") plt.ylabel("State") plt.tight_layout() plt.show() except ImportError: print("\nNeed to install seaborn library to display heatmaps") ``` ## 3.7 Stochastic Differential Equations and Diffusion Processes ### 3.7.1 Introduction to Diffusion Processes **Diffusion processes** are continuous-time, continuous-state Markov processes, typically described by stochastic differential equations (SDEs): $$dX(t) = \mu(X(t), t) dt + \sigma(X(t), t) dW(t)$$ where: - $\mu(x, t)$: Drift coefficient - $\sigma(x, t)$: Diffusion coefficient - $W(t)$: Brownian motion ### 3.7.2 Geometric Brownian Motion **Geometric Brownian motion** is an important diffusion process in financial modeling: $$dS(t) = \mu S(t) dt + \sigma S(t) dW(t)$$ **Solution**: $$S(t) = S(0) \exp\left(\left(\mu - \frac{\sigma^2}{2}\right)t + \sigma W(t)\right)$$ This is the foundation of the Black-Scholes model. ## 3.8 Chapter Summary This chapter covered the core theory of continuous-time Markov processes: **Main Concepts**: 1. **Continuous-Time Markov Property**: Memorylessness in a continuous time framework 2. **Generator Matrix**: Matrix describing instantaneous transition rates 3. **Kolmogorov Equations**: Differential equations describing the evolution of transition probabilities 4. **Poisson Process**: The most important counting process 5. **Birth-Death Process**: Markov process with special structure **Key Formulas**: - Transition probability: $P(t) = e^{Qt}$ - Kolmogorov equation: $\frac{d}{dt}P(t) = QP(t) = P(t)Q$ - Poisson process: $P(N(t) = n) = \frac{(\lambda t)^n e^{-\lambda t}}{n!}$ - Stationary distribution: $\pi Q = 0$ **Practical Applications**: 1. **Queueing Theory**: M/M/1 queueing system modeling 2. **Communication Systems**: Telephone exchange capacity analysis 3. **Population Dynamics**: Biological population growth models 4. **Financial Risk**: Default time and credit risk modeling Continuous-time Markov processes provide an important theoretical foundation for interest rate models, credit risk models, and other applications in finance.