1.2 Distribution of the Poisson process

1.3 Poisson process and exponential random variables

Poisson processes and exponential (and gamma) random variables are intimately related. Given a Poisson process we consider the interarrival times T_1, T_2, \cdots: T_1 is time of the occurence of the first event, T_2 is the time elapsed between the first and second event, and so on.

1.4 Poisson process and uniform distribution

Let us start with a special case and assume assume that N_t=1, that is exactly one event has occured in [0,t]. Since the Poisson process has independent increments it seems reasonable the event may have occur with equal probability at any time on [0,t]. Indeed for s\le t we have, using the independence of increments. \begin{aligned} P\{T_1 \le s | N_t=1\} & = \frac{P\{T_1 \le s , N_t=1\}}{P\{N_t=1\}} = \frac{P\{\textrm{1 event in }(0,s] \textrm{ and no event in (s,t]} \}}{P\{N_t=1\}} \\ &= \frac{ \lambda s e^{-\lambda s} e^{-\lambda(t-s)} }{\lambda t e^{-\lambda t}} = \frac{s}{t} \end{aligned} and thus the density of the arrival time T_1 conditioned on N_t=1 is uniform on [0,t]

We study further the properties of the arrival times S_1< S_2 <\cdots of a Poisson process. The following result tells us that they follow a uniform distribution on [0,t].

1.5 Simulation of Poisson process

The characterization of the Poisson process in terms of exponential random variables suggest immediately a very simple algorithm to simulate N_t.

#| code-fold: true
#| code-summary: "Show the code"


#| label: fig-Poisson
#| fig-cap: "A sample of the Poisson process"

import numpy as np
import matplotlib.pyplot as plt

def poisson_process(rate, T):
    # generate a Poisson process with rate `rate` up to time T
    t = 0
    events = [0]
    while t < T:
        # draw time until next event from an exponential distribution
        dt = np.random.exponential(1/rate)
        t += dt
        if t < T:
            # record time of event
            events.append(t)
    return events

rate=10
T=200
events=np.asarray(poisson_process(rate,T))

plt.figure() 
plt.step(events,np.arange(0,len(events)), label=r'sample of $N_t$' )
plt.plot(events,rate*events, label=r'$E[N_t]=\lambda t$')
plt.legend()
plt.xlabel('t')
plt.ylabel(r'$ N_t $')
plt.show()

1.6 Simulation of a Poisson random variable

• It is not easy to simulate directly a Poisson random variable X from its pdf/cdf but we can do it elegantly using its relation with exponential random variable. To do this generate independent exponential random variable until they sum up to 1 (so as to generate X=N_1) and use the relation between exponential and uniform.
  \begin{aligned} X =n & \iff \sum_{k=1}^n T_k \le 1 < \sum_{k=1}^{n+1} T_k \iff \sum_{k=1}^n - \frac{1}{\lambda} \ln(U_k) \le 1 < \sum_{k=1}^{n+1} - \frac{1}{\lambda} \ln(U_k) \\ & \iff \ln( \prod_{k=1}^n U_k) \ge - \lambda > \ln( \prod_{k=1}^{n+1} U_k) \iff \prod_{k=1}^n U_k \ge e^{-\lambda} > \prod_{k=1}^{n+1} U_k \end{aligned}

• Algorithm to simulate a Poisson random variable with parameter \lambda: Generate random numbers until their product is smaller than e^{-\lambda}.

  • Generate random number U_1, U_2, \cdots.

  • Set X=n if n+1 = \inf \left\{ j \,:\, \prod_{k=1}^{j} U_j < e^{-\lambda} \right\}

1.7 Long time behavior of the Poisson process

• We investigate the behavior of N_t for large time. We prove a CLT type result, namely that
  \lim_{t \to \infty} \frac{N_t-\lambda t}{\sqrt{\lambda t}} = Z \textrm{ in distribution} where Z is a standard normal RV.

