stochastic-processes

# Stochastic Processes — Theory and Applications Covers: **Probability Foundations · Markov Chains · Poisson Processes · Brownian Motion · Stochastic Calculus · Applications** ----- ## Probability Foundations ### Probability Spaces and Random Variables A probability space consists of three components: the sample space Omega representing all possible outcomes, a sigma-algebra F containing all measurable events, and a probability measure P assigning probabilities to events. Understanding this foundation is essential for rigorously defining stochastic processes. A random variable is a measurable function X: Omega -> R that assigns a real number to each outcome. The distribution of X is characterized by its cumulative distribution function F(x) = P(X <= x), or by its probability density function f(x) for continuous variables. **Key Distributions:** | Distribution | Parameters | Mean | Variance | |--------------|-----------|------|----------| | Normal | mu, sigma^2 | mu | sigma^2 | | Exponential | lambda | 1/lambda | 1/lambda^2 | | Poisson | lambda | lambda | lambda | | Uniform | a, b | (a+b)/2 | (b-a)^2/12 | | Bernoulli | p | p | p(1-p) | ### Expectations and Moments ```python import numpy as np from scipy import stats # Calculate expectations def expectation(distribution, func=lambda x: x): """Calculate E[g(X)] for distribution""" if hasattr(distribution, 'expect'): return distribution.expect(func) # Monte Carlo approximation samples = distribution.rvs(100000) return np.mean(func(samples)) # Normal distribution example normal = stats.norm(loc=5, scale=2) E_x = normal.mean() # 5 Var_x = normal.var() # 4 E_x2 = normal.moment(2) # E[X^2] = Var + E[X]^2 = 4 + 25 = 29 # Conditional expectation # E[X|Y=y] - expectation of X given Y=y def conditional_expectation(joint_samples, x_idx, y_idx, y_value): """Calculate E[X|Y=y] from samples""" mask = joint_samples[:, y_idx] == y_value return np.mean(joint_samples[mask, x_idx]) # Moment generating function def mgf_normal(t, mu, sigma): """MGF of normal distribution""" return np.exp(mu * t + 0.5 * sigma**2 * t**2) # Characteristic function def cf_normal(t, mu, sigma): """Characteristic function of normal""" return np.exp(1j * mu * t - 0.5 * sigma**2 * t**2) ``` ----- ## Markov Chains ### Definition and Properties A stochastic process {X_n} is a Markov chain if it satisfies the Markov property: given the present, the future is independent of the past. Formally, P(X_{n+1} = j | X_n = i, X_{n-1} = i_{n-1}, ..., X_0 = i_0) = P(X_{n+1} = j | X_n = i). The transition matrix P with entries P_{ij} = P(X_{n+1} = j | X_n = i) characterizes the chain. A Markov chain is homogeneous if these transition probabilities are independent of n. **Key Properties:** - **Irreducible** — Every state can be reached from every other state - **Recurrent** — Expected return time to any state is finite - **Transient** — Some states may never be revisited - **Periodic** — Returns to states occur at regular intervals - **Aperiodic** — No regular return pattern ### Markov Chain Implementation ```python import numpy as np from dataclasses import dataclass from typing import Dict, List, Optional import matplotlib.pyplot as plt @dataclass class MarkovChain: """Homogeneous Markov chain""" states: List[str] transition_matrix: np.ndarray initial_distribution: Optional[np.ndarray] = None def __post_init__(self): n = len(self.states) if self.initial_distribution is None: self.initial_distribution = np.ones(n) / n # Validate assert np.allclose(self.transition_matrix.sum(axis=1), 1) assert len(self.states) == self.transition_matrix.shape[0] def simulate(self, n_steps: int, seed: int = None) -> List[str]: """Simulate chain for n_steps""" if seed is not None: np.random.seed(seed) path = [] current = np.random.choice(len(self.states), p=self.initial_distribution) path.append(self.states[current]) for _ in range(n_steps - 1): current = np.random.choice( len(self.states), p=self.transition_matrix[current] ) path.append(self.states[current]) return path def transition_prob(self, from_state: str, to_state: str, steps: int = 1) -> float: """Probability of transitioning from one state to another in n steps""" i = self.states.index(from_state) j = self.states.index(to_state) if steps == 1: return self.transition_matrix[i, j] return (self.transition_matrix ** steps)[i, j] def stationary_distribution(self) -> np.ndarray: """Find stationary distribution pi such that pi = pi * P""" eigenvalues, eigenvectors = np.linalg.eig( self.transition_matrix.T ) # Find eigenvector with eigenvalue closest to 1 idx = np.argmin(np.abs(eigenvalues - 1)) stationary = np.real(eigenvectors[:, idx]) return stationary / stationary.sum() def verify_ergodicity(self) -> Dict: """Check if chain is ergodic""" n = len(self.states) # Check irreducibility reachable = np.linalg.matrix_power( self.transition_matrix + np.eye(n), n - 1 ) > 0 irreducible = np.all(reachable) # Check aperiodicity eigenvalues = np.linalg.eigvals(self.transition_matrix) eigenvalues = np.round(eigenvalues, 10) eigenvalues = eigenvalues[np.abs(eigenvalues) > 0.01] is_aperiodic = not any( np.isclose(np.abs(e), 1) and np.isclose(np.abs(e ** n - 1), 0) for n in range(2, n + 1) for e in eigenvalues ) stationary = self.stationary_distribution() return { "irreducible": irreducible, "aperiodic": is_aperiodic, "has_stationary": True, "stationary_distribution": dict(zip(self.states, stationary)) } # Example: Weather model weather_states = ["sunny", "cloudy", "rainy"] weather_transitions = np.array([ [0.7, 0.2, 0.1], # Sunny -> sunny, cloudy, rainy [0.3, 0.4, 0.3], # Cloudy -> sunny, cloudy, rainy [0.2, 0.3, 0.5] # Rainy -> sunny, cloudy, rainy ]) weather_chain = MarkovChain(weather_states, weather_transitions) # Simulate path = weather_chain.simulate(30) print("Weather sequence:", " -> ".join(path[:10])) # Stationary distribution stationary = weather_chain.stationary_distribution() print("\nStationary distribution:") for state, prob in zip(weather_states, stationary): print(f" {state}: {prob:.3f}") # Verify ergodicity print("\nErgodicity:", weather_chain.verify_ergodicity()) ``` ### Continuous-Time Markov Chains ```python class ContinuousTimeMarkovChain: """Continuous-time Markov chain with rate matrix Q""" def __init__(self, states: List[str], rate_matrix: np.ndarray): self.states = states self.rate_matrix = rate_matrix # Verify row sums to 0 assert np.allclose(rate_matrix.sum(axis=1), 0) # Transition matrix for embedded chain n = len(states) self.embedded_matrix = np.zeros((n, n)) for i in range(n): if rate_matrix[i].sum() > 0: self.embedded_matrix[i] = rate_matrix[i] / -rate_matrix[i, i] def simulate_jump_times(self, n_jumps: int, initial_state: int = 0, seed: int = None): """Simulate path with jump times""" if seed is not None: np.random.seed(seed) path = [initial_state] times = [0.0] current = initial_state for _ in range(n_jumps): # Time to next jump (exponential) rate = -self.rate_matrix[current, current] dt = np.random.exponential(1 / rate) times.append(times[-1] + dt) # Next state (embedded chain) current = np.random.choice( len(self.states), p=self.embedded_matrix[current] ) path.append(current) return { "states": [self.states[s] for s in path], "times": times, "path": path } ``` ----- ## Poisson Processes ### Definition and Properties A Poisson process with rate lambda counts events occurring randomly in time. It satisfies: the number of events in any interval follows a Poisson distribution, events occur independently of past events, and the process has stationary increments. The interarrival times (times between consecutive events) are independently and exponentially distributed with mean 1/lambda. This property allows simple simulation of Poisson processes. ```python class PoissonProcess: """Poisson process with rate lambda""" def __init__(self, lambda_rate: float): self.lambda_rate = lambda_rate def simulate(self, T: float, seed: int = None) -> Dict: """Simulate Poisson process up to time T""" if seed is not None: np.random.seed(seed) # Generate interarrival times interarrivals = [] t = 0 while True: dt = np.random.exponential(1 / self.lambda_rate) t += dt if t > T: break interarrivals.append(dt) # Arrival times arrival_times = np.cumsum(interarrivals) counts = np.arange(1, len(arrival_times) + 1) return { "arrival_times": arrival_times, "counts": counts, "n_events": len(arrival_times) } def probability_k_events(self, k: int, T: float) -> float: """P(N(T) = k)""" from scipy.special import gammainc mu = self.lambda_rate * T return np.exp(-mu) * (mu ** k) / np.math.factorial(k) def arrival_time_distribution(self, n: int) -> np.ndarray: """Distribution of arrival times given N(T) = n""" # Order statistics of uniform return np.sort(np.random.uniform(0, 1, n)) # Compound Poisson process class CompoundPoissonProcess: """Compound Poisson process: N(t) with i.i.d. jump sizes""" def __init__(self, lambda_rate: float, jump_distribution): self.lambda_rate = lambda_rate self.jump_distribution = jump_distribution def simulate(self, T: float, seed: int = None) -> Dict: """Simulate compound Poisson process""" if seed is not None: np.random.seed(seed) # Generate Poisson arrivals poisson = PoissonProcess(self.lambda_rate) arrivals = poisson.simulate(T) # Generate jump sizes n = arrivals["n_events"] jump_sizes = self.jump_distribution.rvs(n) # Compound process values values = np.cumsum(jump_sizes) times = arrivals["arrival_times"] return { "times": times, "values": values, "jump_sizes": jump_sizes, "final_value": values[-1] if n > 0 else 0 } ``` ----- ## Brownian Motion ### Definition and Properties Standard Brownian motion B(t) is a stochastic process with: B(0) = 0, stationary increments, independent increments, and continuous paths almost surely. The increments B(t) - B(s) ~ N(0, t-s). Brownian motion is a martingale, has infinite variation, is nowhere differentiable, and exhibits fractal behavior. These properties make it both mathematically interesting and practically useful for modeling random fluctuations. ```python class BrownianMotion: """Standard Brownian motion""" def __init__(self, drift: float = 0, diffusion: float = 1): self.drift = drift self.diffusion = diffusion def simulate(self, T: float, n_steps: int, seed: int = None) -> Dict: """Simulate Brownian motion""" if seed is not None: np.random.seed(seed) dt = T / n_steps t = np.linspace(0, T, n_steps + 1) # Increments dW = np.random.normal(0, np.sqrt(dt), n_steps) # Path W = np.zeros(n_steps + 1) W[1:] = np.cumsum(dW) # Add drift and diffusion if self.drift != 0 or self.diffusion != 1: W = self.drift * t + self.diffusion * W return { "time": t, "values": W, "final_value": W[-1] } def simulate_bridge(self, T: float, n_steps: int, start: float, end: float, seed: int = None): """Brownian bridge from start to end""" if seed is not None: np.random.seed(seed) dt = T / n_steps t = np.linspace(0, T, n_steps + 1) # Standard Brownian bridge W = self.simulate(T, n_steps).values # Transform to bridge