• Recall that the characteristic function of a Poisson RV Y with parameter \mu is E[e^{i\alpha Y}]= e^ {\mu (e^{i\alpha} -1)}. Therefore E\left[ e^{i \alpha \frac{N_t - \lambda t}{\sqrt{\lambda t}}} \right] = e^{\lambda t \left( e^{i \frac{\alpha}{\sqrt{\lambda t}}} -i \frac{\alpha}{\sqrt{\lambda t}} -1\right)} Expanding the exponential we have \lambda t \left(e^{i \frac{\alpha}{\sqrt{\lambda t}}} -i \frac{\alpha}{\sqrt{\lambda t}} -1\right) = - \frac{\alpha^2}{2} + O( \frac{1}{\sqrt{\lambda t}}) and thus \lim_{t \to \infty} E\left[ e^{i \alpha \frac{N_t - \lambda t}{\sqrt{\lambda t}}} \right]=e^{-\frac{\alpha^2}{2}}.

• The same computation shows also that, for any fixed t, \lim_{\lambda \to \infty} \frac{N_t-\lambda t}{\sqrt{\lambda t}} = Z \textrm{ in distribution} since rescaling the parameter is equivalent to rescaling time.

1.8 Sampling a Poisson process

We can sample or split a Poisson process. Suppose that every event of a Poisson process (independently of the other events) comes into two different types, say type 1 with probability p and type 2 with probability q=1-p.

1.9 The coupon collecting problem

1.10 Poisson process with variable rate

1.11 Queueing model with infinitely many servers

• Assume that the the flow of customers entering an onine store follows a Poisson process N_t with rate \lambda. The time S spent in the store for a single customer (browsing around, checking out, etc..) is given by tis CDF G(t)=P\{ S \le t\} and we assume that the customers are independent of each other.

• To figure out how to allocate ressources one wants to figure out what is number of customers, M_t, which are still in the sytem at time t.

• To find the distribution of M_t let us consider one of the customer by time t. If he arrived at time s\le t then he will have left the system at time t with probability G(t-s) and will still be in the system by time t with probability 1-G(t-s). Since the arrival time of that customer is uniform on [0,t] the distribution of M_t is Poisson with mean E[M_t] = \int_0^t (1- G(t-s))ds = \lambda \int_0^t (1 - G(s))ds,.
  For large t, we see that E[M_t] \approx \lambda E[S].

1.12 Compound Poisson process

2.2 Definition of a continuous time Markov chain

2.3 Solving the Komogorov equation

For finite state space S finding the transition probability reduces to a linear algebra problem. For example if the matrix A is diagonalizable (e.g. if the matrix symmetric or if all the eigenvalues are distinct). let us denote the eigenvalues of A are \lambda_1, \lambda_2, \lambda_N of A with corresponding the eigenvectors f_1, \cdots, f_N. Note that 0 is always an eigenvalue with eigenvector \begin{pmatrix} 1 & 1 & \cdots &1 \end{pmatrix}^T (the sum of the rows is equal to 0).

Then we have D = Q^{-1} A Q = \begin{pmatrix} \lambda_1 & 0 & \cdots & 0 \\ 0 & \lambda_2 & 0 & \cdots \\ \cdots &\cdots & \ddots & \cdots \\ \cdots & \cdots & 0 & \lambda_N \end{pmatrix} We find then that P_t = e^{At} = Q E^{Dt} Q^{-1} = Q \begin{pmatrix} e^{\lambda_1 t} & 0 & \cdots & 0 \\ 0 & e^{\lambda_2 t} & 0 & \cdots \\ \cdots &\cdots & \ddots & \cdots \\ \cdots & \cdots & 0 & e^{\lambda_N t} \end{pmatrix} Q^{-1}

2.4 Stochastic simulation algorithm

In this section we propose an alternative descritpion which is more probabilistic in nature and allows us to construct the paths of the Markov chains.
To any pair of states (i,j) we associate a “clock” T(i,j) which is an exponential random variable with paramter (rate) \alpha(i,j). All the random variables used are assumed to be independent.

2.5 Stochastic Simulation algorithm (or Gillepsie algorithm)

The description of Markov chain using the rates is essentially a pseudo-code to generat trajectories of X_t.

#| code-fold: true
#| code-summary: "Show the code"

import numpy as np
import matplotlib.pyplot as plt

# SSA function with A matrix
def stochastic_simulation_algorithm(A, initial_state, max_time):
    num_states = A.shape[0]
    current_state = initial_state
    current_time = 0.0
    times = [current_time]
    states = [current_state]

    while current_time < max_time:
        rates = A[current_state, :].copy()
        rates[current_state] = 0
        total_rate = rates.sum()

        if total_rate <= 0:
            break

        delta_t = np.random.exponential(1 / total_rate)
        current_time += delta_t
        if current_time > max_time:
            break

        probabilities = rates / total_rate
        next_state = np.random.choice(num_states, p=probabilities)

        times.append(current_time)
        states.append(next_state)
        current_state = next_state

    return times, states

# Plot function
def plot_trajectory(times, states):
    plt.step(times, states, where='post', linewidth=2)
    plt.xlabel('Time')
    plt.ylabel('State')
    plt.title('Markov Chain Trajectory (SSA)')
    plt.grid(True)
    plt.tight_layout()
    plt.show()

# Example rate matrix A
A=np.array([
    [-3, 3, 0, 0],
    [1, -4, 3, 0],
    [0, 2, -5, 3],
    [0, 0, 2, -2]
])

initial_state = 0
max_time = 30.0

# Run simulation
times, states = stochastic_simulation_algorithm(A, initial_state, max_time)

# Plot the trajectory
plot_trajectory(times, states)

2.6 Example: Uniformizable chain

• Consider a Markov chain Y_n with transition matrix Q and we assume that Q(i,i)=0. Then we pick rate \alpha(i)=1 for all states i. The times at which the Markov chain has transition is thus a sum of IID exponential, that at time t is then described by a Poisson process N_t. In other terms we have X_t = Y_{N_t} where N_t is a Poisson process with rate 1.

• In this case we can compute the transition matrix quite explicitly: \begin{aligned} P\{X_t=j|X_0=i\}& = P\{ Y_{N_t}=j| X_0=i\} \\ &= \sum_{n=0}^\infty P\{ Y_{N_t}=j, N_t=n| X_0=i\} \\ &= \sum_{n=0}^\infty P\{ Y_n=j| X_0=i\} P\{N_t=n\} \\ &= \sum_{n=0}^\infty \frac{e^{- t} t^n}{n!} Q^n(i,j) = e^{t(Q-I)} \end{aligned}

and the generator is given by A = (Q-I).

3.1 Stationary distributions and detailed balance

• A probability vector \pi is a stationary distribution for the Markov chain with generator A if
  \pi P_t = \pi \quad \textrm{ for all } t > 0 0 = \frac{d}{dt}\pi P_t = \pi A P_t \implies \pi A =0.

• In terms of the rate \alpha(i,j) we see that \pi is stationary if and only if \sum_{i \not= j} \pi(i) \alpha(i,j) = \pi(j)\alpha(j) = \sum_{i \not=j} \pi(j) \alpha(j,i) which we can interpret as balance equation. The quantity \pi(i)\alpha(i,j) is the rate at which the chain in a state \pi changes from i to j and the stationarity equation means that \text{ flow of probability away from state } i = \text{ flow of probability into state } i \quad \quad \text{ holds for all states} i

• As for discrete time we say that a Markov chain satisfies detailed balance if \pi(i) \alpha(i,j) = \pi(j) \alpha(j,i) \textrm{ for all } i \not= j and clearly detailed balance implies stationarity.

3.2 Convergence to the stationary distribution

If the state space is finite and the chain is irreducible then we have convergence to equilibrium.

3.3 Transient behavior

To study the transient behavior in continuous time we can use similar ideas as in discrete time.

• Absorption probabilities: the absorption probabilities do not depend on the time spent in every state so they can be computed using the transition matrix Q(i,j) for the embedded chain Y_n and the formula in Section [Absorption probabilities]

• Expected hittiing time: For example we have the following result

3.4 Explosion

• For a continuous time Markov chain X_t let us consider the time of the successive jumps S_1=T_1, S_2=T_1+T_2, S_3=T_1+T_2+T_3, \cdots Here the T_i are independent exponential but, in general, no identically distributed with parameters \alpha_i which depends on the state being visited. We have then E[S_n]= \sum_{i=1}^n \frac{1}{\lambda_i}

• Explosion: To see what can happen consider a Markov chain with rate \alpha(i,i+1)=(n+1)^2 and all other rates equal to 0. Then the Markov chain moves up by 1 at every jump like a Poisson process but at accelerated pace. We have then E[S_\infty]= \sum_{n=1}\frac{1}{n^2} < \infty so S_\infty < \infty with probability 1. So there are infinitely many jumps in finite time and X_t=+\infty after a finite time. This is called explosion.

3.5 Transience, recurrence, and positive recurrence.

• The recurrence and transience for a continuous time Markov chain can be defined as for discrete time. A an irreducbile Markov chain X_t is recurrent if the chain starting from a state i returns to i with probability 1 and transient if X_t never return to i with some non-zero probability.

• Transience/recurrence has nothing to do with the actual time spent in every state and so we have immediately X_t \textrm { with rates } \alpha(i,j) \textrm{ is } \begin{array}{l} \textrm{transient} \\ \textrm{recurrent}\end{array} \iff Y_n \textrm{ with transition } Q(i,j)=\frac{\alpha(i,j)}{\alpha(i)} \textrm{ is } \begin{array}{l} \textrm{transient} \\ \textrm{recurrent}\end{array}

• For positive recurrence we need to define the first return to a state i. This is most easily defined in terms of the return for the embeded Markov chain Y_n and its first return time \tau(i). The time until X_0 returns to i will the sum of the time spent in each state. If S_n denotes the time of the jumps and \tau(i)=n after visiting the state j_0, j_1, \cdots j_{n-1}, j_n=i then the return time is S_n= \sum_{k=0}^{n-1}(S_{k+1}- S_k) \quad \quad \text{ with }S_0=0 where S_{k+1}-S_k is exponential with parameter \alpha(j_k). The first return time to state i for the Markov chain X_t is then given by
  \Sigma(i) = \sum_{k=0}^{\tau(i)-1}(S_{k+1}- S_k)

3.6 Stationary distribution for recurrent chains

3.7 Ergodic theorem for positive recurrent Markov chains

We have the following theorem which is the exact counterpart of the discrete time case (and is proved very similarly so we will omit the proof).

3.8 Exercises

4.1 Birth and death process

• A general birth and death process is a continuous time Markov chain with state space \{0,1,2,\cdots\} and whose only non-zero transition rates are \begin{aligned} &\lambda(n) = \alpha(n,n+1) & = \textrm{ birth rate for a population of size } n & (\textrm{ for } n\ge 0) \\ &\mu(n)=\alpha(n,n-1) & = \textrm{ death rate for a population of size } n & (\textrm{ for } n\ge 1) \end{aligned}

• The Kolmogorov equations for the distribution of X_t are, for n\ge 1 \frac{d}{dt}p_t(n) = \underbrace{ \mu(n+1) p_t(n+1)}_{\textrm{ increase due do death} \atop \textrm{in a population of size } n+1} + \underbrace{\lambda_{n-1} p_t(n-1)}_{\textrm{ increase due do birth} \atop \textrm{in a population of size } n-1} - \underbrace{(\lambda(n) + \mu(n))p_t(n)}_{\textrm{ decrease due do birth/death} \atop \textrm{in a population of size} n}. For n=0, the eqation reads \displaystyle \frac{d}{dt}p_t(0)= \mu(1) p_t(1) - \lambda(0) p_t(0).

• The generator has the form \begin{matrix} 0 \\ 1 \\ 2\\ \vdots \end{matrix} \begin{pmatrix} -\lambda(0) & \lambda(0) & 0 & 0 & \dots \\ \mu(1) &-\lambda(1) -\mu(1) & \lambda(1) & 0 & \dots \\ 0 & \mu(2) &-\lambda(2) -\mu(2) & \lambda(2) & \dots \\ \vdots & \ddots & \ddots & \ddots & \\ \end{pmatrix}

4.2 Examples

• Poisson process: There is no death and the birth rate is constant \mu(n)=0 \quad \textrm{ and } \quad \lambda(n)=\lambda

• Queueing models: Imagine a Poisson process (with rate \lambda) describing the arrival of customers. In the systems there are k service station at which the service time is exponential with rate \mu. If one service station is empty, upon arrival a customer immediately enter a service station and stay in the system until being served. If all stations are currently busy serving other customers, the customer wait until one stay becomes free. These models are called M/M/k queues where M statnds for Markovian and k is the number of queues, the first M describing the arrival process and the second M the service time (which is exponential).
  

• M/M/1 queue: If there is only one service station \lambda(n)=\lambda \quad \textrm{ and } \quad \mu(n) =\mu

• M/M/k queue: If there is k service station then the rates are \lambda(n)=\lambda \quad \textrm{ and } \quad \mu(n) = \left\{ \begin{array}{cl} n \mu & n \le k \\ k \mu & n > k \end{array} \right.

• M/M/\infty queue: If there is k service station then the rates are \lambda(n)=\lambda \quad \textrm{ and } \quad \mu(n) =n\mu

4.3 Transience of birth/death chains

To study transience we use the embedded Markov chain and ?@thm-criteriontransience. Choosing the reference state 0 we look for a solution a(n) with a(0)=1 and 0 < a(n)< 1 for n\ge 1 of the system of equations \lambda(n) a(n+1) + \mu(n) a(n-1) = (\lambda(n) + \mu(n)) a(n) This leads to a(n+1) - a(n) = \frac{\mu(n)}{\lambda(n)} (a(n) - a(n-1)) = \cdots = \prod_{j=1}^n\frac{\mu(j)}{\lambda(j)} (a(1) - 1) and thus, by telescoping, a(n+1)= 1 + \sum_{k=0}^{n} (a(k+1) -a(k)) = 1 + \left[1+ \sum_{k=1}^{n} \prod_{j=1}^k\frac{\mu(j)}{\lambda(j)}\right] (a(1) - 1) Taking n \to \infty we must have a(n)\to 0 and \sum_{k=1}^\infty \prod_{j=1}^k\frac{\mu(j)}{\lambda(j)} < \infty and we find the solution a(n)=\frac{\sum_{k=n}^\infty \prod_{j=1}^k\frac{\mu(j)}{\lambda(j)}}{1+ \sum_{k=1}^\infty \prod_{j=1}^k\frac{\mu(j)}{\lambda(j)}}

4.4 Positive recurrence of birth/death chains

To study postive recurrence we simply solve for the stationary distribution \pi A=0. Since the embedded Markov chain satisfies detailed balance it is natural tro try to solve the detailed balance equations which amounts \pi(n)\lambda(n) = \pi(n+1)\mu(n+1) which is eaisly solved to find \pi(n) = \prod_{j=1}^n \frac{\lambda(j-1)}{\mu(j)} \pi(0) and thus we have a stationary distribution if and only if \sum_{k=1}^\infty \prod_{j=1}^k \frac{\lambda(j-1)}{\mu(j)} <\infty\,.

4.5 Speed of convergence via coupling M/M/\infty queue

We can analyze the M/M/\infty queue in more details by the following argument. Suppose X_0=j is the size of the queue at time o. Then at time t the people in queue are of two different types: either they were in queue at time 0 and are still in queue at time t or they arrive after time between times 0 and time t and are still in queue.

The probability that a person in queue at time 0 is still in queue at time t is e^{-\mu t}.

Since the arrival of a Poisson process condiioned on N_t are uniformly distributed on [0,t], the probability that someone not present at time 0 is still present at time t is q(t)=\frac{1}{t}\int_0^t (1- e^{-\mu(t-s)})= \frac{1- e^{-\mu t}}{\mu t} ,. By combining this we see that X_t = N_t + Y_t where N_t \text{ has a binomial distribution with parameters } j \text{ and } e^{-\mu t} Y_t \text{ has a Poisson distribution with parameters } \lambda q(t)t = \frac{\lambda}{\mu}(1 - e^{-\mu t})

4.6 Queueing networks

• A queueing network is an interconnected netwrok of service facilities called nodes. Each node has its own queueing rules and there are probabilistic rules to move between nodes and exit the system. Queueing netwroks can be either closed or open.

• A closed network has a fixed number K of nodes and a fixed number M of customers. No agents enters or exit the systems. Imagine for example a company with M trucks which can be either in a repair shop or assigned to a variety of tasks.

• A open network has a fixed number K of nodes but a variables number of customers. Agents can enters or exit the system at some of the nodes in the systems. Imagine for example data packets being routed in a computer network.

• To describe the process in a network with M nodes we consider the queue length process {\bf Q_t} = \left(Q_{1,t}, \cdots, Q_{K,t}\right) \quad \text{ where } Q_{i,t} = \text{ number of customers at node } i \text{ at time} t.

• The state space of the process for the queue length is process is given by \begin{aligned} S_M & = \left\{ {\bf n} = \left(n_1, \cdots, n_K\right) \,:\, n_i \ge 0\,,\, \sum_{i=1}^K n_i=M \right\} \quad \quad \text{ closed networks} \\ S & = \left\{ {\bf n} = \left(n_1, \cdots, n_K\right) \,:\, n_i \ge 0\,\right\}\quad \quad \quad \quad \quad \quad\quad\quad\quad \text{ open networks} \\ \end{aligned}

:::

4.7 Closed network of queues

To describe a closed queuing network we need to describing how the M agents move along the K nodes of the networks. This is described in terms of routing matrix P with entries P_{il} = \text{ Probability that an agent exiting node $i$ moves to node $l$ } It is possible that P_{ii}>0 describing the case where an agent may need two services. The matrix P is a stochastic matrix with P_{il}\ge 0 \quad \text{ and } \sum_{l=1}^K P_{il} =1 We assume that P is irreducible and we will denote by \pi the corresponding stationary distribution, i.e. \pi P=\pi.

We assume that at node i there is a single service station and that when at node i an agent is served with a rate b_i so that service times at nodes i are IID exponentials with parameters b_i.

Traffic equation: If the system is in equilbrium then the rate r_j at which customers leave a node j should be the to the rate at which customer leaves a node j. This leads to the equation \underbrace{r_j}_{\text{exit rate}} = \underbrace{\sum_{i} r_i P_{ij}}_{\text{entrance rate}} and therefore the stationary distribution should depend on the stationary distribution \pi for P.

4.8 Open queueing (Jackson) networks

In an open queuing network the number of customers in the system can be arbitrary and customers can both enter and leave the systems (both must occur to obtain a stable system.) We need the following ingredients

• Arrival rates a_i \ge 0 at each node which describing customers which enter the system at node i according to a Poisson process with rate a_i.

• Waiting times at each node with rate b_i>0 which describe the service time at node i.

• Routing matrix Q which describe the transition between nodes and the exits from the system. We assume Q is sub-stochastic, i.e. Q_{ij} \ge 0 and \sum_{j} Q_{i,j}\le 1 with at least one row have a sum \sum_{j} Q_{i,j} <1. We interpet Q_{i,j} = P( \text{ go to node } j \text{ after service at node } i ) q_i = 1-\sum_{j}Q(i,j) = P( \text{ exit the system } \text{ after service at node } i ) We can think of Q as describing the transition probablities of a set of transient states in a Markov chain with an absorbing state which corresponds to being out of the system. In particular we know that I-Q is invertible.