Sets up a chain on a countable state space through its transition matrix, derives the Chapman-Kolmogorov equations that make the n-step probabilities matrix powers, classifies states by recurrence, mean recurrence time and period, and characterises the stationary distribution as the left eigenvector of the transition matrix with eigenvalue one.
A random process is in general an unmanageable object: the joint law of X0,X1,…,Xn needs a number for every path. A Markov chain is the case in which that joint law collapses into a single matrix of one-step probabilities, because the chain's memory of its past is entirely contained in where it is now. Everything in this unit is a consequence of that collapse: n-step probabilities become matrix powers, long-run behaviour becomes a statement about eigenvalues, and the equilibrium of the chain becomes a left eigenvector.
The limit theorem that says the chain actually converges to that equilibrium is the subject of Convergence of Markov chains; the walk on Zd, the canonical example with an infinite state space, is treated in Random walks. Here we fix the objects and the bookkeeping.
Definition(The Markov property, transition matrices and n-step probabilities)AdaptedAdapted from Grimmett & Stirzaker (2020), §6.1, Definitions (6.1.1), (6.1.4) and (6.1.6), and Theorem (6.1.5), pp. 239-241
Let S be a countable set, the state space, whose elements we may and do take to be integers. A sequence X={Xn:n≥0} of random variables taking values in S, defined on a probability space (Ω,F,P), is a Markov chain if
for every n≥1 and all states s,x1,…,xn−1 for which the conditioning event has positive probability. Grimmett and Stirzaker record two equivalent forms: one may condition on the values at any finite increasing set of past epochs n1<n2<⋯<nk≤n and only the last one survives, and one may ask for Xm+n given X0,…,Xm and only Xm survives. The conditional probabilities here are those of Conditional expectation, taken with respect to the discrete conditioning events.
The chain is homogeneous if P(Xn+1=j∣Xn=i) does not depend on n. Its transition matrix is then P=(pij)i,j∈S with
pij=P(Xn+1=j∣Xn=i).
Every chain in this unit is homogeneous. The matrix P is stochastic: its entries are non-negative and each of its rows sums to 1, since from state i the chain must go somewhere. Conversely every stochastic matrix is the transition matrix of some chain.
The n-step transition probabilities are pij(n)=P(Xm+n=j∣Xm=i), which by homogeneity and the Chapman-Kolmogorov statement below do not depend on m; we write Pn for the matrix (pij(n)) and set pij(0)=δij. Finally, write μi(n)=P(Xn=i) and let μ(n) be the row vector of these numbers.
A warning on conventions: pij is the probability of a jump fromitoj, so that probability mass is propagated by multiplying a row vector on the right by P. Some older texts, and much of the physics literature on master equations, use the transposed convention in which the transition matrix acts on column vectors from the left. Grimmett and Stirzaker note that the transposed convention was formerly in use and warn the reader about it. We use the row convention throughout the handbook.
Definition(Chapman-Kolmogorov: n-step probabilities are matrix powers)AdaptedAdapted from Grimmett & Stirzaker (2020), §6.1, Theorem (6.1.7) and Lemma (6.1.8), pp. 241-242
Let X be a homogeneous Markov chain with transition matrix P. Then for all states i,j and all integers n,r≥0,
pij(n+r)=k∈S∑pik(n)pkj(r),
the Chapman-Kolmogorov equations. In matrix form Pn+r=PnPr, and therefore
Pn=Pn,n≥0.
Moreover the distribution of the chain evolves by μ(m+n)=μ(m)Pn, and in particular
μ(n)=μ(0)Pn.
Thus the law of the whole process is determined by the pair (μ(0),P).
Definition(Recurrence, mean recurrence time, period, irreducibility)AdaptedAdapted from Grimmett & Stirzaker (2020), §6.2, Definitions (6.2.1), (6.2.7), (6.2.8), (6.2.10), (6.2.11) and the first-passage probabilities (6.2.2), pp. 246-249; Grimmett & Stirzaker (2020), §6.3, Definition (6.3.1) and Theorem (6.3.2), pp. 250-251
Write Pi for P(⋅∣X0=i) and Ei for the corresponding expectation, and let
Tj=min{n≥1:Xn=j}
be the first passage time to j, with Tj=∞ if the chain never visits j after time 0. Put fij(n)=Pi(Tj=n) and fij=∑n≥1fij(n)=Pi(Tj<∞).
State i is recurrent (Grimmett and Stirzaker also write persistent) if fii=1, that is, if return to i from i is certain; otherwise fii<1 and i is transient.
The mean recurrence time of i is μi=Ei(Ti)=∑nnfii(n) when i is recurrent, and μi=∞ when i is transient. A recurrent state is positive (or non-null) if μi<∞ and null if μi=∞; a recurrent state can be null, so certainty of return does not imply a finite expected return time.
The period of i is d(i)=gcd{n≥1:pii(n)>0}. The state is aperiodic if d(i)=1 and periodic otherwise. A state that is recurrent, positive and aperiodic is called ergodic.
Say ileads toj, written i→j, if pij(m)>0 for some m≥0, and that i and jintercommunicate, i↔j, if each leads to the other. Intercommunication is an equivalence relation, and intercommunicating states share their period, their transience or recurrence, and their nullity or positivity. The chain is irreducible if S is a single equivalence class.
Grimmett and Stirzaker record the following for a finite state space S, with no irreducibility assumed: some state of the chain must be recurrent, and nullity is impossible, every recurrent state being positive. Combining this with the class properties, every state of a finite irreducible chain is positive recurrent.
Definition(When a state is recurrent, and the stationary distribution)AdaptedAdapted from Grimmett & Stirzaker (2020), §6.2, Theorem (6.2.3) and Corollary (6.2.4), pp. 247-248; Grimmett & Stirzaker (2020), §6.4, Definition (6.4.1) and Theorem (6.4.3) with equation (6.4.11), pp. 254-257; Grimmett & Stirzaker (2020), §6.6, Theorem (6.6.1), p. 269
(a) A summability criterion. Define the generating functions Pij(s)=∑n≥0snpij(n) and Fij(s)=∑n≥0snfij(n) for ∣s∣<1, with fij(0)=0. Then
Consequently state j is recurrent if and only if ∑npjj(n)=∞. If j is recurrent then ∑npij(n)=∞ for every i with fij>0; if j is transient then ∑npij(n)<∞, and hence pij(n)→0, for every i.
(b) Stationary distributions. A row vector π=(πj:j∈S) is a stationary distribution of the chain if πj≥0, ∑jπj=1, and
π=πP,that isπj=i∈S∑πipij for all j.
Iterating gives πPn=π for all n, so a chain started with X0∼π has Xn∼π for every n. Note that π is a left eigenvector of P with eigenvalue 1, not a right one.
(c) Existence and uniqueness. A stationary distribution exists for an irreducible chain exactly when every one of its states is positive recurrent. In that case π is unique and
πi=μi1for every i∈S,
where μi is the mean recurrence time. No assumption of aperiodicity is needed for this.
(d) The finite case. If S is finite with N=∣S∣ and the chain is irreducible with period d, the Perron-Frobenius theorem says that λ1=1 is an eigenvalue of P; that the d complex d-th roots of unity are eigenvalues; that the remaining eigenvalues λd+1,…,λN satisfy ∣λr∣<1; and that there is a unique π with πP=π, all of whose entries are strictly positive.
Idea of the proofOriginal
Two ideas carry the whole unit.
The first is that conditioning on where the chain is at an intermediate time decomposes a long journey into two independent-looking short ones. Summing over the intermediate state is exactly a matrix product, which is why Pn=Pn and why the long-run questions turn into linear algebra.
The second is that a renewal happens at every return to a fixed state k: by the Markov property the chain starts afresh there, so the excursions between consecutive visits to k are independent and identically distributed. The time between visits has mean μk, so in a long run of n steps the chain makes about n/μk visits to k, and the fraction of time spent at k is about 1/μk. That heuristic is the content of πk=1/μk, and it also explains the summability criterion: ∑npkk(n) is the expected number of visits to k, which is infinite exactly when return is certain and the visits never stop.
Proof(Chapman-Kolmogorov, the evolution of μ(n), and the summability criterion)AdaptedAdapted from Grimmett & Stirzaker (2020), §6.1, proofs of Theorem (6.1.7) and Lemma (6.1.8), pp. 241-242; Grimmett & Stirzaker (2020), §6.2, proofs of Theorem (6.2.3) and Corollary (6.2.4), pp. 247-248
Chapman-Kolmogorov. Fix i,j and m,n,r≥0. Partition the event {Xm+n+r=j} according to the state occupied at the intermediate epoch m+n:
which is legitimate because S is countable and the events are disjoint. Each summand factorises as
P(Xm+n+r=j∣Xm+n=k,Xm=i)P(Xm+n=k∣Xm=i),
by the elementary identity P(A∩B∣C)=P(A∣B∩C)P(B∣C). The Markov property, in the form that allows conditioning on any finite set of past epochs, lets us drop the condition Xm=i from the first factor; homogeneity then identifies the two factors as pkj(r) and pik(n). This is the displayed equation, and in matrix language Pn+r=PnPr. Taking r=1 and inducting on n gives Pn=Pn, and in particular the right-hand side does not depend on the starting epoch m.
which is the j-th entry of the row vector μ(m)Pn. Setting m=0 gives μ(n)=μ(0)Pn. □
The summability criterion. Fix i and j. Decompose a visit to j at time m≥1 according to the time r of the first visit to j. These m possibilities are disjoint and exhaust the event, so by the Markov property applied at the first-passage time,
pij(m)=r=1∑mfij(r)pjj(m−r),m≥1.
The right-hand side is a convolution. Multiplying by sm and summing over m≥1, and recalling pij(0)=δij, turns it into the product of generating functions (see Generating functions):
Pij(s)−δij=Fij(s)Pjj(s),
which is both displayed identities at once. Taking i=j and solving,
Pjj(s)=1−Fjj(s)1(∣s∣<1).
Now let s↑1. Both Pjj and Fjj have non-negative coefficients, so Abel's theorem lets us pass to the limit termwise: Fjj(s)→fjj and Pjj(s)→∑npjj(n), the two limits being finite or infinite together with the sums. Hence ∑npjj(n)=∞ precisely when fjj=1, that is, precisely when j is recurrent. The statements for i=j follow from Pij(s)=Fij(s)Pjj(s) by letting s↑1 again: the left side blows up iff Pjj does and fij=Fij(1)>0. □
Remark(Where each hypothesis is used)Original
Countability of S. It enters at the first line of the Chapman-Kolmogorov proof, where the event is split into a countable disjoint union and summed. For an uncountable state space the sum must be replaced by an integral against a transition kernel, and the result is the same equation in integral form.
Homogeneity. Used only to identify P(Xm+n+r=j∣Xm+n=k) with pkj(r). Without it the Chapman-Kolmogorov equations still hold in the form P(m,m+n+r)=P(m,m+n)P(m+n,m+n+r), but the matrices depend on the epoch and are no longer powers of a single P, so the eigenvalue analysis is lost.
The Markov property. Used exactly once in each proof: to delete the condition Xm=i in Chapman-Kolmogorov, and to restart the chain at the first passage time in the convolution identity. If it fails, the convolution identity fails and the generating-function criterion with it.
Irreducibility in part (c). Without it the conclusion is false in both directions. The two-state chain with p11=p22=1 is positive recurrent in every state yet has infinitely many stationary distributions, namely every π=(a,1−a); so uniqueness genuinely needs irreducibility. And a chain with an absorbing state 0 reachable from a transient part has a stationary distribution concentrated on 0 while other states are transient, so "has a stationary distribution" does not by itself make every state positive recurrent.
Positive recurrence in part (c). Null recurrence is not an exotic pathology: the symmetric simple random walk on Z is irreducible and recurrent, every μi is infinite, and it has no stationary distribution at all, only the stationary measure xi≡1, which cannot be normalised (see Random walks). This is where the finiteness of S does real work: by the lemma above, a finite irreducible chain cannot be null recurrent, so it always has exactly one π.
Aperiodicity. Deliberately not assumed anywhere in this unit. The chain on {1,2} with p12=p21=1 is irreducible and positive recurrent, has the unique stationary distribution π=(21,21), and yet p11(n) alternates between 1 and 0 and converges to nothing. Stationarity is an algebraic statement; convergence to π is a separate analytic statement that does need aperiodicity, and it is proved in Convergence of Markov chains.
Method(Analysing a chain given its transition matrix)Original
Write P down and check it is stochastic. Every row must sum to 1. An arithmetic slip here invalidates everything after it.
Draw the directed graph with an arrow i→j whenever pij>0, and find its communicating classes. The chain is irreducible exactly when the graph is strongly connected.
Classify. For a finite irreducible chain, stop: all states are positive recurrent (lemma above). Otherwise identify the closed recurrent classes and the transient states; from a transient state the chain eventually enters one of the closed classes and stays there, so analyse each closed class as a chain in its own right.
Find the period of one state in each class; all states in a class share it. In practice: if pii>0 for some i in the class, the class is aperiodic. Otherwise take the gcd of the lengths of the cycles through i in the graph.
Solve πP=π with ∑jπj=1. Because P is stochastic, the all-ones column vector satisfies P1=1, so P−I is singular and the N equations π(P−I)=0 are linearly dependent. Discard any one of them and append the normalisation ∑jπj=1 in its place.
Read off the mean recurrence times as μi=1/πi.
For n-step probabilities, diagonalise P if you can: with P=B−1ΛB one has Pn=B−1ΛnB, and the eigenvalues other than 1 control how fast pij(n) settles down. If P is not diagonalisable, use its Jordan form (Jordan normal form); the extra polynomial factors nk do not change which terms decay.
Example(A three-state chain: classification, stationary distribution, and a three-step probability)Original
A turbulence monitor samples a flow once per minute and records which of three regimes it is in: 1 = laminar, 2 = intermittent, 3 = fully turbulent. The record is modelled as a homogeneous Markov chain with
P=2131021314103143.
Step 1: P is stochastic. The row sums are 21+21=1, 31⋅3=1 and 41+43=1.
Step 2: irreducibility. The arrows are 1→2, 2→1, 2→3, 3→2, together with the loops at 1, 2 and 3. So 1↔2↔3 and the chain is irreducible. The state space is finite, so every state is positive recurrent.
Step 3: aperiodicity.p11=21>0, so d(1)=1; by the class property all three states are aperiodic, hence all three are ergodic.
Step 4: the stationary distribution. Solve πP=π, that is
The first equation gives 31π2=21π1, so π2=23π1. The third gives 31π2=41π3, so π3=34π2=2π1. The second equation is then automatically satisfied, as it must be. Normalising, π1(1+23+2)=1, so π1=92 and
π=(92,31,94).
A check: πP has first entry 92⋅21+31⋅31=91+91=92. Correct.
Step 5: mean recurrence times. By part (c) of the statement, μi=1/πi, so
μ1=29=4.5,μ2=3,μ3=49=2.25minutes.
The monitor returns to full turbulence on average every 2.25 minutes, and to laminar flow every 4.5 minutes.
Step 6: a three-step probability by Chapman-Kolmogorov. First square P. Its top row is
So a flow observed laminar now is fully turbulent three minutes later with probability 19/72. For comparison, the stationary value is π3=4/9≈0.4444: three steps is not yet long enough to forget the laminar start. The eigenvalues of P are 1 and 247±73, numerically 0.64767 and −0.06433; the largest of these in modulus is 0.64767, and since 0.64767n falls below 0.01 only at n=11, the chain needs of order ten minutes to settle. Direct computation confirms this: p13(8)=0.42387 and p13(20)=0.44433, against π3=0.44444.
Application(Kinetic Monte Carlo and Markov state models)Original
Discrete-time chains are the standard coarse-graining of a complicated continuous dynamics onto a manageable set of labels.
In Markov state models of molecular dynamics one partitions the configuration space of a protein into a few hundred metastable sets, runs many short molecular-dynamics trajectories, and estimates pij as the observed fraction of trajectories that are in set j a lag time τ after being in set i. The payoff is exactly the content of this unit: the stationary vector π gives the equilibrium populations of the conformations (and hence their free-energy differences, ΔFij=−kBTlog(πj/πi)), the mean recurrence times 1/πi give residence statistics, and the eigenvalues of P other than 1 give the relaxation timescales −τ/log∣λr∣ of the slow conformational transitions, all from short simulations that individually never see a full folding event.
In lattice models of surface growth and radiation damage, kinetic Monte Carlo evolves an adatom configuration by hopping moves; recording the configuration after each sweep gives a chain whose stationary distribution is the equilibrium ensemble. Whether such a chain actually relaxes to π, and how fast, is the subject of Convergence of Markov chains; the design principle that makes π come out right in the first place is detailed balance, treated there.
The same bookkeeping, with S the set of web pages and P the hyperlink matrix smoothed by a small uniform term, is the PageRank computation: the ranking vector is the left eigenvector π of P with eigenvalue 1, computed by the power iteration μ(n)=μ(0)Pn of the evolution lemma.
Proves by coupling that the n-step probabilities of an irreducible aperiodic chain converge to the reciprocal of the mean recurrence time, shows that a distribution in detailed balance with the transition matrix is automatically stationary and makes the chain reversible, and gives a geometric bound on the remaining total-variation distance in terms of the second-largest eigenvalue modulus.
Note(What is being proved, and what is not)Original
Discrete-time Markov chains produced a stationary distribution π as a solution of the algebraic equation π=πP. Nothing there said that a chain started anywhere else ever approaches π. This unit supplies the missing analytic statement, identifies the extra hypothesis it needs (aperiodicity), and then asks the practical question that follows: how long must one wait.
Three things are proved. First, the limit theorem: pij(n)→1/μj for an irreducible aperiodic chain, with no positivity assumption, so that the transient and null recurrent cases are included as the statement pij(n)→0. Second, a sufficient condition, detailed balance, which both identifies π cheaply and makes the chain reversible. Third, a geometric bound on the error, controlled by the second largest eigenvalue modulus of P.
Definition(Total variation distance, reversibility, detailed balance)AdaptedAdapted from Grimmett & Stirzaker (2020), §4.12, equations (4.12.7) and (4.12.8) with the footnote on normalisation, p. 145; Grimmett & Stirzaker (2020), §6.5, Definition (6.5.2) and the equations (6.5.3), pp. 265-266
Let X be an irreducible homogeneous Markov chain on a countable state space S with transition matrix P=(pij), and write pi⋅(n) for the mass function k↦pik(n).
The total variation distance between two mass functions f and g on S is
Convention. Grimmett and Stirzaker define dTV without a factor of one half, and note in a footnote that some authors halve it; much of the mixing-time literature uses the halved version, for which the supremum form has no factor 2 and the distance never exceeds 1. We follow Grimmett and Stirzaker throughout this unit, so our dTV lies in [0,2]. Halve every displayed bound to convert.
Suppose X has stationary distribution π and X0∼π, so that Xn∼π for every n. For a finite horizon N, the time-reversal of X is Yn=XN−n, 0≤n≤N. The chain is reversible if Y has the same transition matrix as X, which (by the time-reversal statement below) is the same as requiring the detailed balance equations
πipij=πjpjifor all i,j∈S.
More generally, a transition matrix P and a distribution λ are said to be in detailed balance if λipij=λjpji for all i,j; an irreducible chain whose P is in detailed balance with its stationary distribution is called reversible in equilibrium.
Theorem(The Markov chain limit theorem)AdaptedAdapted from Grimmett & Stirzaker (2020), §6.4(B), Example (6.4.19) and Theorem (6.4.20) with remarks (a)-(d), pp. 259-260; Grimmett & Stirzaker (2020), §6.4(A), Theorem (6.4.3), p. 254
Let X be a Markov chain on a countable state space S which is irreducible and aperiodic, with transition matrix P and mean recurrence times μj=Ej(Tj). Then
pij(n)⟶μj1as n→∞, for every i,j∈S,
with the convention 1/∞=0. Four consequences are worth separating.
When the chain is transient, or recurrent with every μj infinite, the limit is pij(n)→0 for all i,j: a limit exists, but it is not a probability distribution.
If the chain is positive recurrent then it has a unique stationary distribution π, with πj=1/μj, and pij(n)→πj.
The limit does not depend on the starting state i. Hence for any initial distribution μ(0),
P(Xn=j)=i∑μi(0)pij(n)⟶μj1,
so the chain forgets where it started.
Aperiodicity cannot be dropped, but it can be worked around. If X is irreducible with period d, then {Xnd:n≥0} is aperiodic and pjj(nd)→d/μj.
Theorem(Detailed balance, reversibility, and a cheap route to π)AdaptedAdapted from Grimmett & Stirzaker (2020), §6.5, Theorem (6.5.1), Definition (6.5.2) and Theorem (6.5.4), pp. 265-266
Let X be irreducible and positive recurrent with stationary distribution π, and let X0∼π.
(a) The reversed chain is a chain.Yn=XN−n is a Markov chain, with transition probabilities
P(Yn+1=j∣Yn=i)=πiπjpji.
In particular X is reversible exactly when πipij=πjpji for all i,j.
(b) Detailed balance implies stationarity. Suppose an irreducible chain with transition matrix P admits some distribution π for which the balance relation πipij=πjpji holds at every pair of states i,j∈S. That π is then stationary, and the chain is reversible in equilibrium.
This is the practical content: detailed balance is a system of two-term equations, one for each edge of the transition graph, whereas π=πP couples every state to every other. When a guess satisfies detailed balance one has found π with no linear algebra at all. The converse fails: a stationary distribution need not be in detailed balance with P, and a chain that is not reversible still has its π.
Theorem(Geometric rate of convergence for a finite reversible chain)AdaptedAdapted from Grimmett & Stirzaker (2020), §6.6, Theorem (6.6.1) and the bound following it, pp. 269-270; Grimmett & Stirzaker (2020), §6.14, Theorem (6.14.9) and the expansion (6.14.11), pp. 333-334
Let S be finite with ∣S∣=N and let X be irreducible and aperiodic with transition matrix P and stationary distribution π. By the Perron-Frobenius theorem the eigenvalues of P may be listed as λ1=1,λ2,…,λN with ∣λr∣<1 for r≥2; order them so that λ2 has the largest modulus among these. The quantity ∣λ2∣ is called the second largest eigenvalue modulus, or SLEM. If the eigenvalues are distinct, so that P=B−1ΛB is diagonalisable, then Pn=B−1ΛnB, every row of Pn tends to π, and there is a constant c<∞ with
∣pij(n)−πj∣≤c∣λ2∣nfor all i,j and all n.
Grimmett and Stirzaker state the bound in this diagonalisable setting; more generally, if λ2 has multiplicity m one must allow a polynomial factor, and the entries of Pn approach their limits with an error of order nm−1∣λ2∣n.
If in addition X is reversible in equilibrium, the constant can be made explicit. Equip RN with the inner product ⟨x,y⟩=∑i∈Sπixiyi; reversibility makes P self-adjoint for this inner product, so its eigenvalues are real and its right eigenvectors v1,…,vN may be taken to form an orthonormal basis, with v1≡1. Writing vr(i) for the i-th entry of vr,
The left-hand side is the total variation distance between pi⋅(n) and π. So convergence is geometric at rate ∣λ2∣, and the number of steps needed to bring the distance below a tolerance ε grows like log(1/ε)/log(1/∣λ2∣). A chain with ∣λ2∣ close to 1 is said to be slowly mixing.
Idea of the proofOriginal
The limit theorem is proved by coupling, which replaces an analytic estimate by a construction. Run two copies of the chain side by side and independently, one started at i and one started at j. If the two copies ever occupy the same state at the same time, then from that moment on they have the same conditional law, because a chain's future depends on nothing but its present position and the matrix P. So the difference pik(n)−pjk(n) can only come from the event that they have not yet met, and it is enough to show that they meet almost surely. The pair is itself a Markov chain on S×S, and it is irreducible exactly because the original chain is aperiodic: if X had period 2, the two copies started at even and odd distance could never be at the same place at the same time. That one sentence is where aperiodicity is used, and it is the only place.
Once the limit is known to be independent of the start, averaging over the start against π identifies it as π itself, because averaging a constant sequence changes nothing.
Detailed balance is a different kind of statement and needs no limit at all. It says that in equilibrium the probability flux from i to j matches the flux back, edge by edge. Stationarity only demands that the total flux into each node matches the total flux out. Local balance obviously implies global balance, and that one-line implication is the whole proof.
Proof(The limit theorem, by coupling)AdaptedAdapted from Grimmett & Stirzaker (2020), §6.4(B), proof of Theorem (6.4.20) including the coupling inequality (6.4.21) and the bounded convergence step (6.4.22), pp. 260-262
Assume first that X is positive recurrent; the transient and null recurrent cases are dealt with at the end.
Step 1: the coupled chain. Let X={Xn} and Y={Yn} be independent chains, each with state space S and transition matrix P, and set Zn=(Xn,Yn). Then Z is a Markov chain on S×S with transition probabilities
P(Zn+1=(k,l)∣Zn=(i,j))=pikpjl,
the product arising from the independence of the two copies.
Step 2: Z is irreducible — this is where aperiodicity enters. For an irreducible aperiodic chain and any pair of states a,b there is an N(a,b) with pab(n)>0 for all n≥N(a,b): the set {n:paa(n)>0} is closed under addition and has greatest common divisor 1, so it contains every sufficiently large integer, and composing with one path from a to b shifts this to b. Hence for any (i,j) and (k,l) and all n beyond the maximum of the two thresholds, pik(n)pjl(n)>0, so Z can get from (i,j) to (k,l). Without aperiodicity this fails: for the chain on {1,2} with p12=p21=1, the pair started at (1,2) visits only (1,2) and (2,1).
Step 3: Z is positive recurrent.X has a unique stationary distribution π. The product νij=πiπj is a distribution on S×S, and it is stationary for Z because
Since Z is irreducible and possesses a stationary distribution, it is positive recurrent.
Step 4: the two copies meet. Fix any s∈S, start the pair at Z0=(i,j), and let
T=min{n≥1:Zn=(s,s)}.
By recurrence of the irreducible chain Z, the state (s,s) is reached from (i,j) with probability one, so P(T<∞)=1.
Step 5: the coupling inequality. Suppose T≤n. At time T the two copies are both at s; from then on, their conditional distributions coincide, since each depends only on P and on the common current position. Therefore, conditionally on {T≤n}, Xn and Yn have the same law, and
The summands are bounded by 1 and tend to 0. To pass to the limit on an infinite S, split the sum at a finite set F: the part over F is a finite sum of terms tending to zero, and the part outside is at most 2∑i∈/Fπi, which enlarging F forces below any prescribed tolerance, because π is a probability distribution. Hence pjk(n)→πk=1/μk.
Step 7: the remaining cases. If X is transient then ∑npij(n)<∞ for all i,j, so pij(n)→0=1/μj. If X is null recurrent the coupled chain Z cannot be positive recurrent: starting from (i,i), the return time of Z to (i,i) is at least the return time of X to i, whose mean is infinite. If Z is transient then pij(n)2=P(Zn=(j,j)∣Z0=(i,i))→0 and we are done; if Z is null recurrent, Step 5 still applies, and Grimmett and Stirzaker complete the argument by a diagonal-selection device which shows that a non-zero subsequential limit would produce a stationary distribution and so contradict null recurrence. □
Proof(Reversal and detailed balance)AdaptedAdapted from Grimmett & Stirzaker (2020), §6.5, proofs of Theorem (6.5.1) and Theorem (6.5.4), pp. 265-266
(a) The reversed chain. Write a trajectory of Y as a trajectory of X read backwards. For states i0,…,in+1,
Because X0∼π and π is stationary, every Xm has distribution π, and the numerator is the product πin+1pin+1inpinin−1⋯pi1i0 while the denominator is πinpinin−1⋯pi1i0. All the factors pimim−1 with m≤n cancel, leaving πin+1pin+1in/πin. This depends on in and in+1 only, which proves both that Y is Markov and that its transition probabilities are (πj/πi)pji. Equality of these with pij is exactly πipij=πjpji.
(b) Detailed balance implies stationarity. Suppose π is a distribution with πipij=πjpji for all i,j. Sum this over i, holding j fixed:
i∈S∑πipij=i∈S∑πjpji=πji∈S∑pji=πj,
the last step because P is stochastic, so row j sums to one. The left side is the j-th entry of πP; hence πP=π and π is stationary. Reversibility then follows from part (a). □
Remark(Where the hypotheses are used, and what fails without them)Original
Aperiodicity is used once, in Step 2, to make the coupled chain irreducible. Its necessity is visible in the simplest example: on S={1,2} with p12=p21=1, the chain is irreducible and positive recurrent, π=(21,21) is the unique stationary distribution, and yet p11(n) equals 1 for even n and 0 for odd n. The Cesaro averages still converge to π, and part 4 of the limit theorem recovers the full statement along multiples of the period. A standard repair is to make the chain lazy: replace P by L=aP+(1−a)I with 0<a<1. Then L is irreducible and aperiodic, has the same stationary distribution, and is reversible whenever P is.
Irreducibility is used in Step 4 (to reach (s,s)) and in Step 3 (uniqueness of π). Without it, different starting states can have genuinely different limits: a chain with two absorbing states converges to whichever one it is started in.
Positivity of π in Step 6 is where the chain's being a probability distribution matters. For a null recurrent chain there is still a stationary measure, but it has infinite total mass and cannot be used to average; that is why Step 7 needs a separate argument.
Independence of the two copies in Step 1 is a choice, not a necessity. The coupling only has to be a joint construction of two processes with the right marginals which meet quickly; cleverer, dependent couplings are how sharp mixing-time bounds are actually obtained.
Reversibility is not needed for the limit theorem at all; it is needed only for the explicit constant in the spectral bound, because it is what makes P self-adjoint for the inner product ⟨x,y⟩=∑iπixiyi and hence gives real eigenvalues and an orthonormal eigenbasis (see The spectral theorem for Hermitian operators). An irreversible P may have complex eigenvalues and may fail to be diagonalisable; the bound c∣λ2∣n survives, with a polynomial correction nm−1 when λ2 has multiplicity m in a Jordan block, but the constant c is no longer explicit.
Finiteness of S is needed for the spectral statements only. The limit theorem itself holds on any countable S, and there the rate can be anything at all: an infinite chain may converge sub-geometrically.
Method(Getting a chain to converge to a prescribed π, and estimating how long it takes)Original
Given a target distribution π on a finite set S:
Choose a proposal matrixH=(hij) that is irreducible and easy to sample from; symmetry (hij=hji) is convenient but not required.
Accept with the Hastings probability. From state i, draw a proposal j with probability hij and accept it with probability
aij=min{1,πihijπjhji};
on rejection stay at i. The resulting matrix is pij=hijaij for i=j, with pii whatever is left over. Only ratios πj/πi are needed, so an unknown normalising constant in π is irrelevant. When H is symmetric this reduces to the Metropolis rule aij=min{1,πj/πi}.
Verify detailed balance, which is immediate from the form of aij: for i=j, πihijaij=min{πihij,πjhji}, an expression symmetric in i and j. By the statement above, π is then stationary.
Check aperiodicity. If any proposal is ever rejected with positive probability then some pii>0 and the chain is aperiodic. If not, make the chain lazy.
Estimate the rate. Compute, or bound, the SLEM ∣λ2∣. The number of steps to reach total variation distance ε is of order log(1/ε)/log(1/∣λ2∣), and the reversible bound above turns this into an inequality with explicit constants.
Discard the burn-in. Estimate ∑θg(θ)πθ by n−1∑r=MM+n−1g(Xr), with M chosen large enough that the bound in step 5 is below the accuracy wanted.
Example(A lazy Ehrenfest chain: exact convergence against the spectral bound)Original
A partition with a small hole separates two chambers, A and B, which hold m=2 gas molecules between them. A fair coin is tossed at each tick of a clock; heads leaves the configuration alone, while tails sends one molecule, chosen uniformly from the two, through the hole. Take Xn to be chamber A's occupancy at time n, so that S={0,1,2}. From state i the chosen molecule is in A with probability i/2, so
pi,i−1=21⋅2i,pi,i+1=21(1−2i),pii=21,
that is
P=2141021212104121.
Stationary distribution by detailed balance. Rather than solve π=πP, look for π with πipi,i+1=πi+1pi+1,i; for a chain that only moves to neighbours these are the only non-trivial detailed balance equations. They read π0⋅21=π1⋅41 and π1⋅41=π2⋅21, giving π1=2π0 and π2=π0. Normalising,
π=(41,21,41),
the binomial distribution (i2)2−2, as physical intuition demands. By the detailed balance theorem π is stationary and the chain is reversible in equilibrium. The chain is irreducible and p00=21>0, so it is aperiodic, and the limit theorem applies: the mean recurrence times are μ0=μ2=4 and μ1=2 ticks.
Eigenvalues and eigenvectors. One checks directly that
P111=1⋅111,P10−1=2110−1,P1−11=0.
For instance the middle one: the first row gives 21⋅1+21⋅0=21, the second 41⋅1+41⋅(−1)=0, the third 21⋅0+21⋅(−1)=−21. So the spectrum is {1,21,0} and the SLEM is ∣λ2∣=21. Normalising in the inner product ⟨x,y⟩=∑iπixiyi: the first vector has ⟨v,v⟩=41+21+41=1; the second has 41⋅1+21⋅0+41⋅1=21, so v2=2(1,0,−1); the third has 41+21+41=1 already. They are mutually orthogonal, as reversibility guarantees.
Exact n-step probabilities. Expanding pik(n)=∑rvr(k)πkλrnvr(i) from state i=0:
valid for n≥1. (Check at n=1: 41+41=21, 21, 0 — the first row of P.) Hence the exact total variation distance from the stationary distribution, started from the empty chamber, is
dTV(p0⋅(n),π)=2−(n+1)+0+2−(n+1)=2−n.
At n=5 this equals 321=0.03125.
Comparison with the bound. The reversible bound gives, with N=3, ∣λ2∣=21 and supr∣vr(0)∣=max{1,2,1}=2,
dTV(p0⋅(n),π)≤32⋅2−n,
which at n=5 is 32/32=0.1326. The bound holds, and overestimates the true distance by the factor 32≈4.24 — it gets the geometric rate exactly right and only the prefactor is loose. To force dTV≤10−3 the bound demands n≥log2(32×103)=12.05, so n=13; the truth is n=10.
Application(Markov chain Monte Carlo and detailed balance in statistical physics)Original
Every Monte Carlo simulation of a thermal system is an instance of this unit. One wants averages ⟨g⟩=∑θg(θ)π(θ) over the Boltzmann distribution π(θ)=Z−1exp(−βE(θ)) on a configuration space far too large to enumerate, and the partition function Z is unknown. The Metropolis recipe in the method above builds a chain whose stationary distribution is exactly π, using only the energy differences E(θ′)−E(θ), because the acceptance probability min{1,e−βΔE} depends on π only through ratios and Z cancels. Detailed balance is the condition physicists impose on the update rule precisely because it guarantees, by the theorem above, that the Boltzmann weight is left invariant; it is a sufficient condition, imposed for convenience, not a law of nature, and schemes that are stationary without being reversible are also in use.
The limit theorem is then the statement that the simulation gives the right answer at all, and aperiodicity is the fine print: a deterministic sweep through the lattice in a model with no self-transitions can be periodic, which is one reason random site selection or a laziness parameter is standard.
The rate statement is where the physics gets hard. Near a continuous phase transition the SLEM of a local-update chain for the Ising model approaches 1 as the system size grows, the relaxation time diverges as Lz with a dynamical exponent z≈2 for Metropolis dynamics, and the simulation suffers critical slowing down. Cluster algorithms are designed to keep detailed balance while drastically reducing ∣λ2∣. The eigenvalue λ2 and its eigenvector also carry physical meaning away from criticality: they identify the slowest collective mode of the dynamics, which in a chain built from molecular dynamics (see Discrete-time Markov chains) is the rate-limiting conformational change.
Finally, the same bound tells an experimentalist when an equilibration protocol is long enough. For the Ehrenfest model of the example, generalised to m molecules, the SLEM of the lazy chain is 1−1/m, so equilibration takes of order mlogm ticks: the characteristic time grows only a little faster than the number of molecules, which is why gases equilibrate fast.
Solves the gambler's ruin problem for a simple random walk between two absorbing barriers, and proves Polya's trichotomy: the symmetric nearest-neighbour walk returns to its starting point with probability one on the line and on the plane, but not in space, where the return probability is about 0.34.
Definition(Simple random walk, on the line and on the lattice)AdaptedAdapted from Grimmett & Stirzaker (2020), §3.9, equations (3.9.1)-(3.9.2) and Lemmas (3.9.3)-(3.9.5), pp. 77-79; Grimmett & Stirzaker (2020), Problem (6.15.9), p. 340
Let X1,X2,… be independent, all carrying the same two-point law, which puts mass p on the value +1 and mass q=1−p on the value −1, with 0<p<1. The simple random walk started at a∈Z is the sequence
S0=a,Sn=a+i=1∑nXi(n≥1).
The walk is symmetric if p=q=21. It is spatially homogeneous (P(Sn=j∣S0=a) depends only on j−a), temporally homogeneous, and has the Markov property, so it is a Markov chain on Z in the sense of Discrete-time Markov chains, with transition probabilities pi,i+1=p, pi,i−1=q and all others zero.
Two variants are needed below.
Absorbing barriers. Fix an integer N≥1 and restrict to the state space {0,1,…,N}, declaring 0 and N to be absorbing: p00=pNN=1. In gambling language, a player holding k units stakes one unit on each of a sequence of independent bets won with probability p, and stops on reaching 0 (ruin) or N (the target).
The lattice walk in d dimensions. Let S0=0 and Sn=∑i≤nϵi, where the steps ϵi are independent and uniform on the 2d unit vectors ±e1,…,±ed of Zd, each with probability 1/(2d). For d=1 this is the symmetric walk above.
Occupation and first-return probabilities. Whichever of these unrestricted walks is under discussion — the walk on Z with parameter p, started at a=0, or the symmetric lattice walk in d dimensions — we write
p0(n)=P(walk is at the origin at time n),f0(n)=P(the first return to the origin occurs at time n),
with p0(0)=1 and f0(0)=0. In every case a single step changes the parity of the sum of the coordinates, so p0(n)=0 for odd n and every state has period 2.
The walk is called recurrent if it returns to its starting point with probability one and transient otherwise; by spatial homogeneity this does not depend on the starting point, and it is the recurrence of the origin as a state of the chain.
Theorem(Gambler's ruin)AdaptedAdapted from Grimmett & Stirzaker (2020), §3.9, Example (3.9.6), the difference equations (3.9.7) and (3.9.8) and their solutions, including (3.9.9), pp. 79-80; Grimmett & Stirzaker (2020), Exercise (3.9.1), p. 80
Take the walk of the definition, confined as described to the states 0,1,…,N and absorbed at whichever end it reaches first, and start it from k, where 0≤k≤N. Write Pk and Ek for probabilities and expectations of the walk so started, write pk for the chance that it finishes at 0 rather than at N (the ruin probability, indexed by the starting state; it has nothing to do with the occupation probabilities p0(n) of the definition, which are indexed by time, and the clash of the two at p0=p0(0)=1 is a coincidence of values), write Dk for the mean number of steps it makes before stopping, and put ρ=q/p. Then
Absorption is certain, and the time T it takes has a geometric tail: writing δ=min(p,q)N>0,
Pk(T>mN)≤(1−δ)m(m=0,1,2,…).
In particular Pk(T<∞)=1 and every moment of T is finite, so Dk=Ek(T) is the finite number displayed above, and 1−pk is the probability of reaching the upper barrier N.
Theorem(Polya's trichotomy: recurrence in one and two dimensions, transience in three)AdaptedAdapted from Grimmett & Stirzaker (2020), §5.3, Theorem (5.3.1) and Corollary (5.3.4), pp. 183-185; Grimmett & Stirzaker (2020), §6.2, Corollary (6.2.4) and Example (6.2.12), pp. 248-249; Grimmett & Stirzaker (2020), Problems (5.12.5), (5.12.6), p. 231, and Problem (6.15.9), p. 340
(a) The criterion, specialised to the walk. By the summability criterion of Discrete-time Markov chains, a state j is recurrent exactly when ∑npjj(n)=∞. Here p0(n)=0 at odd n, so the lattice walk is recurrent exactly when
n≥0∑p0(2n)=∞,
and this sum is the expected number of visits the walk makes to the origin.
(b) One dimension, general p. For the simple random walk on Z, the generating functions P0(s)=∑nsnp0(n) and F0(s)=∑nsnf0(n) of the occupation and first-return probabilities satisfy
so that the chance of a return to the origin ever occurring is
F0(1)=1−∣p−q∣.
Recurrence therefore holds precisely at p=21, and there the mean time until the walk first comes back diverges: the symmetric walk on Z is null recurrent, and every state has period 2. For p=21 every state is transient. Equivalently, p0(2n)∼(4pq)n/πn, which is summable precisely when 4pq<1, that is when p=21.
(c) Two dimensions. For the symmetric walk on Z2,
p0(2n)=[(n2n)4−n]2∼πn1,
so ∑np0(2n)=∞ and the walk is recurrent: it returns to the origin with probability one.
(d) Three dimensions. For the symmetric walk on Z3,
p0(2n)=(n2n)2−2ni+j+k=n∑(3ni!j!k!n!)2≤n3/2C
for a finite constant C, so ∑np0(2n)<∞ and the walk is transient.
Idea of the proofOriginal
Both halves of the unit come from the Markov property, used in two different ways.
For gambler's ruin, condition on the first bet. The walk started at k is, after one step, a walk started at k+1 or at k−1, with no memory of how it got there. So the unknown pk satisfies a linear second-order difference equation in k — the discrete analogue of Laplace's equation — with boundary values 1 at 0 and 0 at N. Its characteristic equation pθ2−θ+q=0 has the roots 1 and q/p. When p=q the two roots are distinct and the powers 1k and (q/p)k span the solutions; in the symmetric case the root 1 is repeated, the second solution is k rather than a second power, and the general solution reads A1+A2k. The expected duration satisfies the same equation with an inhomogeneous term 1, because one step has been used up.
For recurrence, condition instead on the time of first return. That turns the sequence p0(n) into the convolution of itself with the first-return sequence, and convolution is multiplication of generating functions, which inverts to P0=1/(1−F0). Letting s↑1 converts the question "is return certain?" into the question "is ∑np0(n) infinite?" — a question about the size of the walk, not its combinatorics. And the size is set by the central limit scaling: after n steps the walk is spread over a region of diameter n in each direction, so the chance of being at any one named site, the origin included, is of order n−d/2. The series ∑n−d/2 diverges for d=1,2 and converges for d≥3. Everything below is this one-line dimension count made exact.
Proof(Gambler's ruin)AdaptedAdapted from Grimmett & Stirzaker (2020), §3.9, Example (3.9.6) and the solution of (3.9.7) and (3.9.8), pp. 79-80
Absorption is certain and Dk is finite. This has to come first, because the difference equations below are solved in the real numbers, which presupposes that their solutions are finite; the recursions themselves are also satisfied by D≡∞. Let T be the absorption time and δ=min(p,q)N. Whatever the walk has done so far, if it is not yet absorbed it sits at some i with 0<i<N, and if its next N steps are all +1 it reaches N after N−i<N of them and is absorbed. That block of steps has probability pN≥δ, and the steps are independent of the past, so from any position not yet absorbed the chance of absorption within the next N steps is at least δ. Iterating over m disjoint blocks of N steps gives Pk(T>mN)≤(1−δ)m. Hence Pk(T=∞)=0, and summing the tail,
Dk=Ek(T)=n≥0∑Pk(T>n)≤Nm≥0∑(1−δ)m=δN<∞.
(The same bound makes every moment of T finite, which is the content of Exercise (3.9.1) of the source.) Since absorption is certain, pk and 1−pk are the probabilities of the two alternatives, finishing at 0 and finishing at N.
The difference equation. Let 1≤k≤N−1 and condition on the first step. With probability p the walk moves to k+1, and by the Markov property and temporal homogeneity it is then absorbed at 0 with probability pk+1; with probability q it moves to k−1 and is absorbed at 0 with probability pk−1. Hence
pk=ppk+1+qpk−1,1≤k≤N−1,
with p0=1 and pN=0 because those states are absorbing.
Solving it. Substituting the trial solution pk=θk and dividing by θk−1 gives the characteristic equation pθ2−θ+q=0, whose roots are θ=1 and θ=q/p=ρ, since p+q=1.
Case p=21. The roots are distinct, so the general solution is pk=A1+A2ρk. The boundary conditions give A1+A2=1 and A1+A2ρN=0, so A2=1/(1−ρN) and A1=−ρN/(1−ρN), which is the stated formula.
Case p=21. The root θ=1 is repeated, so the general solution is pk=A1+A2k, and the boundary conditions give A1=1, A2=−1/N.
The duration. Condition on the first step again, this time keeping track of the step taken. From k, one step is spent and the walk then needs Dk±1 further steps on average, so
with D0=DN=0. This is the same homogeneous equation with the constant forcing term 1, so the general solution is a particular solution plus A1+A2ρk (or A1+A2k when p=21).
Case p=21. Try the particular solution Dk=ck. Then ck=1+cp(k+1)+cq(k−1)=1+ck+c(p−q), so c=1/(q−p). Imposing D0=DN=0 on Dk=k/(q−p)+A1+A2ρk gives A1+A2=0 and N/(q−p)+A1+A2ρN=0, so A2=−A1 and A1(1−ρN)=−N/(q−p). Substituting back,
Dk=q−p1(k−N1−ρN1−ρk).
Case p=21. Now p−q=0, so a linear particular solution cannot work; try Dk=ck2, giving ck2=1+2c[(k+1)2+(k−1)2]=1+ck2+c, so c=−1. With Dk=−k2+A1+A2k and D0=DN=0 we get A1=0 and A2=N, that is Dk=k(N−k).
Both solutions are legitimate because the tail bound at the head of the proof has already made pk and Dk finite. □
Proof(Recurrence in one and two dimensions, transience in three)AdaptedAdapted from Grimmett & Stirzaker (2020), §5.3, Theorem (5.3.1) and Corollary (5.3.4), pp. 183-185; Grimmett & Stirzaker (2020), §6.2, Corollary (6.2.4), p. 248; Grimmett & Stirzaker (2020), Problems (5.12.5), (5.12.6), p. 231, and Problem (6.15.9), p. 340
The criterion. This part is taken over from the prerequisite rather than reproved: the renewal decomposition and its generating-function consequence are established for a general chain in Discrete-time Markov chains, and we record only what they say for the walk. Splitting a visit to the origin at time n≥1 according to the time r of the first return gives n disjoint possibilities that exhaust the event. On each of them the walk starts afresh at the origin at time r — note that {first return at r} is not an event about the state at a single time but a condition on the whole path up to r, so the Markov property is applied at the fixed time r to each such path separately and the contributions are then summed — and therefore
p0(n)=r=1∑nf0(r)p0(n−r).
Multiplying by sn, summing over n≥1 and using p0(0)=1 turns this convolution into P0(s)=1+F0(s)P0(s), that is P0(s)=1/(1−F0(s)) for ∣s∣<1; letting s↑1 termwise, which Abel's theorem permits because all coefficients are non-negative,
n≥0∑p0(n)=1−f1,f:=n≥1∑f0(n)=P(ever return),
read with the convention 1/0=∞: if f=1 then F0(s)↑1 and P0(s)=1/(1−F0(s))↑∞, so the left-hand side, a series of non-negative terms, diverges. So the sum is infinite exactly when f=1, which is part (a); the left-hand side is also the expected number of visits to the origin. What is left to do, and what the rest of the proof does, is to compute p0(n) in each dimension.
(b) One dimension. The walk is at the origin at time n only if n=2m is even and exactly m of the steps were +1, so
p0(2m)=(m2m)(pq)m.
The generating function of (m2m)xm is (1−4x)−1/2, so P0(s)=(1−4pqs2)−1/2 and therefore F0(s)=1−1/P0(s)=1−(1−4pqs2)1/2. At s=1,
f=F0(1)=1−1−4pq=1−(p+q)2−4pq=1−∣p−q∣,
using p+q=1. Hence f=1 if and only if p=q=21. When p=21, F0(s)=1−(1−s2)1/2 and F0′(s)=s(1−s2)−1/2→∞ as s↑1, so the mean return time is infinite and the walk is null recurrent. The same conclusion follows from the criterion: Stirling's formula m!∼mm+1/2e−m2π gives (m2m)∼4m/πm and hence
p0(2m)∼πm(4pq)m,
which sums to ∞ iff 4pq=1, that is iff p=21.
(c) Two dimensions. The walk is at the origin at time 2n exactly when, for some k, it made k steps east, k west, n−k north and n−k south. The number of orderings is the multinomial coefficient, and each of the 42n step sequences is equally likely, so
p0(2n)=4−2nk=0∑nk!k!(n−k)!(n−k)!(2n)!.
Factor the summand as (n2n)(kn)2, which is an identity between factorials:
By Vandermonde's identity ∑k=0n(kn)2=(n2n), so the sum collapses and
p0(2n)=4−2n(n2n)2=[(n2n)4−n]2∼(πn1)2=πn1,
again by Stirling. Since ∑n1/(πn) diverges, the criterion gives recurrence, and the return probability is 1.
(d) Three dimensions. Now the walk is at the origin at time 2n exactly when, for some i,j,k≥0 with i+j+k=n, it made i steps along each of ±e1, j along each of ±e2 and k along each of ±e3. There are 62n equally likely step sequences, so
where the second form follows by writing 6−2n=4−n9−n and (2n)!/(i!j!k!)2=(n2n)(n!/(i!j!k!))2. The rewriting is the point: the numbers wijk=n!/(3ni!j!k!) are non-negative and sum to 1, because ∑i+j+k=nn!/(i!j!k!)=3n by the multinomial theorem. Therefore
The multinomial coefficient n!/(i!j!k!) is largest when the three parts are as equal as possible, since replacing a pair (i,j) with i≥j+2 by (i−1,j+1) multiplies it by i/(j+1)>1. So the maximum is attained at the nearest-to-equal split, whose parts differ from n/3 by less than 1. For n divisible by 3 that split is (n/3,n/3,n/3), and Stirling applied to n!/((n/3)!)3 gives
For the other two residues each part is n/3+O(1), which changes the Stirling estimate by a factor 1+O(1/n) only, so the same asymptotic holds for every n; numerically the maximum is in fact slightly smaller off the divisible case, nmaxi+j+k=nwijk being 0.8158 at n=121 and 0.8224 at n=120, both below the limiting value 33/(2π)=0.8270.
Combining this with (n2n)2−2n∼1/πn gives
p0(2n)≤n3/2C(1+o(1)),C=2ππ33=0.4666.
The series ∑nn−3/2 converges, so ∑np0(2n)<∞ and the origin is transient. By spatial homogeneity every state is transient. □
Remark(Where the hypotheses are used, and how sharp the conclusions are)Original
Symmetry in parts (c) and (d). Only the symmetric walk is treated there, and for good reason: a lattice walk with any non-zero mean drift is transient in every dimension, because the strong law of large numbers of Laws of large numbers already pushes Sn/n away from the origin. Recurrence in d=1,2 is a statement about the borderline case of zero drift.
Independence of the steps. Used in the counting: every sequence of 2n steps has the same probability. For a self-avoiding walk, or any walk with memory, none of the formulas survive, and even the two-dimensional case becomes a hard open problem.
Finiteness of N in gambler's ruin. Letting N→∞ at fixed k in the ruin formula recovers the one-barrier answer: for p>21 we have ρ<1, so pk→ρk<1 and the player survives forever with positive probability, while for p≤21 we get pk→1 and ruin is certain. The mean duration does not follow that pattern, and the three cases must be kept apart. For p<21 we have ρ>1, so in Dk the ratio N(1−ρk)/(1−ρN) has a numerator linear in N and a denominator geometric in N and therefore tends to 0, leaving
Dk⟶q−pk(p<21):
ruin is certain and quick (at p=52, k=5 the limit is 25 steps, and the finite-N values 19.18, 24.80, 25.00 at N=10,20,50 climb to it). At p=21, by contrast, Dk=k(N−k)→∞: ruin is certain while the expected time to it is infinite, so there "ruin is certain" and "ruin happens quickly" come apart. That single case is the null recurrence of part (b) seen from another side; at p<21 the walk is transient, not null recurrent, and the analogy does not apply to it. (For p>21 also Dk→∞, but for the opposite reason: in the limit absorption at 0 is not certain at all.)
Higher dimensions. The cited treatment stops at d=3, and what follows is our own extension. The central limit scaling of the idea block can be made local: after 2n steps of the symmetric walk on Zd each coordinate has variance 2n/d, the coordinates are uncorrelated, and the Gaussian density at the origin is (4πn/d)−d/2, which has to be shared between the two parity classes of sites, giving
p0(2n)∼2(4πnd)d/2.
For d=1,2,3 this returns 1/πn, 1/(πn) and 0.2333n−3/2, the values computed above, and for every d≥3 it makes ∑np0(2n) converge: all the higher-dimensional walks are transient too, with the exponent d/2 in place of 3/2.
The exponent versus the constant. The estimate in (d) is an upper bound, not an asymptotic. Its constant C=0.4666 is exactly twice the true asymptotic constant 33/(4ππ)=0.2333, because bounding one factor wijk by its maximum throws away a factor 2 in the limit: a direct computation of the exact multinomial sum gives p0(80)=9.136×10−4 against the bound C⋅40−3/2=1.844×10−3. Transience needs only the exponent 3/2, so the slack costs nothing.
How transient is three-dimensional space? The criterion gives the return probability exactly: f=1−1/G with G=∑n≥0p0(n) the expected number of visits. By the Fourier inversion formula for a lattice-valued sum (see Characteristic functions, and Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §30.6, eq. (30.64), p. 1155) for the same representation of the displacement of a general multi-step walk), p0(n)=(2π)−3∫[−π,π]3φ(k)nd3k with single-step characteristic function φ(k)=31(cosk1+cosk2+cosk3), so summing the geometric series gives
Numerical integration gives G=1.51639 and hence f=0.34054: a walker on the cubic lattice comes back to where it started about one time in three, and otherwise never. The partial sums of the exact series agree: ∑n=140p0(2n)=0.44330, and estimating the remainder from the true asymptotic constant C′=33/(4ππ)=0.2333 of the previous paragraph,
n>40∑p0(2n)≈∫40∞n3/2C′dn=402C′=0.0738,
which brings the total to 0.5171 against the exact G−1=0.51639.
Periodicity. Every state of the lattice walk has period 2, so the limit theorem of Convergence of Markov chains applies only along even times — and in any case its conclusion here is the uninformative p0(2n)→0, since the walk is null recurrent or transient and has no stationary distribution.
Method(First-passage and absorption problems by conditioning on the first step)Original
The gambler's ruin calculation is a template. For a Markov chain on a finite or semi-infinite set of integers with nearest-neighbour steps:
Name the unknown as a function of the starting state, hk, and make sure it is a quantity to which the Markov property applies: an absorption probability, an expected hitting time, an expected reward collected before absorption.
Condition on the first step. Each possible first step leads to a state from which the chain starts afresh, so hk=∑jpkjhj for an absorption probability, and hk=1+∑jpkjhj for an expected time (the 1 being the step just taken).
Write the boundary conditions at the absorbing states: h=1 or 0 for a probability, h=0 for a hitting time.
Solve the difference equation. Substitute hk=θk into the homogeneous part to get the characteristic polynomial; distinct roots θ1=θ2 give A1θ1k+A2θ2k, a repeated root θ gives (A1+A2k)θk. For an inhomogeneous equation add a particular solution, raising the degree of the polynomial trial by one each time it is annihilated by the homogeneous operator.
Fix the constants from the boundary conditions, then sanity-check the limits: the symmetric case, the N→∞ case, and the cases k=0 and k=N should all give the obvious answers.
Example(A bounded betting run at the roulette wheel)Original
A player arrives at a European roulette table with 20 chips and decides to bet one chip on red on each spin, stopping when he reaches 50 chips or loses everything. A European wheel has 37 pockets, of which 18 are red, so
p=3718,q=3719,ρ=pq=1819=1.05.
Here k=20 and N=50.
Ruin probability. We need ρ20 and ρ50. From logρ=log(19/18)=0.0540672, we get 20logρ=1.0813444 and 50logρ=2.7033611, so
(The exact rational value is 481.2038, so the six-figure intermediate quantities have cost nothing.) The session lasts about 481 spins on average.
The fair game for comparison. Had the wheel had no zero, p=q=21 would give
p20=1−5020=0.6,D20=20×30=600.
The single green pocket, a house edge of 1/37=2.7% per spin, raises the chance of ruin from 0.600 to 0.860 and shortens the expected session from 600 spins to 481. A useful way to see where the money goes: the expected loss is q−p=37−1 chips per spin, so over 481.204 spins the player expects to lose 481.204/37=13.0055 chips. Independently, the expected final holding is 0×0.860110+50×0.139890=6.9945 chips, which is 20−13.0055. The two calculations agree to five figures, which checks the ruin probability and the mean duration against each other.
Application(Diffusion, polymers, and why chemistry is different in two dimensions)Original
Dimension decides whether diffusing particles find each other. Polya's trichotomy is the reason a diffusion-limited reaction behaves differently on a surface and in a bulk fluid. Two reactants performing independent lattice walks meet exactly when their difference walk, itself a random walk, visits the origin. In three dimensions the difference walk is transient, so a pair that fails to react on a first encounter is likely never to meet again, and the reaction rate is governed by fresh encounters: this is the classical Smoluchowski regime, with a rate constant that settles to a finite value. In two dimensions — a lipid membrane, an adsorbed layer, a thin film — the difference walk is recurrent, so a pair that has once come close will come close again and again, the same partners re-collide many times, and the effective rate constant acquires a logarithmic dependence on time and on system size. Membrane biochemistry and surface catalysis are quantitatively different from bulk chemistry for exactly this reason, and it is a statement about ∑n−d/2.
The same trichotomy controls exciton trapping and defect capture in crystals: a randomly walking excitation in a three-dimensional lattice has a positive chance of escaping a trap forever, while in a one- or two-dimensional conductor capture is certain and only the time matters.
Polymer physics uses the walk as the model of an ideal chain: the positions of the monomers of a flexible polymer in a theta solvent are a random walk with n the degree of polymerisation, giving the end-to-end distance ⟨∣Sn∣2⟩=nℓ for step length ℓ. The characteristic-function representation of the displacement after n independent steps of arbitrary length, used above for the lattice Green function and worked out for a general step density by Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §30.6, eq. (30.64), p. 1155), is precisely the tool used in the polymer literature to get the distribution of the end-to-end vector for a general segment distribution. What the excluded-volume interaction acts on is not recurrence but the number of self-contacts, and the two have different thresholds. The ideal chain in space is transient by part (d) — it comes back to its first monomer only with probability 0.3405 — and yet it touches itself without limit, because there are many pairs of monomers and each pair has its own chance of coinciding. The expected number of coincident pairs is ∑m≥1(n+1−2m)p0(2m), there being n+1−2m pairs of monomers 2m apart along the chain, and since ∑mp0(2m)=G−1 converges this grows linearly in n with slope G−1=0.5164: it is 5.80 at n=20, 39.9 at n=100 and 476 at n=1000, that is 0.29, 0.40 and 0.48 contacts per monomer, climbing towards 0.5164. Whether forbidding those contacts changes the exponent is settled by a different count, Flory's: the order of n2 monomer pairs each interact with strength v while confined to the volume Rd=nd/2 that the ideal chain occupies, so the dimensionless interaction parameter is
z∼Rdvn2=vn(4−d)/2,
which grows with n for d<4, is marginal at d=4 and dies away above. So the swollen exponent ν in ∣Sn∣∼nν replaces the ideal 21 for d<4 (ν=43 in two dimensions, ν≈0.588 in three), while at d=4 the ideal value survives up to logarithmic corrections. Both questions are governed by the same family of sums ∑n−d/2, but they weight it differently and their thresholds differ: d=2 for return to the origin, d=4 for excluded volume.
Gambler's ruin is the first-passage problem for an absorbing boundary, and reappears whenever a fluctuating quantity is watched against two thresholds: the survival of a small population against extinction and a carrying capacity, the sequential probability ratio test in statistics, the voltage of an integrate-and-fire neuron between reset and threshold, and the risk of ruin of an insurer. The formulas above are the discrete, two-barrier ancestors of the exit-time problems for diffusions solved in The Fokker-Planck equation, and the scaling limit that connects them is Brownian motion.
Gives the infinitesimal and interarrival definitions of the Poisson process and shows they describe the same process, derives the Poisson law of the counts and the exponential law of the gaps, proves the conditional property that fixes the number of points and leaves their positions an independent sample, and proves that superposing independent Poisson streams and independently thinning one both return Poisson processes with the obvious intensities.
A Poisson process is the mathematics of events that happen at a steady average rate and otherwise take no notice of one another: clicks on a Geiger counter, photons at a faint detector, calls at an exchange, raindrops on a tile. Three demands pin it down completely. Events are counted, never un-counted. In a short interval the chance of one event is proportional to the length of the interval, and the chance of two or more is smaller than that by an order. And what happens after a given time is unaffected by what happened before it.
Those three demands are stated below as the definition. Everything else in this unit is a consequence: the number of events in a window has a Poisson distribution, the gaps between events are independent and exponential, two independent streams merged give a third stream of the same kind, and a stream from which events are deleted by independent coin tosses is again of the same kind. The last two facts are what make the process so useful in practice, because real detectors miss particles and real signals are mixtures of sources.
Definition(The Poisson process, infinitesimally)AdaptedAdapted from Grimmett & Stirzaker (2020), §6.8, Definition (6.8.1), p. 277
Let λ>0. A Poisson process with intensity (or rate) λ is a family N={N(t):t≥0} of random variables taking values in {0,1,2,…} such that:
N(0)=0, and N(s)≤N(t) whenever s<t;
for every n≥0, as h↓0,
P(N(t+h)=n+mN(t)=n)=⎩⎨⎧λh+o(h),o(h),1−λh+o(h),m=1,m>1,m=0;
with the o(h) terms uniform in the base time t — which is what
time-homogeneity is meant to encode, and what lets the proof below re-read the relation at a shifted
base point rather than only pointwise in t.
whenever s<t, the count N(t)−N(s) of arrivals falling in (s,t] is independent of everything the process did on [0,s], the instants of those earlier arrivals included.
The three clauses are, in order: the process counts; it counts one at a time, at a rate that does not vary; and it has no memory. We call N(t) the number of arrivals (or events, or counts) by time t. The arrival times and interarrival times are
T0=0,Tn=inf{t:N(t)=n},Un=Tn+1−Tn(n≥0),
so that Tn=U0+U1+⋯+Un−1 and, in the other direction, N(t)=max{n:Tn≤t}. Each of the two descriptions therefore determines the other: to know the whole path of N is to know the whole sequence U0,U1,…, and conversely. The second description is the one used in simulation, and it is the subject of part 2 of the statement below.
Definition(Poisson processes on a region)AdaptedAdapted from Grimmett & Stirzaker (2020), §6.13, Definitions (6.13.1) and (6.13.4) and the remark at (6.13.2), pp. 319-320
The time axis plays no special role, and for superposition and thinning it is cleaner to count points in sets. Let d≥1 and let λ:Rd→R be a non-negative measurable function whose integral
Λ(A)=∫Aλ(x)dx,A∈Bd,
over the Borel σ-field Bd of Rd is finite for every bounded A. A random countable subset Π⊆Rd is a non-homogeneous Poisson process with intensity functionλ if the counts N(A)=∣Π∩A∣ satisfy:
N(A) has the Poisson distribution with parameter Λ(A), for every A∈Bd, where a Poisson distribution of parameter +∞ is to be read as the statement P(N(A)=∞)=1. That convention is the one Grimmett and Stirzaker print on p. 319 for their constant-intensity definition; we carry it over to a varying intensity function, where they do not repeat it. It is needed because the hypothesis on λ makes Λ finite on bounded sets only: an unbounded A may well have Λ(A)=∞, and then A holds infinitely many points almost surely;
N(A1),…,N(An) are independent whenever A1,…,An∈Bd are disjoint.
Λ is the mean measure of Π. The constant-intensity case λ(x)≡λ, in which Λ(A)=λ∣A∣ is λ times the volume of A, is the homogeneous Poisson process. Grimmett and Stirzaker record that in the constant-intensity case with d=1 the present definition and the previous one describe the same object: if Π is a Poisson process on R with constant intensity λ then M(t)=N([0,t]) satisfies the infinitesimal definition, and conversely every process satisfying the infinitesimal definition arises this way from some Π. For a genuinely varying intensity function no such comparison is available, since clause 2 of the infinitesimal definition fixes one rate λ for all times.
Definition(The counts and the gaps)AdaptedAdapted from Grimmett & Stirzaker (2020), §6.8, Theorem (6.8.2), p. 277, Theorem (6.8.10), p. 279, the converse remark and the calculation at the top of p. 280, and Theorem (6.8.11), p. 280
Let N be a Poisson process with intensity λ, with arrival times Tn and interarrival times Un as above.
The counts are Poisson. For each t≥0,
P(N(t)=j)=j!(λt)je−λt,j=0,1,2,…,
so EN(t)=varN(t)=λt.
The gaps are independent and exponential. The random variables U0,U1,U2,… are independent, each with the exponential distribution of parameter λ.
The converse. If U0,U1,… are independent, each exponential with parameter λ, and N is built from them by Tn=∑i<nUi and N(t)=max{n:Tn≤t}, then N is a Poisson process with intensity λ. Consequently Tn has the gamma distribution Γ(λ,n), with density λ(λx)n−1e−λx/(n−1)! on x>0, mean n/λ and variance n/λ2.
Conditional uniformity. Conditional on {N(t)=n}, the vector (T1,…,Tn) is distributed as the order statistics V(1)≤⋯≤V(n) of n independent random variables V1,…,Vn, each uniform on [0,t].
Statements 1-3 say that the infinitesimal description and the exponential-gaps description are two presentations of one process; statement 4 says that, once the number of arrivals in a window is fixed, their positions carry no further structure.
Definition(The conditional property)AdaptedAdapted from Grimmett & Stirzaker (2020), §6.13, Theorem (6.13.11) and the remark following it, p. 323
Take Π as in the region-based definition, with intensity function λ and mean measure Λ, and fix a Borel set A with 0<Λ(A)<∞. Write
Q(B)=Λ(A)Λ(B),B⊆A,
for the mean measure restricted to A and normalized to total mass 1; its density at a point x of A is λ(x)/Λ(A). Given that N(A)=n, the configuration Π∩A is then an independent sample of size n from Q. For a constant intensity Q is uniform on A, so in that case one conditions on the number of points and gets n positions thrown independently and uniformly into A.
Part 4 of the previous statement is the case d=1, A=[0,t], λ constant, rewritten in terms of the arrival times; it is proved separately below, directly from parts 1 and 2, so that the one-dimensional theory does not lean on the region-based definition. The general form stated here is what the thinning proof needs, because there the intensity is allowed to vary with position.
Definition(Superposition and thinning)AdaptedAdapted from Grimmett & Stirzaker (2020), §6.13, Superposition theorem (6.13.5) and Colouring theorem (6.13.14), pp. 320 and 324, and §6.8, Exercises (6.8.1) and (6.8.2), p. 287
Superposition. Take two independent Poisson processes Π′ and Π′′ on Rd, with intensity functions λ′ and λ′′ respectively. Their union Π=Π′∪Π′′ is then itself Poisson, with intensity function λ′+λ′′. The same argument covers countably many independent processes with mean measures Λi, provided ∑iΛi(A)<∞ for every bounded A.
Thinning (the colouring theorem). Now take a single Poisson process Π on Rd with intensity function λ. Give the point sitting at x the colour green, with probability γ(x), or else scarlet, with probability σ(x)=1−γ(x), and make that choice separately for each point, independently of every other choice. Writing Γ for the green points and Σ for the scarlet ones, the two sets Γ and Σ are independent Poisson processes, with intensity functions γ(x)λ(x) and σ(x)λ(x).
On the time axis these read: two independent Poisson streams of rates λ and μ merge into a Poisson stream of rate λ+μ; and if each arrival of a rate-λ stream is retained with probability p, independently of the others, the retained arrivals form a Poisson stream of rate λp, independent of the stream of discarded ones. Grimmett and Stirzaker set these two one-dimensional statements as exercises (flies and wasps landing on a plate; green insects in the soup) and prove the general versions as theorems.
Idea of the proofAdaptedAdapted from Grimmett & Stirzaker (2020), §6.8, the two solution methods for Theorem (6.8.2) and the proof of Theorem (6.8.10), pp. 277-279, and §6.13, the route of the proofs of (6.13.5) and (6.13.14), pp. 320 and 324; Riley, Hobson & Bence (2006), §30.8.2, the binomial limit on p. 1175
Why the counts are Poisson. Condition on the state a short time earlier. Clause 2 of the definition says that in (t,t+h] the process either stays put or moves up one, with probabilities 1−λh and λh to first order, whatever count it has reached by time t — the clause is stated conditionally on that count, which is all this first argument needs of the hypotheses. That turns pj(t)=P(N(t)=j) into a chain of differential equations, each coupling pj to pj−1 only. Packing them into the generating function G(s,t)=E(sN(t)) collapses the chain to one first-order equation in t, whose solution eλt(s−1) is recognisably the generating function of a Poisson variable with mean λt. Riley, Hobson and Bence reach the same formula from the other direction, by cutting [0,t] into n slices each carrying one event with probability p=λt/n and letting n→∞; the binomial limit and the differential-equation route are two views of the same first-order bookkeeping.
Why the gaps are exponential. No arrival in [0,t] is exactly the event {N(t)=0}, whose probability is e−λt by statement 1; so the first gap is exponential. For the later gaps, the point is that the event {U0=t0} is decided by what happens up to time t0 while {U1>t} is decided by what happens after it, and clause 3 makes these independent. Induction then peels off one gap at a time.
Why superposition and thinning work. Both are statements that a certain count is Poisson and that counts in disjoint sets are independent, so in both cases there are exactly two things to check. For superposition, the sum of two independent Poisson variables is Poisson with the summed parameters, and independence across disjoint sets is inherited from the two ingredients. For thinning, condition on the total count in a set: by the conditional property the points are then an independent sample from the normalized mean measure, so their colours are independent coin tosses and the counts of the two colours are jointly binomial; multiplying the binomial probability by the Poisson probability of the total makes the joint distribution factorise, and a product of two Poisson factors is exactly what independence of the two colours means.
Step 1: a differential-difference system. Split the event {N(t+h)=j} according to the value of N(t), which produces the conditional probabilities P(N(t+h)=j∣N(t)=i). Clause 2 of the definition is stated in exactly that conditional form, so it evaluates them directly, with no appeal to clause 3: 1−λh+o(h) for j=i, λh+o(h) for j=i+1, and o(h) for j>i+1. Since N is non-decreasing, no term with i>j contributes. Hence, for j≥1,
pj(t+h)=λhpj−1(t)+(1−λh)pj(t)+o(h),
and for j=0 the first term is absent, so p0(t+h)=(1−λh)p0(t)+o(h). Subtracting pj(t), dividing by h and letting h↓0 gives
pj′(t)=λpj−1(t)−λpj(t)(j≥1),p0′(t)=−λp0(t),
with pj(0)=δj0 because N(0)=0.
The limit just taken is one-sided, so what it delivers is the right-hand derivative of pj, whereas step 2 differentiates in t and solves an ordinary differential equation, which asks for the two-sided one. The gap closes because the error terms in clause 2 depend on the length h of the window and not on where the window sits — the rate λ is the same at all times, and this is the same reading of the clause that step 2 of the next proof uses. Granting that, the displayed relation holds with t replaced by t−h for any t>0 and small h>0:
pj(t)=λhpj−1(t−h)+(1−λh)pj(t−h)+o(h).
Read as it stands it gives pj(t)−pj(t−h)→0, so each pj is continuous from the left as well as from the right; divided by h it then gives λpj−1(t)−λpj(t) for the left-hand derivative too. So pj is differentiable on (0,∞) and right-differentiable at 0, which is what step 2 uses.
Step 2: collapse the system with a generating function. Put G(s,t)=∑j≥0pj(t)sj=E(sN(t)), which converges for ∣s∣≤1. Fix r<1 and work on the disc ∣s∣≤r. Multiply the jth equation of step 1 by sj and add over j≥0. Summing the left-hand sides gives ∂G/∂t, and this is the one step that needs an argument, since it differentiates a power series term by term: the equations themselves bound the derivatives, ∣pj′(t)∣≤λ(pj−1(t)+pj(t))≤2λ for every j and every t, so on ∣s∣≤r the differentiated series is dominated term by term by the convergent numerical series 2λ∑jrj, with a bound independent of t. The differentiated series therefore converges uniformly in t, which is exactly the hypothesis that licenses ∑jpj′(t)sj=∂G/∂t. On the right-hand sides, ∑j≥1λpj−1(t)sj equals λsG(s,t) and ∑j≥0λpj(t)sj equals λG(s,t), so
∂t∂G=λ(s−1)G,G(s,0)=1,
the boundary value being 1 because N(0)=0 almost surely. Since r<1 was arbitrary, this holds for every ∣s∣<1.
Step 3: solve and read off the distribution. For fixed s this is a linear first-order equation in t with constant coefficient λ(s−1), so G(s,t)=eλt(s−1). Expanding the exponential,
G(s,t)=e−λtj=0∑∞j!(λt)jsj,
and matching coefficients of sj against the definition of G gives pj(t)=(λt)je−λt/j!. The mean and variance of a Poisson variable with parameter λt are both λt (see Generating functions). □
Grimmett and Stirzaker also record the slower route: solve p0′=−λp0 to get p0(t)=e−λt, feed it into the equation for p1, and induct. Riley, Hobson and Bence derive the same system (their eq. (30.101)) from the switchboard picture and iterate it in that way.
Proof(2 and 3: the gaps, and the converse)AdaptedAdapted from Grimmett & Stirzaker (2020), §6.8, proof of Theorem (6.8.10), p. 279, and the derivation at the top of p. 280; Riley, Hobson & Bence (2006), §30.9.3, eqs. (30.116) and (30.118), pp. 1190-1191
Step 1: the first gap. The first arrival occurs after time t exactly when there is no arrival in [0,t], so by part 1,
P(U0>t)=P(N(t)=0)=e−λt,t≥0,
which is the exponential law with parameter λ.
Step 2: the second gap, given the first. Condition on U0=t0. By step 1 the conditioning event has probability zero, so what follows is a statement about conditional densities and the conditional probabilities computed from them, in the sense of Conditional expectation: we evaluate P(U1>t∣U0=t0) as a function of t0, and finding that it does not depend on t0 is what independence of U0 and U1 means. The event {U1>t} says that no arrival falls in (t0,t0+t]. The conditioning event is decided by the arrival pattern in [0,t0], and clause 3 of the definition makes the number of arrivals in (t0,t0+t] independent of that pattern. Hence
P(U1>t∣U0=t0)=P(no arrival in (t0,t0+t])=e−λt,
the last equality because the right-hand probability does not depend on where the window sits: clause 2 refers only to the length h of the window, so repeating step 1 on (t0,t0+t] gives e−λt. The answer does not involve t0, so U1 is independent of U0 and has the same law.
Step 3: induction. Suppose U0=t0,…,Un−1=tn−1 — again a null event, and again read through the conditional densities of Conditional expectation — and set T=t0+⋯+tn−1. Then {Un>t} is the event of no arrival in (T,T+t], which by the same two reasons has probability e−λt regardless of the conditioning values. So Un is exponential with parameter λ and independent of U0,…,Un−1, which completes the induction. □
Step 4: the converse. Suppose instead that we are given independent exponential gaps and define N by N(t)=max{n:Tn≤t} with Tn=∑i<nUi. Then Tn, a sum of n independent exponential variables of parameter λ, has the gamma distribution Γ(λ,n); this is the law Riley, Hobson and Bence obtain as their eq. (30.118) by asking for r−1 events in [0,x] and one more in [x,x+dx]. The construction gives the identity
{N(t)≥j}={Tj≤t},
since having at least j arrivals by time t is precisely the jth arrival occurring by then, whence
P(N(t)=j)=P(Tj≤t)−P(Tj+1≤t)=j!(λt)je−λt,
by integrating the two gamma densities by parts.
Step 5: the converse, clause 1.T1=U0>0 almost surely, so no arrival has occurred at time 0 and N(0)=0; and {n:Tn≤t} grows with t, so N is non-decreasing. That N(t) is finite needs Tn→∞: by the strong law of large numbers (Laws of large numbers) Tn/n→1/λ almost surely, and 1/λ>0, so Tn→∞ and only finitely many arrivals fall in [0,t].
Step 6: the converse, the restart at a fixed time. The one fact about the exponential law we need is that it forgets: for u,v≥0,
P(U>u+v∣U>u)=e−λue−λ(u+v)=e−λv,
so the unexpired part of a gap that has already lasted u is again exponential with parameter λ. Now fix t≥0 and condition on the entire arrival pattern on [0,t]: on N(t)=n and on T1∈dt1,…,Tn∈dtn with tn≤t. As at steps 2 and 3, this conditioning event has probability zero, and the dtk are the usual shorthand: the statement being proved is one about the conditional law of the gaps that follow time t given the vector (T1,…,Tn) and the event {N(t)=n}, read through the conditional densities of Conditional expectation, and its content is that this conditional law does not depend on n or on the values t1,…,tn. Written in terms of the gaps, the conditioning event is
It pins down U0,…,Un−1, constrains Un only through the one-sided event {Un>t−tn}, and says nothing whatever about Un+1,Un+2,…. Since the gaps are independent, two things follow: Un+1,Un+2,… are untouched by the conditioning and remain independent exponentials of parameter λ; and, by the lack-of-memory calculation just made, the residual
R=Tn+1−t=Un−(t−tn)
is exponential with parameter λ and independent of them. So, conditionally on the pattern on [0,t], the sequence R,Un+1,Un+2,… is a sequence of independent exponential variables of parameter λ.
Step 7: the converse, clauses 3 and 2. The arrivals after time t are at t+R, t+R+Un+1, and so on, so the shifted process s↦N(t+s)−N(t) is built from the sequence R,Un+1,… by exactly the recipe that built N from U0,U1,…. By step 6 that sequence has the same law whatever the conditioning pattern was, so the shifted process has the law of N and is independent of the arrival pattern on [0,t]. That is clause 3, and it also shows that the law of an increment depends only on the length of its window. Clause 2 follows: averaging over the patterns compatible with N(t)=n,
P(N(t+h)=n+m∣N(t)=n)=P(N(h)=m)=m!(λh)me−λh,
by step 4. Expanding the exponential gives 1−λh+o(h) for m=0 and λh+o(h) for m=1, while the whole tail m≥2 carries probability 1−e−λh−λhe−λh=21(λh)2+O(h3)=o(h). All three lines of clause 2 hold, so N is a Poisson process with intensity λ. □
Proof(4: conditional uniformity, and the conditional property on a region)AdaptedAdapted from Grimmett & Stirzaker (2020), §6.8, proof of Theorem (6.8.11), pp. 280-281, and §6.13, proof of Theorem (6.13.11), the partition calculation (6.13.13), pp. 323-324
Part 4: the arrival times, given their number. Fix n≥1 and t>0.
Step 1: the joint density of the first n arrival times. By part 2 the gaps U0,…,Un−1 are independent exponentials of parameter λ, so their joint density is λne−λ(u0+⋯+un−1) on (0,∞)n. Now apply the multivariate density-transformation rule of Joint distributions and transformations: under a smooth bijective change of variables, the new density is the old one composed with the inverse map and divided by the modulus of the Jacobian determinant. The map (u0,…,un−1)↦(t1,…,tn) given by tk=u0+⋯+uk−1 is linear and triangular with 1s on the diagonal, so its Jacobian determinant is 1; it carries (0,∞)n bijectively onto the ordered region {0<t1<⋯<tn}, and it turns the sum in the exponent into tn. Hence T=(T1,…,Tn) has density
fT(t1,…,tn)=λne−λtnon {0<t1<⋯<tn}.
Step 2: impose the condition N(t)=n. Let C be a Borel subset of {0<t1<⋯<tn≤t}. Given the values of T, the event {N(t)=n} says precisely that the next gap outlasts the time remaining, Un>t−tn; and Un is independent of U0,…,Un−1, hence of T, so that event has conditional probability e−λ(t−tn) given T=(t1,…,tn). Integrating that conditional probability against the density of T, which is how a probability is recovered from a conditional probability given a continuous vector (Conditional expectation),
The conditional density of T is the constant n!/tn on the ordered region, and 0 off it. Note that λ has disappeared.
Step 4: recognize the answer. Let V1,…,Vn be independent and uniform on [0,t]; their joint density is t−n on [0,t]n. The ordered vector (V(1),…,V(n)) lives on the ordered region, and each point of that region receives the mass of the n! rearrangements of its coordinates, each contributing t−n; so the order statistics have joint density n!/tn there. That is the density found in step 3, which proves part 4. □
The conditional property on a region. Let 0<Λ(A)<∞.
Step 1: count over a partition. Let A1,…,Ak be disjoint Borel sets with union A, and let n1+⋯+nk=n. The counts N(Ai) are independent, each Poisson with parameter Λ(Ai), and their sum is N(A), which is Poisson with parameter Λ(A)=∑iΛ(Ai). Dividing the joint mass function by that of the sum,
the exponential factors cancelling because the Λ(Ai) add up to Λ(A).
Step 2: read the right-hand side. It is the multinomial probability that n points thrown into A independently, each with law Q, land n1 in A1, ..., nk in Ak.
Step 3: a classical input. What remains is to pass from the partition probabilities of step 1 to the law of the whole configuration, and generating a σ-field is by itself not enough for that: two probability measures can agree on a generating family and differ on the σ-field it generates. The standard repair is to work with a family closed under finite intersections.
Quoted classical input: Dynkin's π-λ theorem. Call a family C of events a π-system if E,F∈C implies E∩F∈C. Then (i) two probability measures that agree on a π-system agree on the σ-field it generates, and (ii) if C1,…,Ck are π-systems, each containing Ω, with P(E1∩⋯∩Ek)=∏iP(Ei) for every choice of Ei∈Ci, then the σ-fields σ(C1),…,σ(Ck) are independent. The condition Ω∈Ci is not decoration. Without it the single full-intersection identity does not give independence once k≥3: on Ω={1,…,8} with the uniform measure, the one-element families C1={A}, C2={B}, C3={C} with A={1,2,3,4}, B={1,2,3,5}, C={1,5,6,7} are π-systems satisfying P(A∩B∩C)=81=P(A)P(B)P(C), yet P(A∩B)=83=41, so the three σ-fields are not independent. With Ω available, the identity for any subcollection follows from the full one by padding the missing slots with Ω, and that is what the proof of independence needs. This theorem is proved in no unit of this handbook and is not in the sources followed here; it is used, in both parts, as a quoted classical result, in the way the σ-field machinery of Sigma-algebras and measures is used generally.
Step 4: a π-system of count events. For a Borel set A write
(a) It is a π-system, by construction: an intersection of two finite intersections of events {N(B)=m} is again a finite intersection of such events, and so is again in C(A) — including when it is empty, as for {N(B)=0}∩{N(B)=1}. It also contains Ω, as the case r=1, B1=∅, m1=0, since N(∅)=0 identically; this is the hypothesis part (ii) above needs.
Arbitrary Borel sets are used here, and not the cells of one partition, for a definite reason. The cylinders of step 1 — events {N(B1)=n1,…,N(Br)=nr} in which B1,…,Br are the cells of a single finite partition of A — are the natural family to try, and it is tempting to argue that they are closed under intersection because two partitions have a common refinement. That argument is wrong. Take two partitions {B1,B2} and {B1′,B2′} of A in general position, so that all four refinement cells Bj∩Bk′ carry positive mean measure. Then
{N(B1)=1,N(B2)=1}∩{N(B1′)=1,N(B2′)=1}
is the union of exactly two cylinders over the refinement {Bj∩Bk′} — one point in B1∩B1′ and one in B2∩B2′, or one point in B1∩B2′ and one in B2∩B1′ — and both have positive probability, so the intersection is a union of two distinct refinement cylinders and is therefore not one of them. Passing to a common refinement makes the two generated σ-fields sit inside a third; it does not turn an intersection of cylinders into a cylinder. Rather than try to repair that route we avoid it: C(A) is closed under intersection by its very definition, and by (b) it generates the same σ-field.
(b) It generates the σ-field of the configuration. That σ-field is σ(N(B):B⊆A Borel), and each N(B) takes values in {0,1,2,…}∪{∞}, so it is generated by the events {N(B)=m} alone, every one of which lies in C(A).
(c) The two laws agree on C(A). Fix B1,…,Br⊆A and m1,…,mr, and let P be the finite Borel partition of A that they generate: the nonempty sets B1±∩⋯∩Br±, where Bj+=Bj and Bj−=A∖Bj. Each Bj is the union of the cells of P it contains, so N(Bj)=∑D⊆BjN(D). Under either law all n points lie in A, so
j=1⋂r{N(Bj)=mj}=⨄{N(D)=nD for all D∈P},
a disjoint union over the count vectors (nD)D∈P with ∑DnD=n and ∑D⊆BjnD=mj for each j. Every event in that union is of the step-1 form, and step 1 gives it the same probability under the conditioned process as under the independent Q-sample; summing the finitely many terms gives the same probability for the left-hand side.
By (a), (b), (c) and part (i) of the quoted theorem, the conditioned process and the independent Q-sample have the same law. □
In the constant-intensity case Q(B)=∣B∣/∣A∣ is the uniform distribution on A, which for d=1 and A=[0,t] returns part 4 once the sample is written in increasing order.
Proof(Superposition and thinning)AdaptedAdapted from Grimmett & Stirzaker (2020), §6.13, proof of the Superposition theorem (6.13.5), p. 320, and proof of the Colouring theorem (6.13.14), p. 324; Riley, Hobson & Bence (2006), p. 1179, additivity of independent Poisson variables
Superposition. Fix A∈Bd and put N′(A)=∣Π′∩A∣, N′′(A)=∣Π′′∩A∣. These are independent, and Poisson with parameters Λ′(A) and Λ′′(A).
Step 1: the count in one set. The sum of two independent Poisson variables is Poisson with the sum of the parameters — the generating functions multiply, eΛ′(A)(s−1)eΛ′′(A)(s−1)=e(Λ′(A)+Λ′′(A))(s−1) — so S(A)=N′(A)+N′′(A) is Poisson with parameter Λ′(A)+Λ′′(A)=∫A(λ′+λ′′). This is the fact Riley, Hobson and Bence state for sums of independent Poisson variables.
Step 2: independence across disjoint sets. If A1,A2,… are disjoint then the families {N′(Ai)} and {N′′(Ai)} each consist of independent variables, and the two families are independent of each other; so the sums S(Ai) are independent.
Step 3: the union really has S(A) points.∣Π∩A∣ equals S(A) unless some point of Π′ coincides with a point of Π′′, in which case the union loses a point. Grimmett and Stirzaker show that almost surely no such coincidence occurs, by covering a bounded A with the O(2nd) dyadic boxes B of side 2−n that meet it and bounding the chance that one box holds a point of each. Since Π′ and Π′′ are independent and P(N′(B)≥1)≤Λ′(B), that chance is at most Λ′(B)Λ′′(B). It is the sum over the boxes that has to go to zero, and it does, because the sum factorises against its largest factor:
where A is the union of the boxes. The second factor is a fixed finite number, while the first tends to 0 as n→∞ by absolute continuity of the integral Λ′′(B)=∫Bλ′′ over sets of volume 2−nd→0. So the chance of a coincidence anywhere in A is zero. They describe this step as technical and set it aside on a first reading. □
Thinning. Fix A with Λ(A)<∞, and let Ng, Ns count the green and scarlet points of Π in A. If Λ(A)=0 then N(A) is Poisson with parameter 0, so N(A)=0 and Ng=Ns=0 almost surely and the conclusion below holds trivially, with both colour parameters equal to 0; so assume 0<Λ(A)<∞, which is the range in which the conditional property was proved.
Step 1: condition on the total. Given ∣Π∩A∣=n, the conditional property proved above says these n points are distributed as n independent points of A, each drawn from the probability measure Q(B)=Λ(B)/Λ(A). This is the general form of the property, with the intensity allowed to vary over A; the uniform special case of part 4 would not be enough here.
Step 2: the colours are then independent coin tosses. Because the points are independent and each is coloured independently of the others, the conditional probability that a given one is green is the position-averaged colouring probability
γˉ=∫Aγ(x)dQ=Λ(A)1∫Aγ(x)λ(x)dx,
a single number, to be kept apart from the function γ(⋅) of the statement; write σˉ=∫Aσ(x)dQ=1−γˉ for its scarlet counterpart. The colours of the n points are independent, so, conditionally on the total being n,
P(Ng=g,Ns=sN(A)=n)=g!s!n!γˉgσˉs,g+s=n.
Step 3: multiply by the Poisson weight and watch it factorise. With n=g+s,
using γˉ+σˉ=1 to split e−Λ(A). The joint mass function is a product of two Poisson mass functions, so Ng and Ns are independent and Poisson, with parameters
γˉΛ(A)=∫Aγ(x)λ(x)dx,σˉΛ(A)=∫Aσ(x)λ(x)dx,
which are the mean measures of the intensity functions γ(x)λ(x) and σ(x)λ(x) claimed in the statement.
Step 4: disjoint sets. Let A1,…,Ak be disjoint, each of finite mean measure. The first thing to establish is that the whole configurations Π∩A1,…,Π∩Ak, and not merely the counts N(A1),…,N(Ak), are independent across i; clause 2 of the region-based definition says only the latter. The passage from the one to the other is the π-system argument set out at steps 3 and 4 of the proof of the conditional property above, and it is worth doing in full, because the obvious family of events is the wrong one.
Take for each i the family C(Ai) of finite intersections {N(B1)=m1}∩⋯∩{N(Br)=mr} with B1,…,Brarbitrary Borel subsets of Ai. By (a) and (b) of step 4 of the conditional-property proof, with Ai in place of A, it is a π-system and it generates σ(Π∩Ai). The cylinders over the cells of one fixed finite partition of Ai are not the family to use: as the two-partition example given there shows, an intersection of two of them is in general a union of refinement cylinders rather than a cylinder, so the common-refinement argument does not supply the π-system hypothesis and the quoted theorem would be applied without it.
Now fix Ei∈C(Ai) for each i, say Ei={N(Bi1)=mi1}∩⋯∩{N(Biri)=miri}, and let Pi be the finite Borel partition of the union Bi1∪⋯∪Biri generated by those sets: the nonempty cells Bi1±∩⋯∩Biri±, where B+=B and B−=Ai∖B, with the all-minus cell discarded. It is the partition of the union of the Bij, and not of Ai, that is wanted here. In (c) of the conditional-property proof the partition was taken over the whole of A, and there the conditioning ∑DnD=n held the leftover cell's count down; in the present step the total count over Ai is not conditioned, so the leftover cell Ai∖(Bi1∪⋯∪Biri) would carry a count on which Ei imposes no constraint whatever, and the decomposition below would run over infinitely many terms. Dropping that cell removes the difficulty. Each Bij is a union of cells of Pi, so N(Bij)=∑D∈Pi,D⊆BijN(D), and splitting Ei over the values of the cell counts gives
Ei=c⨄c,c={N(D)=nD for all D∈Pi},
a disjoint union over the vectors (nD)D∈Pi of non-negative integers satisfying ∑D⊆BijnD=mij for every j≤ri. Every cell of Pi lies in at least one Bij, by construction, so each nD is bounded by maxjmij and only finitely many vectors qualify: this union is genuinely finite. The cells of P1,…,Pk taken together are pairwise disjoint, because the Ai are, so clause 2 makes the whole collection of counts {N(D):D∈P1∪⋯∪Pk} independent. Hence the probability of an intersection c1∩⋯∩ck of one cylinder from each Ei is ∏iP(ci), and expanding each Ei and multiplying out,
where ci runs over the cylinders making up Ei; the middle step is the distributive law, legitimate because each Ei is a union of finitely many cylinders, and the outer ones use that the cylinders making up a given Ei are disjoint. So the kπ-systems are independent, and part (ii) of the quoted π-λ theorem makes the σ-fields σ(Π∩A1),…,σ(Π∩Ak) independent. That is the independence of the configurations.
Now the colour of each point is decided by a toss made independently of every other toss and of all the positions. The pair (Ng(Ai),Ns(Ai)) is a function of the configuration Π∩Ai and of the tosses attached to its points, so the pairs are independent of one another. That is more than independence of the green counts among themselves and of the scarlet counts among themselves, and it is what the statement needs: combined with step 3, which makes Ng(Ai) and Ns(Ai) independent within each pair, it says that the whole family {Ng(Ai)}∪{Ns(Ai)} is independent.
Step 5: sets of infinite mean measure. Clause 1 of the region-based definition quantifies over every Borel set, while steps 1-4 kept the mean measure finite on the sets in play; the sets with Λ(A)=∞, on which the green or the scarlet mean measure may be finite or infinite, are reached from the bounded ones by exhaustion. Let A be an arbitrary Borel set and put Am=A∩{x:∣x∣≤m}, so Λ(Am)<∞ because Λ is finite on bounded sets, and Am↑A. Write Λγ(B)=∫Bγλ. Then Ng(Am)↑Ng(A) pointwise, and Λγ(Am)↑Λγ(A) by monotone convergence. If Λγ(A)<∞, the Poisson mass functions of parameter Λγ(Am) converge term by term to those of parameter Λγ(A), so the laws of the Ng(Am) converge to the Poisson law of that parameter; an almost-sure limit has the limiting law, so Ng(A) is Poisson with parameter Λγ(A). If Λγ(A)=∞, then for each fixed k,
P(Ng(A)≤k)≤P(Ng(Am)≤k)⟶0(m→∞),
so Ng(A)=∞ almost surely, which is clause 1 read with the infinite-parameter convention. The same argument serves for Ns. The independence statements of steps 3 and 4 survive the limit as well: disjoint Borel sets A1,…,Ak have disjoint bounded pieces A1∩{∣x∣≤m},…,Ak∩{∣x∣≤m}, for which those steps apply, and independence is preserved under almost-sure limits. Hence Γ and Σ each satisfy both clauses of the region-based definition, with the intensity functions found in step 3, and the two processes are independent of each other. □
Remark(Where the hypotheses are used, and what fails without them)AdaptedAdapted from Grimmett & Stirzaker (2020), §6.8, Definition (6.8.15) and Theorem (6.8.17), pp. 281-282, and Example (6.8.28), p. 286, and §6.13, Definition (6.13.1) with the convention for an infinite mean measure, p. 319, Definition (6.13.4) and the Superposition theorem (6.13.5), p. 320, and the Colouring theorem (6.13.14), p. 324
Clause 3 (independence of the future from the past) is what makes the gaps independent. It is used twice, both times in the proof of part 2: at step 2 and at step 3, where in each case the conditioning event is a statement about the past and the event being evaluated is a statement about the future. Part 1 does not need it at all — clause 2 is itself stated conditionally on the present count, so it evaluates the conditional probabilities thrown up by splitting {N(t+h)=j} over the value of N(t) without any further hypothesis. In the converse it is not available and has to be recovered, which is the work of steps 6 and 7. Its region-based counterpart is clause 2 of the definition on a region, independence of the counts in disjoint sets, and that is what carries the independence through the superposition and thinning proofs.
What fails without it is step 2 of the proof of part 2, and with it the conclusion. Clauses 1 and 2 speak only about the next instant given the present count; neither says anything about how the increment after time t depends on when the earlier arrivals occurred. The counting processes with unit jumps and no simultaneous arrivals are a far larger family than the Poisson processes — every renewal process whose gap distribution is continuous and not exponential belongs to it — and what picks the Poisson process out of that family is exactly the demand that the past be forgotten. Such a renewal process violates clause 2 as well, since its instantaneous arrival rate depends on how long the current gap has already run; that is the lack-of-memory property of step 6 read backwards, as a characterization of the exponential law. So it is not a counterexample in the strict sense, and we offer none: clauses 1 and 2 speak only about the next instant conditional on the current count, never about the arrival times behind that count, and whether they can be met while clause 3 fails is a question this unit does not settle. What the renewal family does establish is the weaker point the proofs need — the conclusions above are not consequences of the unit-jump structure of clause 1 alone, and clause 3 is genuinely used at the two places just listed.
The constancy of λ in clause 2 is what makes the gaps identically distributed. If the rate is allowed to depend on the number of arrivals so far, λ↦λn, the definition becomes that of a birth process, and the gaps are still independent but exponential with different parameters λn. Then the counts need not be Poisson, and worse can happen: Grimmett and Stirzaker's Theorem (6.8.17) says that such a process, with all its birth rates strictly positive, stays finite for all time if and only if ∑nλn−1=∞, so for rapidly growing rates infinitely many arrivals occur in a bounded time. This is the point at which the Poisson process stops being the whole story and Continuous-time Markov chains take over.
The o(h) in the second line of clause 2 forbids simultaneous arrivals, and it is what makes N increase by one at a time. Consider a process that jumps by 2 whenever a Poisson process of rate λ has an arrival. It satisfies clause 1 and clause 3, and it satisfies the m=0 line of clause 2, since the chance of no jump in (t,t+h] is e−λh=1−λh+o(h). It fails the other two lines, and in opposite directions: the chance of exactly one further arrival is 0, not λh+o(h), and the chance of two is λh+o(h), not o(h). The second failure is the one at issue here: the o(h) line is by itself enough to rule the process out, whether or not the m=1 line is also consulted. How much the clause buys can be read off Grimmett and Stirzaker's (Grimmett & Stirzaker (2020), §6.8, Example (6.8.28), p. 286) Example (6.8.28). Suppose all that is known of a process M is this: it starts at the origin; its paths never decrease and are right-continuous; every value it takes is a whole number; the increments it makes over disjoint time intervals are independent of one another, and the law of an increment depends on the interval only through its length; and it never moves up by more than one at a time. Those demands already force M to be a Poisson process. Unit jumps together with stationary independent increments is therefore a characterization, not merely a convenient hypothesis.
For the superposition theorem, independence of the two processes is essential, and the coincidence step is not decoration. Superposing a Poisson process with a translated copy of itself gives a point set whose counts in disjoint sets are strongly dependent. The coincidence step is the one place where the fact that the mean measures have no atoms — guaranteed here by Λ(A)=∫Aλ, an integral of a function — is used; a "Poisson process" allowed to put positive mass on a single point would overlap its independent twin with positive probability, and the union would then have fewer points than the sum of the counts.
For thinning, the colours must be chosen independently of one another and of the positions, apart from the dependence on position allowed through γ(x). Step 2 fails otherwise. A detector whose dead time makes it more likely to miss a particle shortly after recording one does not produce a Poisson stream of recorded particles, however nearly Poisson the incident stream is; the recorded counts are under-dispersed relative to the Poisson law, which is why var/E is the standard diagnostic.
Finiteness of Λ on bounded sets keeps the counts in bounded regions finite, and that is all it does. It does not keep Λ(A) finite for unbounded A, which is why clause 1 of the region-based definition carries the convention Grimmett and Stirzaker state on p. 319 for a constant intensity: a Poisson distribution of infinite parameter means P(N(A)=∞)=1. The process is still defined on such an A; it simply has infinitely many points there. Every result above that conditions on N(A), the conditional property and the thinning proof included, assumes Λ(A)<∞ explicitly for just this reason; the thinning conclusion is then carried to the sets of infinite mean measure by exhausting them with bounded ones, which is step 5 of that proof.
Method(Modelling a stream of events as a Poisson process)Original
Check the three clauses against the physics, not against the data. Events must be countable and non-simultaneous; the instantaneous rate must not depend on the past (no dead time, no refractory period, no queueing feedback); and the rate must be constant in time. If the rate varies with time but nothing else fails, keep going with the non-homogeneous version, replacing λt by Λ(A)=∫Aλ(x)dx everywhere below.
Fix the intensity λ in units of events per unit time. If several independent sources contribute, add their rates: that is superposition.
Apply any thinning. If each event is recorded with probability p, independently of the others, replace λ by λp; the recorded stream is again Poisson, and is independent of the missed stream.
Answer count questions with the Poisson law. The number of events in a window of length t is Poisson with parameter λt; for a non-homogeneous process, with parameter Λ of that window.
Answer waiting-time questions with the gamma law. The time to the nth event is Γ(λ,n), with mean n/λ and variance n/λ2; for n=1 this is the exponential law. Convert between the two languages with {N(t)≥n}={Tn≤t}, which turns any statement about waiting times into one about counts.
Answer "where did they land" questions with conditional uniformity. Given the number of events in a window, their positions are an independent uniform sample of that window.
Test the model with the index of dispersion. For Poisson counts, varN(t)=EN(t). A sample variance well above the sample mean points to clustering (the rate itself fluctuates); well below it points to a dead time or some other inhibition.
Example(A two-source counting experiment with an imperfect detector)Original
Two weak radioactive sources sit in front of one detector. Source A emits particles towards the detector at a mean rate of 1.2s−1, source B at 0.8s−1, and the two emission streams are independent Poisson processes. The detector registers each incident particle with probability p=0.6, independently of everything else. Let R be the process of registered counts.
Step 1: superpose. The incident particles form the union of two independent Poisson processes, so by the superposition theorem they form a Poisson process of rate
λinc=1.2+0.8=2.0s−1.
Step 2: thin. Colour an incident particle "registered" with probability 0.6 and "missed" with probability 0.4, independently of the others. By the colouring theorem R is a Poisson process of rate
λ=pλinc=0.6×2.0=1.2s−1,
independent of the process of missed particles.
Step 3: a count. The number of registered counts in a 3s window is Poisson with parameter λt=1.2×3=3.6, so
P(R(3)=2)=2!3.62e−3.6=6.48e−3.6=0.17706.
Step 4: a waiting time. The time T5 to the fifth registered count is Γ(1.2,5), with
ET5=1.25=4.1667s,sd(T5)=1.25=1.8634s.
Step 5: a waiting time answered as a count. What is the chance that we have to wait more than 4s for the third registered count? Use {T3>4}={R(4)≤2}. Here R(4) is Poisson with parameter 1.2×4=4.8, so
P(T3>4)=e−4.8(1+4.8+24.82)=17.32e−4.8=0.142539.
That is the number to check: P(T3>4s)=0.1425 to four decimal places.
Step 6: a consistency check on the model. The missed particles form an independent Poisson process of rate 0.4×2.0=0.8s−1, and superposing the registered and missed streams must return the incident stream: 1.2+0.8=2.0s−1, as it does. Also, had we thinned first and superposed afterwards — thinning each source separately to rates 0.72 and 0.48 — we would have reached the same λ=1.2s−1. The two operations commute, as they must, since both are read off from the intensity function alone.
Application(Shot noise, dead time, and counting statistics)AdaptedAdapted from Grimmett & Stirzaker (2020), §6.13, the shot-noise sum over the points of a Poisson process, p. 326
Every counting instrument in the physical sciences is calibrated against this process.
Photon and particle counting. A photomultiplier viewing a steady faint source, a Geiger counter beside a long-lived isotope, a cosmic-ray telescope: in each case the arrival stream is Poisson to good accuracy, and the characteristic consequence is that a measurement of n counts carries an uncertainty of about n, because varN(t)=EN(t)=λt. The familiar rule that one must quadruple the integration time to halve the fractional error, λt/(λt)=(λt)−1/2, is this identity and nothing more.
Detector efficiency is thinning. A detector of quantum efficiency p turns an incident Poisson stream of rate λ into a recorded Poisson stream of rate pλ. This is why photon-counting statistics are not degraded in form by an inefficient detector — only the rate falls, and the recorded counts remain Poisson. It is also why departures from Poisson statistics in recorded data are informative: they cannot have been manufactured by a memoryless loss mechanism, so they come either from the source (bunching, clustering, a varying rate) or from a detector with memory.
Dead time is the standard failure of the model. After each count, a real counter is insensitive for a short interval τ. The retained events are then not an independent thinning, since whether a particle is recorded depends on how recently the last one was; the recorded process is a renewal process whose gaps are τ plus an exponential variable. Those gaps are less variable than exponential ones of the same mean, and the counts they produce are under-dispersed: varR(t) falls below ER(t). Measuring an index of dispersion below 1 is the usual way to detect and then correct for dead time.
Shot noise in a circuit. If each arrival at rate λ contributes a pulse r(t−Ti) of fixed shape to a current, the total current is I(t)=∑ir(t−Ti), a sum over the points of a Poisson process; Grimmett and Stirzaker call the one-dimensional version of such sums shot noise. The superposition theorem says that two independent sources of such pulses give a current of the same type with the rates added, which is why shot noise from independent channels adds in power rather than in amplitude.
Spatial versions. The region-based definition is used directly for the positions of impurities in a crystal, of stars in a survey field, or of raindrops on a surface. There the thinning theorem is the statement that a magnitude-limited sample from a Poisson field of stars is still a Poisson field, with intensity reduced by the detection probability at each position — the usual justification for modelling a flux-limited catalogue as a Poisson process with a position-dependent intensity.
Defines a continuous-time Markov chain through its transition semigroup and its generator, shows that the chain is an exponential holding time in each state followed by a jump governed by a discrete chain, derives the Kolmogorov forward and backward equations and the solution by a matrix exponential, and gives the stationary condition and the birth-death criterion.
Note(What changes when time becomes continuous)Original
A discrete-time chain is run by one matrix: the transition matrix P, with Pn giving the n-step probabilities (Discrete-time Markov chains). In continuous time there is no smallest step, so there is no single matrix to iterate. There is instead a whole family {Pt:t≥0}, one matrix per elapsed time, and the family multiplies: Ps+t=PsPt. A family of matrices that multiplies like that is governed by its derivative at t=0, and that derivative is the generatorG. Everything in this unit is an unpacking of the slogan Pt=etG and of the places where it needs care.
The generator also has a direct probabilistic reading, and it is worth having both pictures from the start. Off the diagonal, gij is the rate at which the chain jumps from i to j. On the diagonal, −gii is the total rate of leaving i, which is the parameter of the exponential time the chain waits there. So a continuous-time chain is a discrete-time chain (where it goes) run on an exponential clock (when it goes): a map of where it can go, together with a speed for each place.
Definition(The chain, its semigroup and its generator)AdaptedAdapted from Grimmett & Stirzaker (2020), §6.9, Definitions (6.9.1), (6.9.2), (6.9.3) and the display (6.9.7), pp. 288-290
Let S be a countable set, which we take to be a set of integers. A family X={X(t):t≥0} of random variables with values in Ssatisfies the Markov property if
for all states and all increasing sequences t1<t2<⋯<tn of times. A sample path of X is right-continuous if for each t it is constant on some interval [t,t+ϵ); equivalently, since S is discrete, the paths are step functions. We call X a continuous-time Markov chain when both hold: its paths are right-continuous, and the Markov property is satisfied.
Its transition probabilities are pij(s,t)=P(X(t)=j∣X(s)=i) for s≤t, and X is homogeneous if these depend only on t−s; we assume homogeneity throughout and write pij(t), and Pt for the matrix (pij(t)). The family {Pt:t≥0} is the transition semigroup of the chain.
Suppose X(0)=i, and write U0=inf{t:X(t)=i} for the instant at which the chain first moves away from i; right-continuity forces U0>0. Write gi∈[0,∞) for the parameter of the exponential distribution of U0 (that it is exponential is part 3 of the statement below), and define the jump matrixY=(yij) by yij=δij if gi=0, and by
yij=Pi(X(U0)=j)if gi>0.
The generator of the chain is the matrix G=(gij:i,j∈S) with
gij={giyij,−gi,j=i,j=i.
Since Y has row sums 1, G has row sums 0: ∑jgij=0 for every i, equivalently ∑j=igij=gi. More generally, any matrix with non-negative off-diagonal entries and row sums 0 is called a generator on S. The off-diagonal entries are rates, in the sense that for j=i
P(X(t+h)=j∣X(t)=i)=gijh+o(h),
together with pii(h)=1+giih+o(h) on the diagonal. Grimmett and Stirzaker state this as their equation (6.9.9), set the finite case as an exercise, and warn that it must be handled carefully when S is infinite; it is part 1 of the Kolmogorov statement below and is proved there, since everything in that statement rests on it.
Let Tm be the instant of the mth move (T0=0), put Um=Tm+1−Tm, and set T∞=limmTm. We call the Um the holding times. The successive values taken by X form a discrete-time chain Y={Yn}, governed by the matrix Y; this is the jump chain. Following the source, the single letter Y carries both meanings — the matrix Y=(yij) and the chain Y={Yn} it runs — and we keep that economy here; the subscripts tell the two apart. The chain explodes from i if Pi(T∞<∞)>0, and one says explosion occurs if it explodes from some state. When explosion is possible one adjoins a state ∞ to S and sets X(t)=∞ for t≥T∞; the process so defined is the minimal process, because it runs as a chain in S for the shortest possible time.
Definition(The architecture of a continuous-time chain)AdaptedAdapted from Grimmett & Stirzaker (2020), §6.9, Theorem (6.9.4), p. 289, Theorem (6.9.5) and Theorem (6.9.6), p. 290, Theorem (6.9.11), p. 291, the construction (6.9.12) and Proposition (6.9.14), p. 292, and Theorem (6.9.17), p. 293
Let X be a homogeneous continuous-time Markov chain with state space S.
The Markov property extends. Fix t>0. Conditioning on {X(t)=j} together with any event decided by the chain strictly before t assigns any event decided strictly after t the same probability as conditioning on {X(t)=j} alone: this is the extended Markov property. In its strong form, for a stopping time T and conditional on {T<T∞}∩{X(T)=i}, the shifted process X∗(u)=X(T+u) is again a continuous-time Markov chain with the transition probabilities of X, started at i and independent of {X(s):s<T}. Grimmett and Stirzaker state both and omit the proof, referring to Norris; we do the same, and flag below the three places where the extended form rather than the bare Markov property is what is being used.
The semigroup property.P0=I; each Pt has non-negative entries and row sums 1; and the Chapman-Kolmogorov equations hold,
Ps+t=PsPt,s,t≥0.
(For the minimal process built from a generator this still holds, because the chain can reach the adjoined state ∞ but cannot come back.)
Holding times are exponential. Started at X(0)=i, the variable U0 defined above is exponentially distributed, with some parameter gi∈[0,∞), and is a stopping time. The value gi=0 means Pi(U0=∞)=1, that is, i is absorbing; the value gi=∞ is impossible, because right-continuity forces Pi(U0>0)=1.
Destination and duration are independent, and the chain restarts. If gi>0 then X(U0) is independent of U0, and, conditional on X(U0)=j, the shifted process X∗(s)=X(U0+s) is again a chain of the same kind, carrying the transition probabilities of X, begun at j, and independent both of U0 and of the starting state i.
The chain is exactly this description. Conversely, given the two ingredients of the construction — a transition matrix Y with yii=0, together with non-negative constants gi, which is the same data as a generator matrix G in the sense defined above — the recipe "sit in Y0 for an exponential time of parameter gY0, jump to Y1, sit for an exponential time of parameter gY1, and so on" defines a process, which is a Markov chain provided it does not exhaust S in finite time.
Sufficient conditions for non-explosion. The chain does not explode from i if any one of the following holds: S is finite; supjgj<∞; or i is recurrent for the jump chain Y. None of these is necessary: a pure birth process with rates λn=n, started from a state i≥1, has an infinite state space, unbounded rates and a transient jump chain, yet it never explodes, because ∑n≥1n−1 diverges. (Started from 0 it is a different matter, since λ0=0 makes 0 absorbing, and then condition (c) holds there anyway.)
The chain is irreducible if for every pair i,j there is some t>0 with pij(t)>0; in that case pij(t)>0 for allt>0, a dichotomy with no analogue in discrete time, where periodicity makes some n-step probabilities vanish.
Definition(The Kolmogorov equations)AdaptedAdapted from Grimmett & Stirzaker (2020), §6.9, equation (6.9.9), p. 291 and Exercise (6.9.11), p. 296, and §6.10, equations (6.10.1)-(6.10.5), pp. 297-298, and Theorem (6.10.7), p. 299
Let X be a chain as above, with generator G and semigroup {Pt}, and let Pt′ denote the matrix whose entries are the derivatives pij′(t).
The short-time expansion. If S is finite then, as h↓0,
pij(h)=gijh+o(h)(j=i),pii(h)=1+giih+o(h),
the error terms being uniform in i and j; in matrix form Ph=I+hG+o(h), so that G=P0′ is the derivative of the semigroup at the origin. This is the statement that makes the generator a matrix of rates, and everything below depends on it.
Finite state space. If S is finite then {Pt} satisfies the forward equationPt′=PtG,t≥0,
and the backward equationPt′=GPt,t≥0,
both subject to P0=I. Each is a finite system of first-order linear differential equations with constant coefficients, so each has a unique solution, namely
Pt=n=0∑∞n!tnGn=etG.
Countable state space. If S is countable and X is the minimal chain with generator G, then, among all non-negative solutions of the backward system with P0=I, the family {Pt} is the smallest; and the same is true for the forward system. ("Minimal" means pij(t)≤πij(t) for any other non-negative solution {Πt} with the same boundary condition.) Any solution carrying that boundary value has row sums at most 1, that is ∑jpij(t)≤1. When explosion does not occur the row sums equal 1, and then {Pt} is the unique solution of either system.
Both equations come from the same two-step decomposition Pt+h=PtPh=PhPt; the forward equation conditions on the state at time t and looks at the last short interval, the backward equation conditions on the state at time h and looks at the first.
Definition(Stationary distributions, the limit theorem, and birth-death chains)AdaptedAdapted from Grimmett & Stirzaker (2020), §6.9, the definitions (6.9.21) of recurrence and transience and of positive recurrence, with Theorem (6.9.22), p. 295, §6.10, Definition (6.10.13), equation (6.10.14), Theorem (6.10.15) and Theorem (6.10.22), pp. 301-304, and §6.11, conditions (2), equation (3), conditions (4) and (6) and Theorem (6.11.7), pp. 307-308
Evolution of the distribution. If X(0) has distribution μ(0), written as a row vector, then X(t) has distribution μ(t)=μ(0)Pt.
Recurrence. Grimmett and Stirzaker call a state irecurrent for X when Pi({t:X(t)=i} is unbounded)=1, and transient when that probability is 0; when gi>0, the state i is recurrent for X precisely when the jump chain is recurrent at i, and transient for X precisely when the jump chain is transient there. A recurrent i is positive recurrent if either S={i} or the first return time Ri=inf{t>U0:X(t)=i} has finite mean mi=Ei(Ri).
Stationarity. For an irreducible, non-explosive X, a measure π is a stationary measure if π=πPt for all t≥0, and a stationary distribution if moreover ∑jπj=1. Now let X be irreducible with ∣S∣≥2.
Should one state be positive recurrent, a stationary distribution π exists and is unique, it is the only distribution satisfying πG=0, and every state is then positive recurrent too.
Suppose instead that X does not explode and that πG=0 for some distribution π. Then every state is positive recurrent, π is stationary, and πk=1/(mkgk).
So the infinite system π=πPt, one equation for every t, reduces to the single linear system πG=0.
The limit theorem. Assume once more that X is irreducible and does not explode. Where a stationary distribution π exists, it is unique and pij(t)→πj as t→∞ for every pair i,j. Where none exists, pij(t)→0 for every pair. No aperiodicity hypothesis appears, in contrast with the discrete-time theorem of Convergence of Markov chains.
Birth-death chains. A birth-death process is a chain on {0,1,2,…} whose generator has gn,n+1=λn, gn,n−1=μn, gnn=−(λn+μn) and all other off-diagonal entries zero, where λi≥0, μi≥0 and μ0=0. Assume
λi>0(i≥0),μi>0(i≥1),μ0=0,
so that the chain is irreducible. Then:
the measure
xn=μ1μ2⋯μnλ0λ1⋯λn−1x0,n≥1,
satisfies xG=0;
a distributionπ with πG=0 exists if and only if
n=0∑∞μ1μ2⋯μnλ0λ1⋯λn−1<∞
(the n=0 term being read as 1), and then x0 is the reciprocal of that sum;
that π is a stationary distribution of X if, in addition,
n=0∑∞μ1μ2⋯μnλ0λ1⋯λn−1(λn+μn)<∞.
The extra condition 3 is not idle: Grimmett and Stirzaker exhibit rates for which πG=0 has a normalisable solution while the chain is transient, and therefore explosive.
Idea of the proofOriginal
Why there is a generator at all. The semigroup property Ps+t=PsPt with P0=I is the matrix version of the functional equation f(s+t)=f(s)f(t), whose solutions are exponentials. Differentiating the semigroup property at t=0 produces a single matrix G=P0′, and the same differentiation at a general t turns the family into a differential equation driven by G. Which side G appears on depends on which of the two factorisations Pt+h=PtPh or Pt+h=PhPt one differentiates, and that is the only difference between the forward and the backward equations.
Why the holding times are exponential. Starting afresh in state i at time s, given that the chain has not yet left i, is the same as starting in i at time 0. The conditioning event {U0>s} is a statement about the whole of [0,s] and not just about X(s), so what is used here is the extended Markov property of part 1 of the architecture statement, together with homogeneity; the same remark applies to part 4 of that statement, which needs the strong form at the stopping time U0. So Q(t)=Pi(U0>t) satisfies Q(s+t)=Q(s)Q(t), the same functional equation again, now for a scalar; with Q(0)=1, and its only right-continuous solutions are Q(t)=e−git for some gi∈[0,∞).
Why πG=0 replaces πP=π. Stationarity should mean that the distribution does not move, so dtd(πPt)=0. By the backward equation this derivative is πGPt, and since P0=I it vanishes for every t exactly when πG=0 — one linear system in place of a continuum of them.
Why the birth-death measure is a product. A birth-death generator is tridiagonal, so the equation xG=0 links only consecutive components and can be unwound one step at a time; what comes out is the detailed-balance relation xnλn=xn+1μn+1 — the mean rate of upward crossings of the gap between n and n+1 must equal the rate of downward crossings — and iterating that relation gives the product formula.
Where the care is needed. Every step above interchanges a limit with a sum. For finite S the sums are finite and nothing can go wrong. For infinite S the chain may run through infinitely many states in finite time, and then the forward and backward systems may have many solutions; the transition probabilities of the chain are the smallest of them.
Proof(The forward and backward equations, and the matrix exponential)AdaptedAdapted from Grimmett & Stirzaker (2020), §6.9, equation (6.9.9) and Exercise (6.9.11), pp. 291 and 296, and §6.10, the derivation of (6.10.1)-(6.10.5), pp. 297-298; Riley, Hobson & Bence (2006), §30.8.2, eq. (30.101), p. 1176
Take S finite, and fix i,j∈S. Write g=maxk∈Sgk, which is finite because S is finite.
Step 1: the short-time expansion, part 1 of the statement. Start the chain at i. By parts 3, 4 and 5 of the architecture statement the chain sits at i for an exponential time U0 of parameter gi, moves to a state j=i with probability yij, and then restarts there; so all three probabilities below can be computed from the holding times and the jump chain.
(a) Two moves in a short time are unlikely, and this is where finiteness enters first. Conditional on the jump chain, U0 and U1 are independent exponential variables with parameters gi and gY1, both at most g; since P(U>t)=e−g∙t≥e−gt for each of them, the pair is stochastically larger than a pair of independent exponential variables of parameter g. The sum of the latter has the gamma law Γ(g,2) (The Poisson process, part 3), so
Pi(T2≤h)≤1−e−gh−ghe−gh=21g2h2+O(h3)=O(h2),
with a bound depending only on g, hence uniform over the states.
(b) The diagonal entry. Either the chain has not moved by time h, which has probability e−gih, or it has moved at least once. If it has moved and is nevertheless at i at time h, it has moved at least twice: the first move cannot end at i, because right-continuity makes the path a step function and U0 is the first instant at which the value differs from i. So pii(h) differs from e−gih by at most Pi(T2≤h)=O(h2), and since e−gih=1−gih+O(h2) and gii=−gi,
pii(h)=1+giih+O(h2).
(c) The off-diagonal entries. Fix j=i and split {X(h)=j} by the number of moves made by time h. The part with two or more moves contributes at most O(h2) by (a). The part with exactly one move requires that move to occur at some instant s≤h, to land on j, and to be followed by no further move in the remaining time h−s; integrating over s, and using the independence of destination from duration and the restart at j,
The two exponentials together lie between e−gh and 1, so the integral lies between gihe−gh and gih, and therefore equals gih+O(h2). Multiplying by yij and recalling gij=giyij,
pij(h)=gijh+O(h2).
All the error bounds are controlled by the single constant g, so they are uniform in i and j, which is what the next two steps need. □
Step 2: condition on where the chain is at time t. By the Chapman-Kolmogorov equations,
pij(t+h)=k∈S∑pik(t)pkj(h).
Step 3: expand the short-time factor. Step 1 gives pkj(h)=gkjh+o(h) for k=j and pjj(h)=1+gjjh+o(h). Substituting,
The o(h) terms may be collected into one because there are finitely many of them, each bounded by the uniform error of step 1. This is the second place finiteness is used, and it is exactly the step Grimmett and Stirzaker flag as needing justification when S is infinite, where the sum is infinite and a limit is being taken through it.
Step 4: divide and let h↓0. We obtain
pij′(t)=k∈S∑pik(t)gkj=(PtG)ij,
which is the forward equation Pt′=PtG. Conditioning instead on the state at time h, that is writing pij(t+h)=∑kpik(h)pkj(t) and expanding the first factor in the same way, gives the backward equation Pt′=GPt. In both cases P0=I because pij(0)=δij.
Step 5: solve. Each system is a linear constant-coefficient system of ∣S∣2 scalar equations, so by the existence and uniqueness theorem for linear systems it has exactly one solution with the given value at t=0. Define Et=∑n≥0tnGn/n!, the matrix exponential of tG; the series converges for every t because ∥Gn∥≤∥G∥n in any submultiplicative norm. Differentiating term by term, which is legitimate because the differentiated series converges uniformly on compact t-intervals by the same bound,
dtdEt=n≥1∑(n−1)!tn−1Gn=EtG=GEt,
the last equality because G commutes with every power of itself. Also E0=I. So Et solves both systems with the right boundary value, and by uniqueness
Pt=etG.□
Step 6: what survives when S is infinite. Steps 1 and 3 break, and with them uniqueness. Grimmett and Stirzaker's Theorem (6.10.7) salvages the following: the transition semigroup of the minimal chain is still a non-negative solution of both systems, and it is the smallest one. Their proof of the backward half starts from the decomposition
obtained by splitting according to whether the first jump happens after t or at some time s≤t; substituting u=t−s, exchanging sum and integral by Fubini's theorem, and differentiating recovers pij′(t)=∑kgikpkj(t). Minimality then follows by showing, by induction on n, that any non-negative solution dominates Pi(X(t)=j,Tn>t), and letting n→∞.
Proof(The stationary measure of a birth-death chain)AdaptedAdapted from Grimmett & Stirzaker (2020), §6.11, the derivation of (3) from (2) and the footnote on detailed balance, p. 307, with the discussion of condition (6) and Example (6.11.8), p. 308
Let G be a birth-death generator satisfying λi>0 for i≥0, μi>0 for i≥1 and μ0=0, and let x be a measure with xG=0. Reading off the columns of G, the equation xG=0 is
Step 1: recognise a telescoping quantity. Define Dn=λnxn−μn+1xn+1 for n≥0. The nth equation above says precisely Dn−1=Dn for n≥1, and the n=0 equation says D0=0. Hence Dn=0 for every n:
λnxn=μn+1xn+1,n≥0.
This is the detailed-balance relation; it holds here not by assumption but because the generator is tridiagonal, so that the only way probability can cross the gap between {0,…,n} and {n+1,n+2,…} is by a single step.
Step 2: unwind it. Since μn+1>0 we may divide, and induction on n gives
Conversely this x satisfies both displayed equations, as substituting back shows, so it is the general solution of xG=0 up to the scalar x0.
Step 3: normalise.x can be scaled to a distribution exactly when ∑nxn<∞, which with x0=1 is the summability condition of the statement; and then x0 is the reciprocal of the sum. □
Why the second condition is needed for stationarity. Having a normalisable solution of πG=0 is not the same as having a stationary distribution, because part 2 of the stationarity statement also requires the chain to be non-explosive. The route Grimmett and Stirzaker take is to pass to the jump chain: the measure νn=xngn=xn(λn+μn) is stationary for the jump chain Y, and if it is summable then Y is recurrent, whence X is recurrent and so non-explosive. That summability is the second condition. Their Example (6.11.8), with λi=λgi and μi=μgi where λ+μ=1, taking gn=2n and λ/μ∈(1,2), satisfies the first condition but not the second; there the jump chain is transient, so X is transient, and since a normalisable π with πG=0 exists, X must be explosive.
Remark(Where the hypotheses are used, and the master equation)AdaptedAdapted from Grimmett & Stirzaker (2020), §6.9, Theorem (6.9.6) and the footnote on (6.9.9), pp. 290-291, Theorem (6.9.17), p. 293, and Example (6.9.16), p. 293, and §6.10, equations (6.10.1)-(6.10.3) and Theorem (6.10.22), pp. 297-298 and 304
Right-continuity is used in three places. It makes the jump times Tm well defined, so that the holding times, the jump chain and the construction of part 5 of the architecture statement all make sense. It forces Pi(U0>0)=1 and hence gi<∞ in part 3 of that statement, which we take from Grimmett and Stirzaker rather than prove: if the paths were not right-continuous the chain could leave i immediately, U0 would be 0 almost surely, and gi would be ∞. And it is what makes the path a step function, so that the first move out of i cannot end at i — the fact used in step 1(b) of the proof to show that leaving i and returning within a window of length h costs O(h2). Grimmett and Stirzaker adopt right-continuity precisely to avoid the measure-theoretic complications of the general theory, and note that without it one must look elsewhere for a treatment.
Finiteness of S is used twice in the proof. In step 1 it bounds the jump rates, g=maxkgk<∞, which is what makes two moves in a window of length h cost O(h2) and makes the errors in the short-time expansion uniform over the states. In step 3 it makes the sum over k finite, so that finitely many uniform o(h) terms can be collected into one and the limit taken through the sum. Everything downstream of those steps — uniqueness of the solution, Pt=etG, the equivalence of the forward and backward systems — is conditional on it. For infinite S the two systems can have different solution sets, and the chain supplies the minimal one.
Explosion is the obstruction, and it is real. For a birth process with rates λn, Grimmett and Stirzaker's criterion is that the chain is non-explosive if and only if ∑nλn−1=∞. So birth rates λn=(1+n)2 give a population that becomes infinite in finite time with positive probability, and then ∑jpij(t)<1, the deficit being the probability of having escaped to infinity by time t. The three sufficient conditions in part 6 of the architecture statement are the ones worth checking in practice, and the first — S finite — covers most applications.
Irreducibility and non-explosion are both needed in the limit theorem, and nothing else is. Aperiodicity, which the discrete-time theorem cannot do without, has no analogue here: by the dichotomy noted above, an irreducible continuous-time chain has pij(t)>0 for every t>0, so there is no parity obstruction for a limit to overcome. A chain that merely satisfies πG=0 with π a distribution, but explodes, is not covered, and the birth-death example above shows the conclusion can then fail.
The master equation. Physicists usually meet the forward equation not as a matrix identity but as a balance of gains against losses. Let μj(t)=P(X(t)=j). Multiplying Pt′=PtG on the left by the initial distribution gives dtdμ(t)=μ(t)G, that is μj′(t)=∑iμi(t)gij. Splitting off the diagonal term and using gjj=−∑k=jgjk,
dtdμj=flow into ji=j∑μi(t)gij−flow out of jμj(t)k=j∑gjk.
This gain-loss form is what is called the master equation of the process; the rates gij are then usually written W(i→j) or Wji. Stationarity, πG=0, says that the two flows balance for every state when summed over the other states; detailed balance, πigij=πjgji for each pair, is the strictly stronger statement that they balance pair by pair, and for a birth-death chain the two coincide, because the only pairs joined by a direct transition are neighbours and step 1 of the proof above showed that the neighbour balances then hold automatically. Note also that the diagonal of G carries no information beyond the off-diagonal entries: it is fixed by the row-sum condition, which is why a master equation is specified by its rates alone.
Method(Setting up and solving a continuous-time chain)Original
List the states and the transitions. For each ordered pair i=j of states between which a single event can carry the system, write down the rate gij at which that event occurs when the system is in i. Rates, not probabilities: gij has dimensions of inverse time.
Assemble the generator. Put the gij off the diagonal, zeros where no direct transition exists, and fill the diagonal with gii=−∑j=igij so that every row sums to zero. Check that every off-diagonal entry is non-negative; a sign error here is the commonest mistake.
Read off the holding times and the jump chain if you want the sample paths. The chain stays in i for an exponential time of mean 1/gi with gi=−gii, and then moves to j with probability yij=gij/gi. This is also the recipe for simulating the chain.
Check for explosion if S is infinite. If S is finite, or supigi<∞, or the jump chain is recurrent, there is nothing to do. Otherwise the transition probabilities are the minimal solution of the Kolmogorov equations and may not sum to 1.
For the equilibrium behaviour, solve the linear system πG=0 together with ∑jπj=1. For a birth-death chain do not solve it as a system: use the product formula, then normalise. If the chain is irreducible and non-explosive and the system has a solution, that solution is the unique stationary distribution and pij(t)→πj.
For the time dependence with S finite, exponentiate. If G has a full set of eigenvectors, diagonalise G=BΛB−1; then Pt=BetΛB−1. A generator need not be diagonalisable, and when it is defective there are two standard ways on: put G in Jordan form, which replaces etθ by etθ times a polynomial in t of degree one less than the block size, or invert the Laplace transform, Pt=L−1[(σI−G)−1], which needs no eigenvectors at all (The matrix exponential). Either way the spectrum is the same, and it always contains the eigenvalue 0, since G has row sums 0 and therefore the right eigenvector (1,1,…,1)T. For a finite irreducible generator the eigenvalue 0 is simple and all the other eigenvalues have negative real parts, and the slowest of them, −∣Reθ∣ with θ the non-zero eigenvalue closest to the imaginary axis, sets the relaxation time 1/∣Reθ∣ towards equilibrium.
For an infinite S with structure, use generating functions. Multiply the master equation for μn(t) by sn and sum, turning the infinite system into one partial differential equation for G(s,t)=∑nμn(t)sn, exactly as for the Poisson process in The Poisson process.
Example(Two instruments, two technicians)Original
A laboratory runs two identical instruments. Each working instrument fails at rate φ=21 per hour; each failed instrument is repaired at rate r=23 per hour, and there are two technicians, so repairs proceed in parallel. All four clocks are independent and exponential. Let X(t)∈{0,1,2} be the number of failed instruments.
Step 1: the generator. From state 0 (both working) either instrument may fail, so the total failure rate is 2φ=1. From state 1 the one working instrument fails at rate φ=21 and the one failed instrument is repaired at rate r=23. From state 2 both repairs are under way, at total rate 2r=3. This is a birth-death chain with λ0=1, λ1=21, μ1=23, μ2=3, so
G=−12301−23021−3,
whose rows sum to 0 and whose off-diagonal entries are non-negative.
Step 2: the stationary distribution, by the product formula. The birth-death recursion xn=(λn−1/μn)xn−1 terminates here after two steps, since the state space is finite:
At t=1 hour, starting from both instruments working,
p00(1)=0.5625+0.375e−2+0.0625e−4=0.6143954586,
and p02(1)=161(1−e−2)2=0.04672782. As t→∞ these tend to 169 and 161, as the limit theorem requires. The slowest non-zero eigenvalue is −2, so the relaxation time is 21 hour.
Step 4: an independent check. Because the two instruments never interact — each has its own clock and its own technician — X(t) is the sum of two independent two-state chains. A single instrument has generator (−1/23/21/2−3/2), hence stays up with probability 43+41e−2t and has stationary up-probability 43. Therefore X(t) should be binomial, and indeed
(43+41e−2t)2=169+83e−2t+161e−4t=p00(t),
while π is the binomial distribution with two trials and success probability 41: (169,166,161). The eigenvalues 0,−2,−4 are 0 and −2 from one factor added to 0 and −2 from the other.
Step 5: a number the laboratory cares about. In equilibrium both instruments are down a fraction π2=161 of the time, that is 24/16=1.5 hours per day. The mean number down is 2×41=0.5. Finally, the mean time to reach the state "both down" starting from "both up" solves the hitting-time system 1+∑kgjkmk=0 for j∈/{2}, with m2=0; this is the system of Grimmett and Stirzaker's Exercise (6.9.6) (Grimmett & Stirzaker (2020), Exercise (6.9.6), p. 296), whose minimal non-negative solution gives the mean hitting times. Here it reads
1−m0+m1=0,1+23m0−2m1=0,
whose solution is m1=5 and m0=6. Starting from a fully working laboratory, the expected wait until both instruments are simultaneously down is exactly 6 hours.
Application(Master equations in physics, chemistry and queueing)AdaptedAdapted from Grimmett & Stirzaker (2020), §6.9, Exercise (6.9.3), p. 296, and §6.11, Example (6.11.12) and Theorem (6.11.14), pp. 308-309, Example (6.11.15) with Theorem (6.11.16), pp. 309-310, and Corollary (6.11.18), p. 310
Rate equations for populations of states. Whenever a physical system occupies one of a discrete set of configurations and hops between them at rates fixed by the configuration alone, its probability distribution obeys a master equation, and the system is a continuous-time Markov chain. Nuclear decay chains, the populations of atomic levels under incoherent pumping and spontaneous emission, the open and closed states of an ion channel, and the conformations of a macromolecule are all of this kind. In each case the physics supplies the generator, the equilibrium populations follow from πG=0, and the relaxation rates are the non-zero eigenvalues of G. Where the rates satisfy detailed balance with respect to a Boltzmann distribution, gije−βEi=gjie−βEj, that Boltzmann distribution is automatically stationary — which is how a kinetic model is made thermodynamically consistent.
Chemical kinetics with small numbers of molecules. In a well-stirred vessel containing few copies of each species, the state is the vector of molecule counts and each reaction channel contributes a transition at a rate proportional to the number of available reactant combinations. The resulting master equation is the chemical master equation; the deterministic rate law is its large-volume limit. A single species that is produced at constant rate λ and destroyed independently per molecule at rate μ is exactly a birth-death chain with λn=λ and μn=nμ, and Grimmett and Stirzaker's Theorem (6.11.14) gives the equilibrium: the molecule number is Poisson with mean ρ=λ/μ. The fluctuation is therefore ρ, so the relative noise in a species held at a copy number of a hundred is about ten per cent — the basic reason that stochastic effects matter in gene expression.
Population growth and extinction. The simple birth-death chain with λn=nλ and μn=nμ models a population in which individuals reproduce and die independently. Its generating function is computed in Grimmett and Stirzaker's Theorem (6.11.16); differentiating that generating function gives the mean population EI(X(t))=Ie(λ−μ)t, which they display just after the theorem rather than inside it, and their Corollary (6.11.18) gives an extinction probability tending to 1 if λ≤μ and to (μ/λ)I if λ>μ — the continuous-time counterpart of the extinction criterion for branching processes.
Queues. Jobs arriving at a processor in the manner of a Poisson process of rate λ, each requiring an independent exponential service time of rate μ, make the number in the system a birth-death chain with λn=λ for all n and μn=μ for n≥1. The product formula gives xn=(λ/μ)n, which is summable exactly when λ<μ; Grimmett and Stirzaker set this as Exercise (6.9.3), with the conclusion that a stationary distribution exists if and only if λ<μ, in which case πn=(1−ρ)ρn with ρ=λ/μ. The mean number in the system is then ρ/(1−ρ), which blows up as the arrival rate approaches the service rate: the familiar fact that a server loaded to ninety per cent of capacity carries, on average, nine jobs rather than one.
Defines strong and weak stationarity and the autocovariance function, proves that an autocorrelation function is non-negative definite and therefore, when continuous at the origin, is the characteristic function of a distribution on the frequency axis, and identifies that distribution's density as the spectral density of the process.
A stationary process has no preferred instant: the statistics you measure do not depend on when you start the clock. That single invariance has a surprisingly rigid consequence. The correlation between the process now and the process a lag t later is a function c(t) of the lag alone, and c cannot be an arbitrary function: it must be the Fourier transform of a positive measure.
That measure is the spectral distribution of the process, and where it has a density, that density is the spectral density. It answers the question an experimentalist actually asks of a noisy signal — how much of the variance sits at each frequency — and it is the object a spectrum analyser estimates. The theorem below, usually named after Wiener and Khinchin, is the statement that the question always has an answer and that the answer is determined by the autocorrelation function alone.
Three warnings before we start. First, there is a second, purely deterministic theorem of the same name, about the Fourier transform of the correlation of two square-integrable functions; it is stated in the remark below so that the two are not confused. Second, several incompatible normalisations of the spectral density are in use; the conversion table is also in the remark. Third, the second name is a transliteration from Russian and is spelt differently by our two sources — Khinchin by Grimmett and Stirzaker, Kinchin by Riley, Hobson and Bence, Khintchine elsewhere again. We write it as the source we are reporting writes it, so both spellings appear below; they are the same man and, where the theorems differ, the difference is in the theorems and not in the name.
Definition(Stationarity)AdaptedAdapted from Grimmett & Stirzaker (2020), §8.2, Definitions (8.2.1) and (8.2.2) with (8.2.3) and the discussion after them, pp. 407-408, with Example (8.2.5), p. 408, and §9.1, Definitions (9.1.1)-(9.1.3), p. 434
Let X={X(t)} be a random process, with t ranging over [0,∞), over R, or over the integers.
X is strongly stationary if, for every n, all times t1,…,tn and every h>0, the families
(X(t1),…,X(tn))and(X(t1+h),…,X(tn+h))
have the same joint distribution. In particular X(t) then has the same distribution for every t.
X is weakly stationary — the alternative names are second-order stationary and covariance stationary — when, for every t1, every t2 and every h>0,
Equivalently, the mean is constant and the autocovariance c(t,t+h)=cov(X(t),X(t+h)) depends only on h.
Neither notion contains the other. A strongly stationary process without a finite second moment need not be weakly stationary, and Grimmett and Stirzaker's Example (8.2.5) — Xn=Acos(λn)+Bsin(λn), where λ∈[0,π] is fixed and A, B are uncorrelated with mean 0 and variance 1 — is weakly stationary but, without further conditions on the joint law of A and B, not strongly so. Following them, we read "stationary", unqualified, as "weakly stationary" from here on, because weak stationarity is both the easier hypothesis to verify and the one the spectral theory needs.
We allow X to take values in C, which costs nothing and is the natural setting. For complex-valued random variables the covariance is defined with a conjugate in the second slot,
so that cov is not symmetric: cov(C2,C1)=cov(C1,C2). Two complex variables with zero covariance are called orthogonal rather than uncorrelated.
Definition(Autocovariance, autocorrelation, spectral distribution and density)AdaptedAdapted from Grimmett & Stirzaker (2020), §9.3, the opening of the section with (9.3.1), Definition (9.3.3), Definition (9.3.9) and Definition (9.3.10), pp. 438-441
Let X={X(t):−∞<t<∞} be a weakly stationary complex-valued process. Its autocovariance function is
c(t)=cov(X(s),X(s+t)),
which by stationarity does not depend on s. Taking t=0 shows that the variance is the same at every time: var(X(t))=c(0). Replacing X by X−EX(0) changes nothing in c, so we may and do assume EX(t)=0 whenever it is convenient.
If c(0)>0, the autocorrelation function is the normalised autocovariance
the correlation between the process and its own value a lag t later. Note ρ(0)=1.
If ρ admits a representation
ρ(t)=∫−∞∞eitλdF(λ)
with F a distribution function, then F is the spectral distribution function of X, and where F has a density, that density f=F′ is the spectral density function. The spectrum of X is the set of real λ with F(λ+ϵ)−F(λ−ϵ)>0 for every ϵ>0: the frequencies that actually carry weight.
The reading of this is that ρ is a mixture of pure oscillations eitλ, the weight attached to frequencies in [λ,λ+dλ] being dF(λ); equivalently, ρ(t)=E(eitΛ) for a random frequency Λ with distribution function F. Because ρ is normalised, so is F: the spectral distribution describes how the fraction of the variance is shared out among frequencies. The variance itself carried by a band B of frequencies is c(0)∫BdF, and where F has a density f the variance per unit frequency — the quantity a spectrum analyser plots — is c(0)f(λ).
Theorem(The autocovariance function is Hermitian and non-negative definite)AdaptedAdapted from Grimmett & Stirzaker (2020), §9.3, Theorem (9.3.2), p. 438, and §9.6, Theorem (9.6.1), p. 464
Let X be weakly stationary with autocovariance function c. Then
c(−t)=c(t) for all t; and
c is non-negative definite: for all real t1,…,tn and all complex z1,…,zn,
j,k∑c(tk−tj)zjzk≥0.
These two properties characterise autocovariance functions. For a real-valued function c:R→R with c(−t)=c(t) and c non-negative definite, Grimmett and Stirzaker construct a real-valued strongly stationary process with autocovariance c, by declaring its finite-dimensional distributions to be multivariate normal with zero means and covariance matrix (c(tk−tj)) and appealing to the Kolmogorov consistency conditions. So nothing is lost and nothing is gained by asking which functions can be autocovariances: exactly these.
Theorem(The spectral theorem for autocorrelation functions)AdaptedAdapted from Grimmett & Stirzaker (2020), §9.3, Theorem (9.3.4), equations (9.3.11)-(9.3.14) and Theorem (9.3.15) with (9.3.16)-(9.3.18), pp. 439-442
Continuous time. Let X={X(t):−∞<t<∞} be a weakly stationary process with strictly positive variance and autocorrelation function ρ. Suppose ρ is continuous at the single point t=0. Then there is a distribution function F whose characteristic function is ρ:
ρ(t)=∫−∞∞eitλdF(λ),t∈R.
This is the result usually named after Wiener and Khinchin. The distribution function F is unique, by the uniqueness theorem for characteristic functions.
Discrete time. If X={Xn} is indexed by the integers, ρ is defined only on Z and cannot be a characteristic function until its domain is extended; the theorem survives in the form
ρ(n)=∫(−π,π]einλdF(λ),n∈Z,
for a distribution function F with F(−π)=0 and F(π)=1; no continuity hypothesis is appropriate here. The restriction to one period is possible because einλ is unchanged by λ↦λ+2π: if ρ on Z is represented by any distribution on the whole of R, the mass of that distribution may be folded into (−π,π] by shifting it through multiples of 2π, and the integrals at integer n are unchanged. If moreover X is real-valued then ρ(n)=ρ(−n) and the representation may be taken in the real symmetric form
ρ(n)=∫[−π,π]cos(nλ)dG(λ)
for a distribution function G of a symmetric distribution on [−π,π]; and conversely, a function ρ with ρ(0)=1 admitting such a representation is the autocorrelation function of some real-valued stationary sequence.
Inversion. Grimmett and Stirzaker's Theorem (9.3.15) states: if the function F above is differentiable with derivative f, then
f(λ)=2π1n=−∞∑∞e−inλρ(n)
and they attach to this the qualification that it holds wherever f is itself differentiable. For a real-valued sequence the formula may be written f(λ)=(2π)−1∑nρ(n)cos(nλ). They record ∑n∣ρ(n)∣<∞ as a sufficient, but not necessary, condition for the spectral density to exist.
Proposition(A clean sufficient condition in discrete time)AdaptedAdapted from Grimmett & Stirzaker (2020), §9.3, the inversion formula (9.3.16) and the sufficient condition (9.3.18), p. 442
Let X={Xn:n∈Z} be weakly stationary with c(0)>0 and autocorrelation function ρ, and suppose
n=−∞∑∞∣ρ(n)∣<∞.
Then the series
f(λ)=2π1n=−∞∑∞ρ(n)e−inλ
converges absolutely and uniformly on (−π,π], its sum f is continuous and non-negative, and it is the spectral density of X:
ρ(n)=∫−ππeinλf(λ)dλ,n∈Z,∫−ππf(λ)dλ=1.
The formula and the summability hypothesis are Grimmett and Stirzaker's (9.3.16) and (9.3.18); what is ours is the packaging of the two into one statement with the conclusions that make it usable — convergence, continuity, non-negativity and the normalisation — and the proof below, which reaches them without Bochner's theorem and without the differentiability hypothesis their Theorem (9.3.15) carries. This is the version used in practice.
Idea of the proofOriginal
Why a correlation function must be a Fourier transform. The variance of any linear combination Z=∑jzjX(tj) is non-negative, and expanding that variance writes it as the double sum ∑j,kc(tk−tj)zjzk. So non-negative definiteness of c is not an extra assumption: it is the statement that variances are non-negative, rewritten. Bochner's theorem says that this property, together with continuity and the normalisation ρ(0)=1, is exactly what characterises characteristic functions. The spectral theorem is therefore Bochner's theorem applied to ρ, and the only work is to check that ρ is continuous everywhere, not merely at the origin.
Why continuity at one point gives continuity everywhere. The increment c(t+h)−c(t) is an expectation of X(0) against the increment X(t+h)−X(t). The Cauchy-Schwarz inequality bounds it by the product of two root-mean-squares, and stationarity makes the second of them 2c(0)−c(h)−c(−h), a quantity that involves h alone. So the modulus of continuity of c at a general t is controlled by its behaviour at 0 — uniformly in t. A function that is non-negative definite and continuous at the origin is automatically uniformly continuous.
Why the discrete-time statement lives on an interval. Sampling at integer times cannot distinguish a frequency λ from λ+2π: the two oscillations agree at every integer. The spectral mass at all these frequencies is therefore indistinguishable and may be piled up at one representative in (−π,π]. That is aliasing, and it is the reason a sampled signal has a spectrum only up to the Nyquist frequency.
Why absolute summability is enough, and why it gives non-negativity for free. With ∑n∣ρ(n)∣<∞ the Fourier series for f converges uniformly, so it may be integrated term by term against einλ, and orthogonality returns ρ(n). Non-negativity of the sum is less obvious, but it comes from the same source as non-negative definiteness: the Cesàro means of the series are, up to a constant, the expected squared modulus of ∑n≤NXneinλ, which cannot be negative.
Proof(Hermitian symmetry and non-negative definiteness)AdaptedAdapted from Grimmett & Stirzaker (2020), §9.3, proof of Theorem (9.3.2), p. 439
Part 1. By stationarity we may compute c(−t) at the time s=t:
c(−t)=cov(X(t),X(0))=cov(X(0),X(t))=c(t),
the middle step being the antisymmetry cov(C2,C1)=cov(C1,C2) of the complex covariance.
Part 2. Fix t1,…,tn and z1,…,zn and set Z=∑jzjX(tj). Covariance is linear in its first slot and conjugate-linear in its second, so
where the last step uses stationarity in the form cov(X(tj),X(tk))=c(tk−tj). A variance is non-negative, so the double sum is non-negative. □
Since ρ=c/c(0) and c(0)=var(X(0))>0, the function ρ inherits both properties and satisfies ρ(0)=1.
Proof(The spectral theorem)AdaptedAdapted from Grimmett & Stirzaker (2020), §9.3, proof of Theorem (9.3.4), p. 439, and §5.7, Theorem (5.7.3) and the statement of Bochner's theorem, p. 204
Assume, as we may, that EX(t)=0, so that c(t)=E(X(0)X(t)).
Step 1: an estimate for the increments of c. For real t and h,
c(t+h)−c(t)=E[X(0)(X(t+h)−X(t))],
by linearity of the expectation. Taking moduli inside and applying the Cauchy-Schwarz inequality,
Step 2: evaluate the second factor by stationarity. Expand the square and use E∣X(u)∣2=c(0), valid for every u, together with cov(X(t),X(t+h))=c(h) and cov(X(t+h),X(t))=c(−h):
E(∣X(t+h)−X(t)∣2)=2c(0)−c(h)−c(−h).
Only h survives on the right. Therefore
c(t+h)−c(t)≤c(0)(2c(0)−c(h)−c(−h))for all t.
Step 3: continuity at 0 upgrades to uniform continuity. By hypothesis ρ, and therefore c=c(0)ρ, is continuous at 0, so c(h)→c(0) and c(−h)=c(h)→c(0) as h→0. The bound of step 2 therefore tends to 0 as h→0, uniformly in t. So c is uniformly continuous on R, and so is ρ=c/c(0).
Step 4: apply Bochner's theorem. By the previous proof, ρ is non-negative definite in the form ∑j,kρ(tk−tj)zjzk≥0; replacing each zj by zj and exchanging the names of the summation indices turns this into ∑j,kρ(tj−tk)zjzk≥0, which is the form in which Bochner's condition is stated. We also have ρ(0)=1 and ∣ρ(t)∣≤1, the latter by the Cauchy-Schwarz inequality, and uniform continuity from step 3. Bochner's theorem — quoted by Grimmett and Stirzaker without proof, as the converse half of their Theorem (5.7.3) — says that any function carrying those three properties arises as the characteristic function of a distribution function F. That is the assertion:
ρ(t)=∫−∞∞eitλdF(λ).□
Uniqueness of F is the uniqueness theorem for characteristic functions (see Characteristic functions).
Proof(The proposition in discrete time)Original
Assume EXn=0 and ∑n∣ρ(n)∣<∞.
Step 1: convergence and continuity. Each term of ∑nρ(n)e−inλ has modulus ∣ρ(n)∣, so by the Weierstrass test the series converges absolutely and uniformly in λ. A uniform limit of continuous functions is continuous, so f is continuous on (−π,π], and 2πf is bounded by ∑n∣ρ(n)∣.
Step 2: the coefficients come back. Uniform convergence permits term-by-term integration, so for each m∈Z
∫−ππeimλf(λ)dλ=2π1n∑ρ(n)∫−ππei(m−n)λdλ=ρ(m),
because ∫−ππeikλdλ equals 2π when k=0 and 0 otherwise: this is the orthogonality of the exponentials over one period that underlies Fourier series, and the display says that ρ(−m), read the other way round, is the mth Fourier coefficient of 2πf. Taking m=0 gives ∫−ππf=ρ(0)=1.
Step 3: non-negativity, through a Cesàro mean. For N≥1 and λ fixed, put ZN=∑n=1NXneinλ and consider
using E(XnXm)=cov(Xn,Xm)=c(m−n). The summand depends on n and m only through k=m−n, and for each k with ∣k∣<N there are exactly N−∣k∣ pairs, so
SN(λ)=2π1∣k∣<N∑(1−N∣k∣)ρ(k)e−ikλ.
Step 4: let N→∞. Each term converges to ρ(k)e−ikλ/2π, and the terms are dominated in modulus by ∣ρ(k)∣/2π, which is summable; so by dominated convergence for series, SN(λ)→f(λ). A limit of non-negative numbers is non-negative, so f(λ)≥0 for every λ.
Step 5: identify f as the spectral density. By steps 1-4, f is a continuous probability density on (−π,π], and by step 2 the distribution function F(λ)=∫−πλf satisfies ρ(n)=∫(−π,π]einλdF(λ) for all n. This is the discrete-time spectral representation, and F is unique, so f is the spectral density. □
Remark(Where the hypotheses are used, conventions, and the other Wiener-Kinchin theorem)AdaptedAdapted from Grimmett & Stirzaker (2020), §9.3, Theorem (9.3.4) and the discussion of (9.3.11)-(9.3.14) and (9.3.18), pp. 439-442, and Examples (9.3.19) and (9.3.20), p. 442; Riley, Hobson & Bence (2006), §13.1.8, eqs. (13.40)-(13.43) and the energy spectrum, pp. 449-450
Strict positivity of the variance is needed only to divide by c(0). A process that is almost surely constant has c≡0 and no autocorrelation function; there is nothing to normalise and no spectral distribution.
Continuity of ρ at the origin is used exactly once, in step 3 of the proof, and it is not a technicality. Two different objects go by the name "white noise" in continuous time, and only the first of them is a process of the kind defined above; it is worth keeping them apart, because the hypothesis fails for one and the definition fails for the other.
A process with orthogonal values. Take c(0)=σ2>0 and c(t)=0 for every t=0. This c is non-negative definite — for distinct times the double sum of the statement above collapses to σ2∑j∣zj∣2 — so by the characterisation quoted there it really is the autocovariance of a process, and that process satisfies the definition of weak stationarity given here: the mean is constant and the covariance depends on the lag alone. What fails is the hypothesis, and with it the conclusion. Its ρ is 1 at the origin and 0 elsewhere, no characteristic function jumps like that, and so this process has no spectral distribution at all — not an unusual one, none. Nor does a formal transform rescue it: ρ vanishes off a single point, so ∫ρ(t)e−iλtdt=0, which is not a density. The lesson is that weak stationarity by itself buys nothing; the spectral theory begins at the hypothesis of continuity.
The white noise of physics. The delta-correlated object c(t)=σ2δ(t) is a different thing and is not covered by the definitions above at all, because c(0) is infinite: it has no variance, so it is not a second-order process and has no autocorrelation function to be continuous or discontinuous. It is this object, not the one in the previous paragraph, that carries the flat spectrum S(ω)=σ2 of infinite total mass; the infinite mass is precisely the infinite variance, seen in the frequency domain. White noise in this sense has to be handled as a generalised process, as the formal derivative of Brownian motion — which is what Brownian motion and The Ito integral and Ito's formula do instead. In the laboratory neither idealisation is met: a real noise source is flat only up to some bandwidth, beyond which S rolls off and the total variance is finite.
Only weak stationarity is used anywhere above. The spectral theory never looks beyond the mean and the covariance. This is the practical reason for preferring the weaker hypothesis: to fit a spectrum one needs to assume only that the first two moments are time-invariant, not that the whole law is.
In discrete time the spectral distribution is unique only on one period. The representation on R is not unique, since einλ is 2π-periodic in λ; uniqueness is restored by folding all the mass into (−π,π]. A sampled process therefore cannot tell a frequency from its aliases, and a spectral density estimated from samples at spacing Δ is defined only for ∣ω∣≤π/Δ.
Absolute summability is sufficient but not necessary for a spectral density to exist. Grimmett and Stirzaker say so explicitly at their (9.3.18). Their Example (9.3.20) shows the opposite extreme: for Xn≡Y with varY=1, one has ρ(n)=1 for all n, which is not summable; the spectral distribution is the unit mass at λ=0 and there is no density. Their Example (9.3.19), discrete white noise with ρ(n)=δn0, is the other extreme: the spectral density is the uniform density on [−π,π], and the spectrum is all of [−π,π].
The hypothesis of the inversion theorem, as printed. Grimmett and Stirzaker's Theorem (9.3.15) assumes that F is differentiable with derivative f, and then asserts the inversion formula "at every point λ at which f is differentiable". We state it exactly as they do. The proposition above avoids the question by assuming absolute summability outright, which yields a continuous f and the formula at every λ.
Conventions. Three normalisations are in common use, and they differ by factors of 2π and by whether the variance is divided out.
Object
Definition here
Relation
Spectral density f (as above)
ρ(t)=∫eitλf(λ)dλ
∫f=1
Power spectral density S (physics)
S(ω)=∫−∞∞c(t)e−iωtdt
S=2πc(0)f, and 2π1∫S(ω)dω=c(0)
One-sided spectrum (engineering)
S+(ω)=2S(ω) for ω>0, S+(0)=S(0)
defined for real processes only
In the handbook's Fourier convention, g^(k)=∫g(t)e−iktdt with the 2π on the inverse transform (see The Fourier transform), the first two lines read f=ρ^/2π and S=c^. Always check which of these a quoted "spectral density" means before comparing numbers.
The other theorem of the same name. (On the spelling, see the third warning in the opening note: Riley, Hobson and Bence write Kinchin where Grimmett and Stirzaker write Khinchin, and we follow whichever source is being reported.) Riley, Hobson and Bence prove, and call the Wiener-Kinchin theorem, a purely deterministic statement: for square-integrable f and g with cross-correlation C(z)=∫f(x)g(x+z)dx, one has C^(k)=2πf^(k)g^(k) in their symmetric convention; specialising to g=f makes the transform of the autocorrelation proportional to ∣f^(k)∣2, which they call the energy spectrum of f. No probability appears, and no stationarity: the two results are different theorems sharing a name and a shape. The probabilistic one above concerns an ensemble average that exists because the process is stationary; the deterministic one concerns a single square-integrable function, which a stationary process's sample path is not. The bridge between them, which is what makes the shared name defensible, is that a time average of ∣∫−TTX(t)e−iωtdt∣2/2T estimates S(ω) — the content of the periodogram, treated in The periodogram.
Method(Finding the spectral density of a stationary process)Original
Check weak stationarity. Confirm that the mean is constant and that cov(X(s),X(s+t)) does not depend on s. If the process is defined by a recursion or a differential equation driven by noise, this is usually a condition on the coefficients: for an autoregression Xn=∑j=1rαjXn−j+Zn the condition is that the roots of zr−α1zr−1−⋯−αr all lie inside the unit circle.
Compute the autocovariance c(t), then ρ(t)=c(t)/c(0). For a process built linearly from orthogonal noise, expand the product E(XnXn+m) and discard every term in which two different noise variables meet.
Check continuity at t=0 in continuous time; in discrete time there is nothing to check.
Transform. In discrete time, if ∑∣ρ(n)∣<∞, use f(λ)=(2π)−1∑nρ(n)e−inλ on (−π,π]; for a real process this is (2π)−1(1+2∑n≥1ρ(n)cosnλ). In continuous time use f(λ)=(2π)−1∫ρ(t)e−itλdt when ρ is integrable, which is Fourier inversion applied to the representation ρ(t)=∫eitλf(λ)dλ; otherwise recognise ρ as a known characteristic function and read off its distribution, which is often quicker than inverting an integral.
Check the answer three ways.f must be non-negative; ∫f must be 1; and f must be symmetric under λ↦−λ if the process is real.
Convert to the convention you need using the table in the remark, and record which one you used.
Read off what you wanted. The fraction of the variance carried by a frequency band B is ∫Bf; a spectral density with a sharp peak at λ0 means a quasi-periodic component of period 2π/λ0; a spectral density that is flat near λ=0 and rolls off at λc means a correlation time of order 1/λc.
Example(A moving average of order two)Original
Let {Zn:n∈Z} be orthogonal random variables, each of mean 0 and variance 1. Define a real-valued process by
Xn=Zn+21Zn−1−41Zn−2,n∈Z.
Write a0=1, a1=21, a2=−41 and ak=0 otherwise.
Step 1: stationarity and the autocovariance. Each Xn has mean 0. Since E(ZjZk)=δjk,
and c(m)=0 for ∣m∣≥3, because then no two of the shifted index sets overlap. Hence
ρ(1)=21/163/8=72,ρ(2)=21/16−1/4=−214.
Step 2: the spectral density. Only five terms of the series are non-zero, so absolute summability is immediate and the proposition applies. For a real process ρ(−n)=ρ(n), so
Step 3: the three checks. Symmetry in λ is visible. The integral is 1, since ∫−ππcos(nλ)dλ=0 for n=1,2 and ∫−ππdλ=2π. For non-negativity, note the general identity for a moving average: with A(z)=∑kakzk,
so that f(λ)=∣A(e−iλ)∣2/(2πc(0))≥0. Here A(z)=1+21z−41z2. As a numerical check at λ=0: the formula gives f(0)=(1+74−218)/(2π)=42π25=0.1894702, while A(1)=45 gives ∣A(1)∣2/(2πc(0))=(25/16)/(2π⋅21/16)=25/(42π), the same. At the other end, f(π)=(1−74−218)/(2π)=1/(42π)=0.0075788, which is the minimum of f and is positive.
Step 4: a number to check. What fraction of the variance of X lies at frequencies ∣λ∣≤π/2, that is, in periods of four sampling intervals or longer? It is
Numerically this is 0.681891. The filter has tilted the flat spectrum of the input noise towards low frequencies: white noise would put exactly half the variance in this band, and the moving average has moved about 18 per cent of the total into it. That is the signature of the positive coefficient on Zn−1, which makes consecutive values of X positively correlated, ρ(1)=2/7>0.
Application(Noise spectra in the laboratory)Original
What a spectrum analyser measures. Point an instrument at a fluctuating voltage, a photocurrent, a seismometer trace, or a sea-surface elevation record, and what comes back is an estimate of the power spectral density S(ω)=∫c(t)e−iωtdt. The theorem above is what guarantees that this is a well-posed thing to measure: provided the record is stationary and the correlation function is continuous at the origin, a non-negative S exists and determines, and is determined by, the correlation function. The estimator itself is the periodogram (The periodogram), and the standard worry about it — that it is noisy however long the record, and must be smoothed — is a statistical matter, not a failure of the theorem.
Lorentzian line shapes. A process whose correlations decay exponentially, c(t)=σ2e−∣t∣/τ, has ρ(t)=e−∣t∣/τ, the characteristic function of a Cauchy distribution scaled by 1/τ; its spectral density is therefore the Lorentzian
This single formula is behind the shape of the pressure-broadened spectral line, the roll-off of the Ornstein-Uhlenbeck velocity spectrum (see Stochastic differential equations), the noise spectrum of a single relaxing dipole, and the Lorentzian fitted to a dynamic light-scattering correlogram to extract a diffusion coefficient. The correlation time is read off as the inverse of the half-width: a spectrum flat up to ω∼1/τ and falling as ω−2 beyond it.
White noise and its colouring. A flat spectral density means an autocorrelation concentrated at zero lag — no memory. Johnson-Nyquist noise in a resistor, and shot noise in a current, are white over the bandwidths usually of interest, and the linear circuit or detector that follows multiplies the spectral density by ∣H(ω)∣2, the squared modulus of the transfer function. That multiplication rule is the continuous-time twin of the identity f=∣A(e−iλ)∣2/(2πc(0)) verified in the example above, and it is the content of Linear filters: filtering a stationary process reshapes its spectrum pointwise and does nothing else.
Transport coefficients from correlation functions. In statistical mechanics the zero-frequency value of a spectral density is a transport coefficient. For the velocity of a tagged particle in three dimensions, the Green-Kubo relation
D=31∫0∞⟨v(0)⋅v(t)⟩dt
expresses the diffusion coefficient as half the zero-frequency spectral density of the velocity, divided by three for the dimensions; the integral converges precisely because the velocity autocorrelation decays, which is the same condition that gives the velocity process a spectral density. Similar expressions give the shear viscosity from the stress autocorrelation and the thermal conductivity from the heat-flux autocorrelation. In every case the quantity actually computed in a molecular-dynamics run is the autocovariance function of a stationary process, and the transport coefficient is read off its transform.
A caution about stationarity. Climate records, financial series and detector baselines drift, and a drifting series is not stationary; its sample "spectrum" will show a spurious rise at low frequency, because a trend looks like an oscillation of period longer than the record. The usual remedies — removing a fitted trend, differencing, or working with increments — are all attempts to recover a process to which this theorem applies.
Defines the Wiener process, shows that it is the scaling limit of a random walk in the sense that the rescaled walk's finite-dimensional distributions converge to it, sketches the proof that a version with continuous paths exists, and establishes the path properties - covariance σ2min{s,t}, quadratic variation σ2t, and total variation infinite on every interval.
Note(What has to be supplied, and what cannot be)Original
A random walk is specified by one number per step. A process in continuous time has uncountably many times, so a list of distributions is no longer enough to pin it down: two processes can have exactly the same distribution at every finite collection of times and yet have completely different paths. One may be continuous and the other wild.
This is the one structural difficulty in defining Brownian motion, and it shapes everything below. The distributional half of the definition is the easy half, and it is strongly suggested - though, as the idea block explains, not forced - by subdividing a time interval and looking at the central limit theorem: the increments ought to be normal, with variance proportional to elapsed time. That is where both of our sources stop arguing and start assuming. The path half - continuity - has to be imposed as a separate requirement as well, and then one must prove that the two halves are compatible. They are, and that is the existence theorem.
Two warnings about names. Grimmett and Stirzaker use Wiener process for the mathematical object and keep Brownian motion for the physical phenomenon it models; Folland uses the two names interchangeably for the process. We follow the first convention, and say "Wiener process" whenever precision matters.
Definition(Gaussian processes, and the Wiener process)AdaptedAdapted from Grimmett & Stirzaker (2020), §13.3, Definition (13.3.1) with the standardisation paragraph, pp. 580-581, and §8.5, pp. 416-417; Folland (1999), §10.5, postulates (10.20)-(10.22), pp. 330-331
A random processX={X(t):t≥0} is a family of random variables on one probability space (Ω,F,P) indexed by time. For a fixed ω∈Ω the function t↦X(t;ω) is a sample path. The finite-dimensional distributions of X are the joint distributions of the vectors (X(t1),…,X(tk)) over all finite sets of times.
X is a Gaussian process if every one of its finite-dimensional distributions is multivariate normal (see The multivariate normal distribution). The process X is said to have independent increments when, for any finite list of pairwise disjoint time intervals, the increments of X across those intervals form an independent family; it is enough to require this for two intervals [u,v] and [s,t] with u≤v≤s≤t and then for finitely many at a time.
Definition. Fix a constant σ2>0. A real-valued Gaussian process W={W(t):t≥0} with W(0)=w is called a Wiener process with variance parameter σ2 when it has all three of these properties.
(Independent increments.) Increments of W over disjoint time intervals are independent.
(Stationary normal increments.) Each increment W(s+t)−W(s), with s,t≥0, is distributed as N(0,σ2t).
(Continuity.) Every sample path of W is a continuous function of time.
W is standard if w=0 and σ2=1. Any Wiener process is reduced to a standard one by W1(t)=(W(t)−W(0))/σ. Because the law of W(s+t)−W(s) depends on t and not on s, the increments are called stationary.
Condition 2 forces E(W(s+t)−W(s))=0: a Wiener process has no drift. Condition 3 is a statement about the functions t↦W(t;ω), not about any distribution, and is the only part of the definition that the finite-dimensional distributions do not see.
Folland reaches the same three conditions by a different route, reading each as a physical postulate about a particle bombarded by molecules: it starts at a known place (his (10.20)); a collision has no lasting memory, so increments over disjoint time intervals are independent (his (10.21)); the mechanism is unchanging in time, so the law of W(t)−W(s) depends only on t−s; and the particle is pushed left and right equally, so the mean increment is zero. Normality and the proportionality of the variance to t−s are not deduced from these. Folland subdivides [s,t] into n equal pieces, says that the central limit theorem "suggests" that the increment should be normal and that additivity of the variance "strongly indicates" that it should be proportional to t−s, and then writes: "we are led to the third assumption", namely his (10.22), that W(t)−W(s) is normal with mean 0 and variance C(t−s). So normality and the variance law are a postulate for him too, as they are condition 2 for us; the subdivision argument recalled in the idea below is its motivation, not its proof.
Definition(Covariance, finite-dimensional distributions, and existence)AdaptedAdapted from Grimmett & Stirzaker (2020), §8.5, Lemma (8.5.1), p. 417, and §13.3, equation (13.3.2), transition density (13.3.4) and Theorem (13.3.19), pp. 581-587; Folland (1999), §10.5, equation (10.23) and Theorem 10.28, pp. 332-334
Let W be a Wiener process with variance parameter σ2 and W(0)=0.
Covariance. For all s,t≥0,
cov(W(s),W(t))=σ2min{s,t}.
Finite-dimensional distributions. For 0≤t1<t2<⋯<tk, the vector (W(t1),…,W(tk)) is multivariate normal with zero mean vector and covariance matrix V with entries vij=σ2min{ti,tj}. Conditions 1 and 2 of the definition therefore determine every finite-dimensional distribution, and so determine the conditional (transition) density of W(t) given W(s)=x, which for the standard process is
f(t,y∣s,x)=2π(t−s)1exp(−2(t−s)(y−x)2),s<t.
Existence. There is a probability space carrying a process that satisfies all three conditions of the definition. Equivalently, the probability measure determined by the finite-dimensional distributions in part 2 can be realised on the space of continuous paths.
Part 3 is the content of the definition that is not free. A Gaussian process with the finite-dimensional distributions of part 2 exists by general constructions; what has to be proved is that one of them may be chosen with every sample path continuous.
Definition(The Wiener process as the scaling limit of a random walk)Original
Let X1,X2,… be independent and identically distributed with E(X1)=0 and var(X1)=σ2, where 0<σ2<∞. Put S0=0 and Sm=X1+⋯+Xm, and for each n≥1 define the rescaled walk
Wn(t)=nS⌊nt⌋,t≥0,
in which space is compressed by n while time is compressed by n.
Claim. For every finite set of times 0≤t1<t2<⋯<tk,
(Wn(t1),…,Wn(tk))D(W(t1),…,W(tk))as n→∞,
where W is a Wiener process with variance parameter σ2 and W(0)=0.
That is: the finite-dimensional distributions of the rescaled walk converge to those of the Wiener process, whatever the step distribution, provided only that its mean is zero and its variance is finite and non-zero. No other scaling works: see the remark below.
Definition(Path properties)AdaptedAdapted from Grimmett & Stirzaker (2020), §13.2, p. 579 (nowhere differentiability), Exercise (13.2.1), p. 580, §8.5, Exercise (8.5.4), p. 418, and §8.8, Definition (8.8.5), p. 427
Let W be a Wiener process with variance parameter σ2, and fix t>0. Write
Qn=j=0∑n−1(W((j+1)t/n)−W(jt/n))2
for the sum of squared increments over the uniform partition of [0,t] into n pieces.
Quadratic variation.Qn→σ2t in mean square as n→∞, and almost surely along the subsequence n=2k. (The mean-square statement is the one Grimmett and Stirzaker set as Exercise (8.5.4); the almost sure statement along the dyadic subsequence is an addition of ours, and both are proved below.)
Unbounded total variation. Almost surely, the sample path of W has infinite total variation on every interval [a,b] with a<b; that is,
supr∑W(sr+1)−W(sr)=∞,
the supremum being over all finite increasing sequences a=s0<s1<⋯<sn=b.
Nowhere differentiability. Almost every sample path is a continuous function of t that is differentiable at no point. Grimmett and Stirzaker record this without proof, and we follow them. Part 2, which we do prove, gives a weaker statement of the same kind by itself: a function that is Lipschitz on an interval has finite total variation there, so almost surely no sample path is Lipschitz on any interval, and in particular none has a derivative that exists and is bounded throughout an interval.
Part 2 is why the Ito integral of The Ito integral and Ito's formula cannot be a Stieltjes integral taken path by path, and part 1 is what replaces it: the correct "dW2" is not negligible but equals σ2dt.
Idea of the proofOriginal
Three separate ideas are at work, one for each statement.
Why the increments are taken to be normal. Split [s,t] into n equal pieces. By stationarity and independence, W(t)−W(s) is a sum of n independent identically distributed increments, for everyn; in a word, its law is infinitely divisible. Two things should be separated here.
What the subdivision really settles is the variance. Write v(h) for the variance of an increment across an interval of length h, assumed finite. Independence and stationarity give v(s+t)=v(s)+v(t) for all s,t≥0, and v is non-decreasing because a longer interval contains a shorter one and the extra variance is non-negative; a non-decreasing solution of that additive equation is v(h)=v(1)h. So once the increments are assumed to have finite variance, proportionality to elapsed time is a genuine consequence of conditions 1 and stationarity, with σ2=v(1).
Normality is not. Being an n-fold convolution of identical pieces for every n, even with finite variance, does not make a law normal: the Poisson and gamma laws, and compound Poisson laws generally, are infinitely divisible and have finite variance. The missing ingredient is a negligibility (Lindeberg) condition saying that no single piece of the subdivision carries a non-negligible share of the sum; a Poisson increment fails it, because a single jump is a macroscopic part of the total. With such a condition the central limit theorem for triangular arrays applies and the limit is normal, and that is the honest content of the heuristic. Folland supplies no such condition and claims none: he says the central limit theorem "suggests" normality and then assumes it. Condition 2 is therefore a postulate - very well motivated, and the physically right one for a particle hit by enormously many enormously small kicks, which is exactly a negligibility hypothesis - but a postulate.
Why the rescaled walk converges. Over disjoint time intervals the rescaled walk accumulates disjoint blocks of steps, so its increments are exactly independent at every n, not just in the limit. Each single increment is a normalised sum of about n(tj−tj−1) steps, so the central limit theorem sends it to N(0,σ2(tj−tj−1)). A joint characteristic function that factorises into pieces, each of which converges, converges to the product of the limits; that product is the characteristic function of the Wiener vector.
Why a continuous version exists. Build the path by successive refinement. Fix its values at the times m/2n for n=1, then n=2, and so on, joining them by straight lines; passing from level n−1 to level n only moves the midpoints of the existing segments. The midpoint displacement at level j is normal with variance 2−j−1, and there are 2j−1 of them, so the largest of them is of order j2−j - small enough that the total displacement over all levels converges. Summable corrections produce a uniformly convergent series of continuous functions, so the limit is continuous; the first Borel-Cantelli lemma turns "summable" into "almost surely".
Why the quadratic variation survives. Each squared increment has mean σ2t/n and variance of the same order squared, so the sum has mean exactly σ2t and variance of order 1/n: adding n independent small terms averages the fluctuation away but not the mean. Since the increments themselves are of order n−1/2 and there are n of them, the sum of their absolute values is of order n, which is why the total variation diverges.
Proof(Covariance, and the scaling limit)Original
Step 1: the covariance. Take 0≤s≤t and W(0)=0. Write W(t)=W(s)+(W(t)−W(s)) and multiply by W(s)=W(s)−W(0):
E(W(s)W(t))=E(W(s)2)+E((W(s)−W(0))(W(t)−W(s))).
The two factors in the second term are increments over the disjoint intervals [0,s] and [s,t], hence independent by condition 1, and each has mean zero by condition 2; so that term vanishes. The first term is var(W(s))=σ2s by condition 2. Since all means are zero, cov(W(s),W(t))=σ2s=σ2min{s,t}, and the same holds with s and t exchanged.
The finite-dimensional distributions follow: (W(t1),…,W(tk)) is multivariate normal because W is a Gaussian process, its mean vector is zero, and its covariance matrix has just been computed.
Step 2: increments of the rescaled walk. Fix 0≤t1<⋯<tk and put t0=0, so Wn(t0)=0. Define the increments
For distinct j these sums use disjoint blocks of the Xi, so Dn,1,…,Dn,k are independent for every n - this is where the independence of the steps is used, and it is exact, not asymptotic.
Step 3: the characteristic function of each increment converges. Write dj=tj−tj−1. We show that for every c∈R
E(eicDn,j)⟶exp(−21c2σ2dj).
The statement allows t1=0, and then d1=0 and the sum defining Dn,1 is empty, so Dn,1=0 for every n and both sides above equal 1; that case is settled, and we assume dj>0 from here. Put mn=⌊ntj⌋−⌊ntj−1⌋, so that mn/n→dj>0 and hence mn→∞. Then
Here Yn is the sum of mn independent copies of X1 divided by mn, which is σ times the sum standardised in the usual way (divided by σmn). Since X1 is centred and has finite non-zero variance, the central limit theorem sends that standardised sum to N(0,1), so YnDN(0,σ2) and, by the direct half of the continuity theorem for characteristic functions (Characteristic functions), E(eiθYn)→e−θ2σ2/2 for each fixed θ.
The factor an is deterministic and tends to dj, but it moves with n, so the convergence just quoted cannot be applied at the argument can directly. Two elementary bounds remove the difficulty. First, Yn has mean 0 and variance exactly σ2 for every n, so E∣Yn∣≤σ by the Cauchy-Schwarz inequality. Second, ∣eiu−1∣≤∣u∣ for real u. Hence
while E(eicdjYn)→e−c2djσ2/2 by the previous paragraph with the fixed argument θ=cdj. Adding the two statements gives the display at the head of this step, and with it Dn,jDN(0,σ2dj).
Step 4: assembling the joint limit. Let θ1,…,θk be real and set cj=θj+θj+1+⋯+θk. Writing each Wn(tj) as the sum of the increments up to j and collecting terms,
j=1∑kθjWn(tj)=j=1∑kcjDn,j,
so by the independence of step 2 and the limit of step 3,
The identical algebra applied to a Wiener process W with parameter σ2, whose increments over [tj−1,tj] are independent and N(0,σ2dj), shows that the right-hand side is exactly the characteristic function of (W(t1),…,W(tk)) evaluated at (θ1,…,θk). Pointwise convergence of characteristic functions to a characteristic function gives convergence in distribution, by the continuity theorem. □
Proof(Existence of a version with continuous paths)AdaptedAdapted from Grimmett & Stirzaker (2020), §13.3, Theorem (13.3.19) and its proof, pp. 587-589
We construct the standard process on [0,1]; the general case follows by scaling (W(t)=σW1(t)+w) and by joining copies on [n,n+1] end to end.
Start from a Gaussian process X on [0,1] with X(0)=0, independent increments and X(s+t)−X(s)∼N(0,t). Such a process exists; the issue is its paths. Let Dn={m2−n:0≤m≤2n} be the level-n dyadic points and D=⋃nDn.
Step 1: the approximating polygons. Let Xn agree with X on Dn and be linear in between, so Xn is a continuous piecewise-linear function with 2n segments. Put X0≡0 and Zj=Xj−Xj−1, so that
Xn=j=1∑nZj.
Since Xj and Xj−1 agree at every point of Dj−1, the function Zj vanishes there; it is piecewise linear with corners only at points of Dj. A piecewise-linear function attains its extreme values at its corners, so
(One level is exceptional, and we follow Grimmett and Stirzaker in passing over it. With the convention X0≡0 the function Z1=X1 does not vanish on D0={0,1}, since Z1(1)=X(1); so at j=1 the supremum is a maximum over {21,1} rather than over the single odd point, and the two values have variances 21 and 1 rather than the 2−j−1=41 that step 2 will give. The slip is harmless, because a single level contributes the one almost surely finite quantity supt∣Z1(t)∣ to the sum in step 4, and the estimates of steps 2 and 3 are needed only for j≥2, where the description above is exact. Taking X0 to be the linear interpolant of X on D0 instead of 0 would remove the exception altogether.)
Step 2: the size of one correction. For odd m, Xj−1 interpolates linearly between the neighbouring even dyadics, so
where A=X(m2−j)−X((m−1)2−j) and B=X((m+1)2−j)−X(m2−j) are increments over adjacent intervals of length 2−j. By conditions 1 and 2 they are independent and each is N(0,2−j), so Zj(m2−j) is normal with mean 0 and variance 41(2−j+2−j)=2−j−1.
Step 3: a tail bound, level by level. For Y∼N(0,v) and x>0 the standard tail estimate gives P(∣Y∣>x)≤(2v/(x2π))e−x2/(2v). With v=2−j−1 this is at most x−12−j/2e−x22j. Combining with step 1 and the union bound over the 2j−1 odd points,
P(0≤t≤1sup∣Zj(t)∣>x)≤x2j/22j−1e−x22j.
Choose the threshold xj=cj2−jlog2 with c>1 fixed. Then xj22j=c2jlog2, so e−xj22j=2−c2j and xj2j/2=cjlog2, whence
P(0≤t≤1sup∣Zj(t)∣>xj)≤cjlog22j(1−c2)−1.
Because c>1 the exponent j(1−c2) is negative and linear in j, so these probabilities form a convergent geometric-type series.
Step 4: Borel-Cantelli, and uniform convergence. By the first Borel-Cantelli lemma (Convergence of random variables), almost surely only finitely many j have supt∣Zj(t)∣>xj. Since ∑jxj=clog2∑jj2−j/2<∞, it follows that
j≥1∑0≤t≤1sup∣Zj(t)∣<∞almost surely.
On that event the series ∑jZj converges uniformly on [0,1] by the Weierstrass criterion, so Xn converges uniformly to a limit Y, and Y is continuous as a uniform limit of continuous functions.
Step 5: Y has the right law. If t∈D then Xn(t)=X(t) for all large n, so Y(t)=X(t). For general t, take dyadics s→t; then E((X(s)−X(t))2)=∣s−t∣→0, so X(s)→X(t) in mean square, while X(s)=Y(s)→Y(t) almost surely by continuity of Y. Mean-square and almost-sure limits agree almost surely, so Y(t)=X(t) almost surely for every fixed t. Hence Y has the same finite-dimensional distributions as X, so conditions 1 and 2 hold for Y.
Step 6: removing the exceptional set. Let Ω′ be the set of ω for which t↦Y(t;ω) is continuous; steps 4 and 5 give P(Ω′)=1. Replace Ω by Ω′ with the induced σ-field and probability measure. No probability changes, and now every path is continuous, so condition 3 holds as stated. □
Proof(Quadratic variation, and the divergence of the total variation)Original
Fix t>0 and write Δj=W((j+1)t/n)−W(jt/n) for j=0,…,n−1, so Qn=∑jΔj2.
Step 1: mean and variance of Qn. The Δj are increments over disjoint intervals of equal length t/n, hence independent and each N(0,σ2t/n). Therefore E(Δj2)=σ2t/n and, using E(Y4)=3v2 for Y∼N(0,v),
var(Δj2)=3(nσ2t)2−(nσ2t)2=n22σ4t2.
Summing over the n independent terms,
E(Qn)=σ2t,var(Qn)=n2σ4t2.
Step 2: mean-square convergence.E((Qn−σ2t)2)=var(Qn)=2σ4t2/n→0, which is convergence in mean square to the constant σ2t. Note that the limit is deterministic although every Qn is random: the fluctuation is averaged away, the mean is not.
and these are summable in k. By the first Borel-Cantelli lemma the event on the left occurs for only finitely many k, almost surely, for each fixed ϵ; taking ϵ=1/i and intersecting the countably many full-probability events gives Q2k→σ2t almost surely. (Summability is what is needed here; mean-square convergence alone would only give a subsequence.)
Step 4: infinite total variation on one fixed interval. Work on [0,t]; the argument on a general [a,b] is the same with the origin shifted, because s↦W(a+s)−W(a) is again a Wiener process with the same parameter, by conditions 1 and 2. Let V be the total variation of the path over [0,t] and bound the quadratic variation by factoring out the largest increment:
Qn=j∑Δj2≤(jmax∣Δj∣)j∑∣Δj∣≤(jmax∣Δj∣)V.
Suppose V<∞ on a path. The path is continuous on the compact interval [0,t] by condition 3, hence uniformly continuous there, so maxj∣Δj∣→0 as the mesh t/n tends to 0. The displayed bound then forces Qn→0. But by step 3, almost surely Q2k→σ2t>0. So V=∞ almost surely.
Step 5: all intervals at once. Step 4 gives, for each fixed interval, a full-probability event on which the variation over it is infinite; part 2 of the statement claims one full-probability event that works for every interval simultaneously, and the two are not the same thing, since there are uncountably many intervals. Let Ωa,b be the event that the variation over [a,b] is infinite, and put
Ω∗=⋂{Ωa,b:a<b,a,b∈Q≥0},
an intersection of countably many events of probability 1, so P(Ω∗)=1. Fix ω∈Ω∗ and let [a,b] be any interval with a<b. Choose rationals a′ and b′ with a≤a′<b′≤b. Total variation is monotone under inclusion of intervals, because every partition of [a′,b′] extends to a partition of [a,b] and adjoining points can only increase the sum of absolute increments; hence the variation of this path over [a,b] is at least its variation over [a′,b′], which is infinite. So on Ω∗ every interval carries infinite variation, which is part 2. □
Remark(Where the hypotheses are used, and what fails without them)AdaptedAdapted from Grimmett & Stirzaker (2020), §13.2, equation (13.2.1) with properties (a)-(d), pp. 578-579, §13.3, p. 581, and §8.7, p. 423
Independent increments is used three times and in three different ways: to kill the cross term in the covariance computation, to make the blocks Dn,j of the rescaled walk independent (which is what lets a product of characteristic functions appear at all), and to make the midpoint corrections Zj of the existence proof normal with a known variance. Drop it and nothing in this unit survives. Fractional Brownian motion with Hurst exponent H=21 is Gaussian and self-similar with stationary increments, and Grimmett and Stirzaker report that its sample paths "may be taken to be" almost surely continuous and nowhere differentiable; but the only value of the exponent at which its increments come out independent is the Wiener value H=21, and its covariance is not min{s,t}.
Finite non-zero variance of the step distribution is where the scaling limit is decided. The compression S⌊nt⌋/n is tuned to σ2<∞: if the steps had infinite variance the central limit theorem would not apply and a different exponent, with a non-Gaussian stable limit, would be needed. Grimmett and Stirzaker make the same point in physical terms for the lattice walk of their §13.2: with lattice spacing δ and inter-jump time τ the limit is non-degenerate only if 31δ2/τ is held at a constant σ2 as δ,τ↓0. If δ2/τ→0 the displacement collapses to the constant 0; if δ2/τ→∞ its variance diverges. The factor 31 there is the probability of stepping along the chosen axis in their three-dimensional lattice; in our one-dimensional statement the corresponding constant is σ2 itself.
Zero mean of the steps keeps the limit centred. With E(X1)=μ=0 the same calculation applied to (S⌊nt⌋−μ⌊nt⌋)/n gives a Wiener process, while S⌊nt⌋/n→μt: the walk is then dominated by its drift, and the Wiener process describes only the fluctuation around it.
Continuity of paths is an extra condition, not a consequence. The finite-dimensional distributions do not determine it. Take a Wiener process W and an independent random time U, uniform on [0,1], and define W(t)=W(t) for t=U and W(U)=W(U)+1. For each fixedt we have P(U=t)=0 and hence P(W(t)=W(t))=1, so every finite-dimensional distribution is unchanged; yet every path of W jumps at U, so no path of it is continuous. This is exactly why part 3 of the existence statement has content, and why the proof of it is the long one.
Normality, not merely independence, fixes the size of the fluctuation. Step 1 of the last proof uses E(Y4)=3v2, a property of the normal distribution; it is the normality supplied by condition 2, not merely independence, that gives var(Qn)=O(1/n).
What the scaling-limit statement does not say. It is a statement about finitely many times at once. Convergence of the whole random path - as an element of the space of continuous functions, which would let one pass to the limit in quantities such as sup0≤s≤1Wn(s) - is Donsker's invariance principle, and it is not treated in the sources used here; our claim should not be read as asserting it.
Method(Computing with a Wiener process)Original
Almost every elementary computation reduces to four moves.
Re-express in increments. Rewrite the quantity in terms of W(t1)−W(0),W(t2)−W(t1),… over disjoint intervals. These are independent and normal, so the joint law is a product of one-dimensional normals.
Use the covariance for anything quadratic. For means, variances and covariances of linear or quadratic functionals, do not integrate densities; use linearity together with cov(W(s),W(t))=σ2min{s,t}. For instance var(W(t)−W(s))=σ2t−2σ2min{s,t}+σ2s=σ2∣t−s∣.
Use linearity to stay Gaussian. Any finite linear combination ∑jθjW(tj) is a one-dimensional normal with mean 0 and variance ∑i,jθiθjσ2min{ti,tj}; a limit in distribution of such combinations (for example a limit of Riemann sums) is normal too, since the normal family is closed under limits in distribution. So probabilities about such quantities are read off a single standard normal table.
Use the scaling law to remove parameters. For standard W and c>0, the process t↦c−1/2W(ct) is again a standard Wiener process, because it is Gaussian with mean zero and covariance c−1min{cs,ct}=min{s,t} and has continuous paths. (For a centred Gaussian process that covariance already forces the other two conditions: increments over disjoint intervals [u,v] and [s,t] with u≤v≤s≤t have covariance min{v,t}−min{v,s}−min{u,t}+min{u,s}=v−v−u+u=0, and jointly normal variables with zero covariance are independent, while var(W(t)−W(s))=∣t−s∣ by move 2.) Use it to reduce a question on [0,T] to one on [0,1].
A fifth move is to use the transition density when a genuine distribution function is wanted: conditionally on W(s)=x, W(t) is N(x,σ2(t−s)) for t>s.
Example(The time-average of a Brownian path)Original
Let W be a standard Wiener process (σ2=1, W(0)=0) and let
A=∫01W(s)ds.
This integral makes sense path by path: every path is continuous by condition 3, hence Riemann integrable on [0,1]. We find the law of A, its correlation with the endpoint W(1), and a numerical probability.
Step 1: A is normal. Let An=n1∑r=1nW(r/n) be the right-endpoint Riemann sum. Each An is a finite linear combination of the W(r/n), hence normal with mean 0 by move 3 of the method. For each path, An→A as n→∞ by Riemann integrability, so An→A almost surely and therefore in distribution. The variances converge, since var(An)=n−2∑r,s≤nmin{r/n,s/n} is a Riemann sum for the double integral evaluated in step 2, and a limit in distribution of N(0,vn) variables with vn→v is N(0,v). So A is N(0,v) with v as computed next.
Step 2: the variance. Using E(W(s)W(u))=min{s,u} and exchanging expectation with the integrals,
where the factor 2 comes from the symmetry of min{s,u} in the two variables. So
A∼N(0,31),sd(A)=1/3=0.577350.
Step 3: correlation with the endpoint. Likewise
cov(A,W(1))=∫01min{s,1}ds=∫01sds=21,
so, since sd(W(1))=1,
corr(A,W(1))=1/3⋅11/2=23=0.866025.
The average over [0,1] and the final value are strongly but not perfectly correlated, as one expects: the average is dominated by the later part of the path, where the path has already moved.
Step 4: a probability. The chance that the path spends its first unit of time above the level 21 on average is
Step 5: a check on the quadratic variation. For the same process on [0,3] with n=1000 equal subintervals, the proof above gives E(Q1000)=σ2t=3 exactly and sd(Q1000)=2σ4t2/n=18/1000=0.134164. A simulated path should return a value of Q1000 within about 0.27 of 3 roughly 95% of the time - a sharp and easily checked prediction, and one that fails for any process with differentiable paths, for which Qn→0.
Application(Diffusion, the heat equation, and measuring molecules)Original
The transition density is the diffusion kernel. Let u(y,t) be the probability density of the position at time t of a particle whose displacement is a Wiener process with parameter σ2, started from a known point. By part 2 of the first statement, u is the Gaussian kernel with variance σ2t, which is exactly the fundamental solution of
∂t∂u=2σ2∂y2∂2u,
the diffusion equation with diffusion coefficient D=σ2/2 (see The heat kernel). Strauss states the correspondence in the other direction: the probability that a particle started at x lies in (a,b) at time t is the integral of the source function over that interval, so a concentration profile transported by Brownian motion obeys the diffusion equation. The general version of this dictionary, for processes with drift and with a position-dependent diffusion coefficient, is The Fokker-Planck equation; the particular case in which the velocity rather than the position diffuses gives the Ornstein-Uhlenbeck process of Stochastic differential equations.
Harmonic functions and the Dirichlet problem. Strauss gives the companion steady-state statement: for Brownian motion in a container D, stopped the moment it reaches the boundary, the probability u(x) that the particle stops on a chosen boundary piece C1 satisfies Δu=0 in D with u=1 on C1 and u=0 on the rest. This is a Dirichlet problem (Harmonic functions), and it is the basis of "walk on spheres" Monte Carlo solvers, which estimate the potential at a single point of a complicated domain by simulating paths rather than meshing the domain.
Einstein's relation and the size of a molecule. The physically measurable content of the model is E((W(t)−W(0))2)=σ2t=2Dt: the mean square displacement grows linearly in time, not quadratically as it would for ballistic motion. For a sphere of radius a in a fluid of viscosity η at absolute temperature T, the Stokes-Einstein relation gives D=kBT/(6πηa). For a bead of radius a=0.5μm in water at 20∘C (η=1.002×10−3Pas, T=293K, kB=1.380649×10−23JK−1),
so after 10 seconds the root-mean-square displacement along one axis is 2Dt=8.567×10−12=2.93μm. Running this backwards - measure the displacements, infer D, and read off kB and hence Avogadro's number - is Perrin's experiment, and it is why a theorem about path regularity has a place in a physics handbook.
Where the irregularity matters. The infinite total variation of part 2 of the path statement is not a curiosity. It means that integrals ∫H(s)dW(s) cannot be defined path by path as Stieltjes integrals, which is why a stochastic integral must be built as a limit in mean square instead, and why its chain rule carries the extra second-order term supplied by the quadratic variation σ2t. That is the subject of The Ito integral and Ito's formula. It also means that a numerical simulation on a time grid of size h reproduces increments of size h, not h: halving the step size reduces the step length only by a factor 2, which is also why the error analysis of a Monte Carlo estimate of such a path functional behaves differently from that of an ordinary quadrature (Monte Carlo methods).
Defines filtrations, martingales, sub- and supermartingales and stopping times, proves that a stopped martingale is a martingale and then the optional stopping theorem in two usable forms, shows by the doubling strategy exactly which hypothesis fails when the conclusion does, and uses the theorem to compute the mean waiting time for a pattern in coin tosses.
A martingale is the mathematics of a fair game: whatever has happened so far, the expected change at the next step is zero. That single property is weaker than independence and weaker than the Markov property, and it is nonetheless enough to prove convergence theorems, concentration inequalities and - the subject of this unit - a rule for exchanging a random stopping time with an expectation.
The use is almost always the same. One wants a quantity attached to a random time: the chance of hitting one barrier before another, the mean number of steps until a pattern appears, the probability that a running statistic ever crosses a threshold. Direct computation means summing over all the ways the history could have gone. Instead one finds a martingale, argues that its mean cannot have changed by the time one stops, and reads off the answer from a single equation. The catch is that the mean can change: a gambler who can borrow without limit and who stops as soon as he is ahead is certain to end ahead. The theorems below say precisely which hypothesis that gambler violates.
Definition(Filtrations, martingales, and stopping times)AdaptedAdapted from Grimmett & Stirzaker (2020), Ch. 12, Definitions (12.1.1), (12.1.8) and (12.1.9), pp. 533-536, and Definition (12.4.1) with the construction of $\mathcal{F}_T$, p. 550
Throughout, (Ω,F,P) is a probability space and time runs over n=0,1,2,….
Filtration, adapted. A filtration is a family F={Fn:n≥0} of sub-σ-fields of F that increases, Fn⊆Fn+1. Read Fn as everything an observer knows at time n. A sequence Y={Yn} is adapted to F if each Yn is Fn-measurable, that is, if Yn is known at time n. The natural filtration of a sequence X={Xn} is Fn=σ(X0,X1,…,Xn).
Martingale. The adapted pair (Y,F) is a martingale if, for every n≥0,
E∣Yn∣<∞,E(Yn+1∣Fn)=Yn.
Here E(⋅∣Fn) is conditional expectation given a σ-field (Conditional expectation); it is a random variable, and the identity, like every identity between conditional expectations in this unit, is asserted up to an event of probability zero. Taking expectations in the defining identity gives E(Yn)=E(Y0) for every n: the mean of a martingale never moves.
When the filtration is the natural one of some observed sequence X, one says that Y is a martingale with respect to X; the condition is then E(Yn+1∣X0,…,Xn)=Yn. A finite family {Yn:0≤n≤N} satisfying the same conditions is also called a martingale.
Submartingale, supermartingale. The adapted pair (Y,F) is a submartingale if E(Yn+)<∞ and E(Yn+1∣Fn)≥Yn for all n, and a supermartingale if E(Yn−)<∞ and E(Yn+1∣Fn)≤Yn for all n, where x+=max{0,x} and x−=−min{0,x}. Note that these moment conditions are one-sided, and so weaker than E∣Yn∣<∞. A sequence is a martingale exactly when it is both; Y is a submartingale exactly when −Y is a supermartingale. The names are the opposite way round from what the inequalities suggest; they come from the theory of harmonic functions, where a subharmonic function lies below the harmonic one with the same boundary values.
Stopping time. A random variable T with values in {0,1,2,…}∪{∞} is a stopping time for F if
{T=n}∈Fnfor every n≥0,
equivalently if {T≤n}∈Fn for every n, and hence also {T>n}∈Fn. The content is that the decision to stop at time n may use the first n observations and nothing later. A stopping time is allowed to be infinite: "never" is a permitted value.
The standard example is a first passage time: if X is adapted and B⊆R, then TB=min{n:Xn∈B} (with TB=∞ if X never enters B) is a stopping time, since {TB=n} is decided by X0,…,Xn. By contrast "the last time X visits B" is not a stopping time: deciding it requires the whole future.
The σ-field at a stopping time.FT is the collection of events A with A∩{T≤n}∈Fn for every n. It is a σ-field, and it is the information available at the random time T.
We write a∧b=min{a,b}, so YT∧n is the process watched until T and then held fixed.
Definition(Stopping a martingale)AdaptedAdapted from Grimmett & Stirzaker (2020), §12.4, Theorem (12.4.5) and the optional sampling theorem (12.4.11), pp. 551-552, and §12.5, Theorems (12.5.1) and (12.5.9) with Example (12.5.12), pp. 554-556
Throughout this statement T is a stopping time for the filtration F, and (Y,F) is a martingale.
A. The stopped process.Zn=YT∧n defines a martingale (Z,F); if (Y,F) is only a submartingale, the same holds with "submartingale" in place of "martingale". In particular
E(YT∧n)=E(Y0)for every finite n.
If moreover P(T≤N)=1 for some deterministic N, then T∧N=T and so E(YT)=E(Y0): for a bounded stopping time no further hypothesis is needed. This last statement is the optional sampling theorem, and it has a sharper form that is where the σ-field FT earns its place: if T1≤T2≤⋯ are stopping times, each bounded by a deterministic constant, then the sampled pairs (YTj,FTj) again form a martingale. Watching a fair game at a sequence of bounded, non-anticipating times gives a fair game, with FTj as the information available at the jth observation. The boundedness is essential: with Y the symmetric simple random walk and T1 its first passage to level 1, a single unbounded sampling time already gives E(YT1)=1=0=E(Y0).
B. Optional stopping theorem.E(YT)=E(Y0) provided all three of
P(T<∞)=1,
E∣YT∣<∞,
E(Yn1{T>n})→0 as n→∞.
C. A sufficient condition that is easier to check.E(YT)=E(Y0) provided
P(T<∞)=1 and E(T)<∞, and
there is a constant c with E(∣Yn+1−Yn∣∣Fn)≤c for all n<T.
D. Wald's equation. Let X1,X2,… be independent and identically distributed with E∣X1∣<∞ and mean μ, let Sn=X1+⋯+Xn, and let T be a stopping time with E(T)<∞ for the filtration Fn=σ(X1,…,Xn). Then
E(ST)=μE(T).
Grimmett and Stirzaker state this with the filtration generated by the martingale Yn=Sn−nμ rather than by the data, that is with σ(Y1,…,Yn) in place of σ(X1,…,Xn). The two are the same σ-field: Yn is a function of X1,…,Xn, and conversely Sn=Yn+nμ and Xn=Sn−Sn−1 recover the data from the Y's. So nothing turns on the choice, and we use the natural filtration of the data throughout.
Version C is the one used in practice: both of its hypotheses are statements about the step sizes and about T alone, and neither mentions YT.
Idea of the proofOriginal
Stopping does not break fairness. Between times n and n+1 the stopped process either moves exactly as Y does, on the event {T>n}, or does not move at all. Which of the two happens is decided by time n, because T is a stopping time; so, conditionally on Fn, the increment is either a fair increment or zero, and in both cases has conditional mean zero. Nothing is needed beyond the definition of a stopping time, which is exactly why this part is free.
Why that is not yet the theorem. Part A controls Y at the time T∧n, not at T. The two differ only on the event {T>n}, whose probability tends to zero; but "the event is unlikely" is not "the contribution is small", because Yn may be huge precisely there. The one-line identity
YT=YT∧n+(YT−Yn)1{T>n}
isolates the whole difficulty in the two terms E(YT1{T>n}) and E(Yn1{T>n}). The first dies automatically under hypothesis 2, being the tail of a convergent sum; the second is hypothesis 3, and the examples show it is a genuine extra demand.
Why bounded increments and a finite mean suffice. Hypothesis 3 is awkward to verify, so one looks for a single dominating random variable. The total movement of the stopped process never exceeds the sum of all the increments up to T, so ∣YT∧n∣≤∣Y0∣+∑i≤T∣Yi−Yi−1∣. Bounding each conditional increment by c and summing over i with the weights P(T≥i) turns the right-hand side into cE(T), which is finite. A dominated sequence may be passed to the limit, and that is version C.
Proof(A: the stopped process is a martingale)AdaptedAdapted from Grimmett & Stirzaker (2020), §12.4, proof of Theorem (12.4.5), p. 551
Let Zn=YT∧n.
Step 1: Zn is Fn-measurable. Split according to the value of T:
Zn=t=0∑n−1Yt1{T=t}+Yn1{T≥n}.
Each Yt with t≤n is Fn-measurable because Y is adapted and the filtration increases; each indicator is Fn-measurable because T is a stopping time, using {T≥n}={T≤n−1}c. Hence so is the sum.
Step 2: Zn satisfies the moment condition. Here the two cases must be kept apart, because the definitions impose different moment conditions. If Y is a martingale, the display shows ∣Zn∣≤∑t=0n∣Yt∣, a finite sum of integrable variables, so E∣Zn∣<∞. If Y is only a submartingale we may not assume E∣Yt∣<∞ - a submartingale is required to satisfy only the one-sided condition E(Yt+)<∞ - and the bound above is unavailable. But Zn equals one of Y0,…,Yn pointwise, so Zn+≤∑t=0nYt+, whence
E(Zn+)≤t=0∑nE(Yt+)<∞,
which is exactly the one-sided condition a submartingale must satisfy. This is the bound Grimmett and Stirzaker use, and it is the reason their theorem is stated for submartingales in the first place.
Step 3: the increment. Subtracting the displays for n+1 and n, every term cancels except on the event that T has not yet occurred:
Zn+1−Zn=(Yn+1−Yn)1{T>n}.
The indicator is Fn-measurable, so it may be taken outside the conditional expectation:
E(Zn+1−Zn∣Fn)=1{T>n}E(Yn+1−Yn∣Fn)=0,
the last step by the martingale property of Y. (If Y is only a submartingale the final bracket is non-negative, and since the indicator is non-negative the product is too, which gives the submartingale version.) Taking expectations in E(Zn+1∣Fn)=Zn gives E(YT∧n)=E(Z0)=E(Y0).
Step 4: bounded stopping times. If P(T≤N)=1 then YT∧N=YT almost surely, and step 3 at n=N gives E(YT)=E(Y0). □
Proof(B and C: the optional stopping theorem)AdaptedAdapted from Grimmett & Stirzaker (2020), §12.5, proofs of Theorems (12.5.1) and (12.5.9), pp. 554-556
Part B. Assume hypotheses 1-3. Because P(T<∞)=1, the variable YT is defined almost surely, and splitting according to whether T exceeds n,
YT=YT∧n+(YT−Yn)1{T>n}.
Every term is integrable - YT by hypothesis 2, Yn and YT∧n by step 2 of the previous proof - so expectations may be taken term by term, and part A replaces E(YT∧n) by E(Y0):
E(YT)=E(Y0)+E(YT1{T>n})−E(Yn1{T>n}).
This holds for every n, so it holds in the limit. The third term tends to 0: that is hypothesis 3. For the second, decompose over the value of T,
E(YT1{T>n})=k=n+1∑∞E(YT1{T=k}),
which is the nth tail of the series ∑k≥0E(YT1{T=k}). Hypothesis 2 makes that series absolutely convergent (its absolute terms sum to E∣YT∣), and hypothesis 1 makes it sum to E(YT), so its tails tend to 0. Letting n→∞ leaves E(YT)=E(Y0).
Part C. Now assume instead that E(T)<∞ and that the conditional increments are bounded by c before T. Write Ri=∣Yi−Yi−1∣ for i≥1 and
R=i=1∑TRi=i=1∑∞Ri1{T≥i}.
Step 1: R has finite mean. The event {T≥i}={T≤i−1}c lies in Fi−1, so the indicator comes outside the conditional expectation and hypothesis 2 applies on it:
E(Ri1{T≥i}∣Fi−1)=1{T≥i}E(Ri∣Fi−1)≤c1{T≥i}.
Taking expectations and summing over i, with monotone convergence (Monotone convergence and Fatou's lemma) to exchange the sum and the expectation of non-negative terms,
E(R)≤ci≥1∑P(T≥i)=cE(T)<∞.
Step 2: domination. For every n, the stopped process has moved from Y0 by at most the total of the increments it has used, so ∣YT∧n∣≤∣Y0∣+R, and the dominating variable ∣Y0∣+R is integrable and does not depend on n.
Step 3: pass to the limit. Since E(T)<∞ we have T<∞ almost surely, so YT∧n→YT almost surely. Dominated convergence (Dominated convergence) applies by step 2, so E(YT∧n)→E(YT), while part A says every term of that sequence equals E(Y0). Hence E(YT)=E(Y0).
Part D. Put Yn=Sn−nμ with Y0=0 and Fn=σ(X1,…,Xn). Then E(Yn+1−Yn∣Fn)=E(Xn+1)−μ=0 by independence, so (Y,F) is a martingale; and
E(∣Yn+1−Yn∣∣Fn)=E∣Xn+1−μ∣=E∣X1−μ∣≤2E∣X1∣<∞,
a constant, so hypothesis C2 holds with that constant. With E(T)<∞, part C gives E(YT)=0, which is E(ST)=μE(T). □
Remark(Where the hypotheses are used: the doubling strategy)AdaptedAdapted from Grimmett & Stirzaker (2020), §12.1, Example (12.1.3), pp. 533-534, and §12.4, the symmetric random walk stopped at level 1 as an instance of failed optional sampling, p. 552
One example breaks every version of the theorem at once, and it is worth tracking exactly which hypothesis it violates.
A gambler stakes 1 on an even-money bet. After each loss he doubles: having lost the first n plays, he stakes 2n on play n+1. He stops at T, the first play he wins. Let Yn be his accumulated gain after n plays, with Y0=0.
(Y,F)is a martingale: before stopping, the next play adds +2n or −2n with probability 21 each, so the conditional mean increment is zero; after stopping it is zero outright.
T is a stopping time, P(T=n)=2−n for n≥1, so P(T<∞)=1 and even E(T)=∑nn2−n=2<∞. Hypotheses B1 and C1 hold.
YT=2T−1−(1+2+⋯+2T−2)=1 with probability 1, so E∣YT∣=1<∞ and hypothesis B2 holds.
And yet E(YT)=1=0=E(Y0).
So hypothesis B3 is what fails, and it fails by exactly the right amount. On {T>n} the gambler has lost n times, so Yn=−(1+2+⋯+2n−1)=−(2n−1) there, while P(T>n)=2−n. Hence
E(Yn1{T>n})=−(2n−1)2−n=−(1−2−n)⟶−1,
not 0; substituting into the identity in the proof of part B gives E(YT)=0+0−(−1)=1, the right answer. A vanishingly unlikely event carrying an exponentially large debt contributes a finite amount in the limit, which is the whole phenomenon.
Hypothesis C2 fails too, and only C2. On {T>n} the next stake is 2n, so E(∣Yn+1−Yn∣∣Fn)=2n, which no constant c bounds. Since C1 holds, this isolates the bounded-increment condition as the one doing the work in version C. The practical reading is familiar: the strategy is a certainty only for a gambler with unbounded credit, and Grimmett and Stirzaker note that his mean accumulated loss just before the winning play is infinite.
A second counterexample, where C2 holds and C1 fails. Let Y be a symmetric simple random walk started at 0 and let T=min{n:Yn=1} (Random walks). The walk is recurrent, so T<∞ almost surely and B1 holds; YT=1, so B2 holds; and the increments are bounded by 1, so C2 holds with c=1. Yet E(YT)=1=0=E(Y0). What fails is E(T)<∞, that is C1, and with it B3: since YT∧n has mean 0 and equals 1 on {T≤n}, we get E(Yn1{T>n})=−P(T≤n)→−1. Put beside the doubling strategy, where C1 held and C2 failed, this shows that neither of the two hypotheses of version C can be dropped in favour of the other.
Condition B1 is of a different kind. It is not traded against B2 and B3; it is what makes the conclusion mean anything, because if P(T=∞)>0 then YT is undefined on an event of positive probability. The gap is not always repairable by a convention. Take X1,X2,… independent with P(Xi=1)=p<21 and P(Xi=−1)=q=1−p, let Sn=X1+⋯+Xn with S0=0, and let Mn=Sn−(p−q)n, a martingale with M0=0. Let T be the first passage of S to level 1. This walk drifts downwards, and the first-passage probability is P(T<∞)=p/q<1, so {T=∞} has probability 1−p/q>0. Now M is itself a random walk whose steps Xi−(p−q) are bounded, non-degenerate and of mean zero, and such a walk oscillates: almost surely limsupnMn=+∞ and liminfnMn=−∞. Hence on {T=∞} the sequence Mn has no limit, no value can be assigned to MT there by any limiting convention, and the statement E(MT)=E(M0) cannot even be formulated, let alone rescued. B1 is a condition of well-posedness, while B2, B3, C1 and C2 are conditions that make a well-posed conclusion true.
What is never needed. Nothing in part A: for a bounded stopping time the conclusion is unconditional. Every counterexample must therefore use a stopping time that is unbounded, even if it is finite with probability one.
Method(Using optional stopping to compute a hitting quantity)Original
To find a hitting probability, a mean hitting time, or the mean value of some functional at a random time:
Name the filtration and the stopping time. Usually Fn is generated by the observations up to time n, and T is a first passage time, which is automatically a stopping time.
Find a martingale whose value at T encodes the answer. For a hitting probability, look for a function ψ of the state with E(ψ(Xn+1)∣Fn)=ψ(Xn) - a harmonic function for the dynamics - so that ψ(XT) takes one known value on each possible stopping set. For a mean time, look for a compensated quantity such as ψ(Xn)−αn, whose value at T contains T itself.
Verify one sufficient condition. Check E(T)<∞ and a uniform bound c on the conditional absolute increments before T (version C). The most common way to get E(T)<∞ is a geometric tail: exhibit m and ϵ>0 such that, from any history, the event occurs within the next m steps with probability at least ϵ; then P(T>km)≤(1−ϵ)k and so E(T)<∞.
Write E(YT)=E(Y0) and expand the left side over the possible values of YT. One martingale gives one linear equation; if there are two unknowns, find a second martingale.
Sanity-check the answer in a limiting or symmetric case where it is known independently, and - if the hypotheses were only sufficient - against a direct recursion.
Example(How long until HH appears?)Original
A fair coin is tossed repeatedly, giving independent outcomes C1,C2,…, each equal to H or T with probability 21. Let
T=min{n≥2:Cn−1=Cn=H}
be the first time the block HH is complete. We compute E(T).
Step 1: build the martingale. Imagine a casino that pays fair odds: a stake of s on a correctly guessed letter returns 2s, and nothing otherwise. Before toss k, and only while the game is still running (that is, if k≤T), a new gambler Gk arrives, stakes 1 on "Ck=H", and thereafter follows one rule: if she wins, she stakes her entire holding on the next letter of the pattern HH; if she loses, she is broke and plays no more; if she has matched the whole pattern, she keeps her holding and retires.
Let Fn=σ(C1,…,Cn), let Hn be the total holding of all gamblers after toss n (so H0=0), and set
Yn=Hn−(T∧n),
the gamblers' total wealth minus the total 1 per gambler they have paid in.
Y is a martingale. The event {T>n} lies in Fn, since T is a first passage time and hence a stopping time, so the two cases may be treated separately. On {T≤n} no new gambler arrives, no stake is placed and nothing is paid in, so Hn+1=Hn and T∧(n+1)=T∧n: both terms of Y are frozen and the increment is zero. On {T>n} gambler Gn+1 adds 1 to the holdings and 1 to the amount paid in, retired gamblers' holdings are unchanged, and each active gambler's holding is doubled or wiped out with probability 21 each, so its conditional mean is unchanged. Hence
E(Hn+1−Hn∣Fn)=1{T>n}=E((T∧(n+1))−(T∧n)∣Fn),
the two displayed quantities being equal because T∧(n+1)−T∧n is the indicator of {T>n} outright. Subtracting, E(Yn+1∣Fn)=Yn, with Y0=0. (The unconditional identity E(Hn+1∣Fn)=Hn+1 would be wrong: it holds only while the game is still running.)
Step 2: check the hypotheses of version C. At most two gamblers can hold anything at once (one who has matched H, holding 2, and one who has just matched HH, holding 4), so 0≤Hn≤6 and ∣Yn+1−Yn∣≤7 always; condition C2 holds with c=7. For C1, split the tosses into disjoint pairs: each pair is HH with probability 41 independently of the others, so P(T>2m)≤(3/4)m and E(T)<∞. Also E∣Yn∣≤6+n<∞.
Step 3: evaluate at T. At the moment HH is complete, gambler GT−1 guessed both letters and holds 22=4; gambler GT guessed the first letter H correctly (it is the second H of the block) and holds 21=2; every other gambler is broke. So HT=4+2=6, a constant, and T∧T=T. Optional stopping gives
0=E(Y0)=E(YT)=E(HT)−E(T)=6−E(T),
so
E(T)=6
Step 4: the same calculation for HT. Let T′ be the first completion of the block HT. Now GT′−1 has matched both letters and holds 4, but GT′ staked on "CT′=H" and the toss was T, so she is broke: HT′=4 and E(T′)=4. Two patterns of the same length have different mean waiting times, and the difference is caused entirely by whether the pattern overlaps itself - HH does, HT does not. Conditioning on the first step gives the same answers, but needs a separate recursion for each pattern; the martingale reads the answer straight off the overlap structure.
Step 5: check against a direct recursion. Track the length of the current match: let e0 and e1 be the mean remaining number of tosses when the current match has length 0 and 1. For HH,
e0=1+21e1+21e0,e1=1+21⋅0+21e0,
whose solution is e1=4, e0=6, confirming step 3. The corresponding pair for HT gives e0=4.
Step 6: a biased check. If the coin shows H with probability p, fair odds return s/p on a correct guess of H and s/q on a correct guess of T, where q=1−p. The argument is unchanged, and the holdings at T become p−2+p−1 for HH and (pq)−1 for HT, so
E(THH)=p21+p1,E(THT)=pq1.
At p=21 these give 6 and 4 as before. At p=31 they give 9+3=12 and 1/(31⋅32)=4.5, and the recursion of step 5 with p=31 returns e0=12 and e0=4.5 respectively. The ranking of the two patterns is not universal. At p=32 the formulas give 49+23=3.75 for HH against 1/(32⋅31)=4.5 for HT, so with a coin biased towards heads the self-overlapping pattern is the faster of the two.
Application(Anytime-valid testing, and fair games)AdaptedAdapted from Grimmett & Stirzaker (2020), §12.8, Example (12.8.7), p. 568, §12.4, Exercise (12.4.4), p. 553, and §12.5, Example (12.5.5), p. 555
A non-negative martingale rarely gets large, ever. Let (Y,F) be a non-negative martingale with E(Y0)=1 and let a>0 (the interesting case being a>1). Put T=min{n:Yn≥a}, a stopping time. By part A, E(YT∧n)=1 for every n; since Y≥0 we may throw away the part of that expectation on {T>n} and bound the rest below by a:
1=E(YT∧n)≥E(YT1{T≤n})≥aP(T≤n).
For each fixed n the event {T≤n} is {max0≤m≤nYm≥a}, which is the form in which Grimmett and Stirzaker record the bound, with E(Y0)/x on the right, for non-negative supermartingales. Letting n→∞ and using that {T≤n} increases to {T<∞},
P(Yn≥afor some n≥0)=P(T<∞)≤a1.
This is the martingale form of Markov's inequality (Markov, Chebyshev, Jensen and Chernoff inequalities): it bounds the chance of ever reaching the level a, over all time, by what Markov's inequality would give at a single time.
The level must be attained on the left: the event {supnYn≥a} is strictly larger, since a path may climb towards a without ever reaching it, and the two probabilities can differ. (Take Y0=1 and let Y either jump to 0 and stay there or survive, with the surviving value at time n equal to cn=2−2−n and the survival probability from time n chosen as cn/cn+1 to keep the mean at 1. The total survival probability is ∏n≥0cn/cn+1=c0/limcn=21, and on the survival event Yn<2 for every n while supnYn=2; so with a=2 we have P(T<∞)=0 but P(supnYn≥2)=21.) The supremum form of the bound is nevertheless true, by one extra line: for 0<b<a the event {supnYn≥a} is contained in {Yn≥b for some n}, whose probability is at most 1/b, and letting b↑a gives P(supn≥0Yn≥a)≤1/a as well. The example above, with P(supnYn≥2)=21=1/a, shows that this bound is attained.
Why this matters for sequential experiments. Suppose data X1,X2,… arrive one at a time and are independent and identically distributed with density either q (the null hypothesis) or p. Assume q>0 wherever p>0, so that the ratio p/q is defined and finite at every point the data can reach under either hypothesis, and the set {q=0} is null under both. The likelihood ratio
Ln=i=1∏nq(Xi)p(Xi),L0=1,
is a martingale under the null, because the Xi are independent and, almost surely under q,
the cancellation being legitimate exactly because the integral runs over {q>0}, and the last equality holding because p vanishes off that set.
It is also non-negative, so the bound above applies: if the rule is "reject the null the first time Ln≥100", then the probability of ever rejecting a true null is at most 0.01no matter how long one watches, and no matter how the decision to stop is made. This is the formal content of the modern "anytime-valid" or "safe" testing methods, and it is why they are immune to the usual objection to repeatedly peeking at accumulating data: a fixed-sample test loses its guarantee when the stopping rule is data-dependent, and a martingale test does not.
A concrete instance: to detect a unit shift in the mean of standard normal noise, q=N(0,1) and p=N(1,1), so logLn=∑i≤nXi−n/2, and the rule "stop and reject when ∑i≤nXi≥n/2+log100=n/2+4.6052" has false-alarm probability at most 0.01 over an unlimited run. The same quantity, watched under p instead, grows at rate 21 per observation, so the detection delay is about 9.2 observations.
Elsewhere. Martingales are the natural language wherever a quantity is conserved in conditional mean. For a Markov chain with a bounded harmonic function ψ (one satisfying ∑jpijψ(j)=ψ(i)), the sequence ψ(Xn) is a martingale, and optional stopping at a first passage time shows that on an irreducible recurrent chain such a ψ must be constant - the discrete Liouville theorem, which is Grimmett and Stirzaker's Example (12.5.5); see Discrete-time Markov chains. In continuous time the Wiener process and its compensated square W(t)2−t are martingales, and the stochastic integral of The Ito integral and Ito's formula is built so as to be one; the pricing of a financial derivative, and the "fluctuation-dissipation" bookkeeping of a driven physical system, are both statements that a suitably discounted quantity is a fair game.
Builds the Ito integral against a Wiener process by extending the elementary sum over predictable step functions along the Ito isometry, shows that the resulting integral is a martingale with second moment equal to the mean integrated square of the integrand, and states and sketches Ito's formula, the chain rule with its extra second-derivative term.
Note(What goes wrong with the ordinary chain rule)Original
Write δW=Wt+δt−Wt for an increment of a standard Wiener process. The whole of this unit turns on one fact about that increment: it is of size δt, not δt, because E[(δW)2]=δt.
In ordinary calculus a second-order Taylor term is negligible because (δx)2 is small compared with δx. Here (δW)2 is comparable with δt, so the second-order term of a Taylor expansion survives the limit and contributes a drift. Everything below is the bookkeeping that makes this precise: a definition of ∫ψdW that is stable under the limit, and a chain rule that carries the extra term.
A reader who wants only the operational rule can take away one line, to be used inside any smooth function of an Itô process: (dW)2=dt, dWdt=0, (dt)2=0.
Definition(Admissible integrands, predictable step functions, and the integral)AdaptedAdapted from Grimmett & Stirzaker (2020), §13.8, pp. 605-608, including (13.8.1), (13.8.2), (13.8.3) and the definition of $I(\psi)$
Let W be a standard Wiener process on a complete probability space (Ω,F,P). For t≥0 let Ft be the smallest sub-σ-field of F that contains every null event and with respect to which each of Ws, 0≤s≤t, is measurable, and write F={Ft:t≥0} for the filtration so obtained.
A process ψ is measurable if the map (t,ω)↦ψt(ω) is measurable for the product σ-field B⊗F, where B is the Borel σ-field on [0,∞); a measurable ψ is adapted if ψt is Ft-measurable for every t. Let
With ∥⋅∥ as its norm, A is a Hilbert space, and therefore complete. (Strictly, ∥⋅∥ is a norm on the equivalence classes of processes that agree with probability one, since ∥ψ∥=0 only forces ψ to vanish almost everywhere.)
A predictable step function is a process of the form
where each Cj has a finite second moment and is Faj-measurable — the height on an interval is decided before the interval begins. Its integral is the obvious sum,
I(φ)=j=0∑n−1Cj(Waj+1−Waj),(3)
and I is linear on predictable step functions.
For general ψ∈A the Itô integralI(ψ) is the mean-square limit of I(φ(n)) along any sequence {φ(n)} of predictable step functions with ∥φ(n)−ψ∥→0. The statement below says that such a sequence exists and that the limit does not depend on which one is used. We write
∫0∞ψsdWs:=I(ψ),∫0tψsdWs:=I(ψ1(0,t]).
Finally, a continuous process X is an Itô process, or a diffusion, with infinitesimal mean μ(t,x) and infinitesimal variance σ(t,x)2 if
Xt=X0+∫0tμ(s,Xs)ds+∫0tσ(s,Xs)dWs,(4)
for which the shorthand is dXt=μ(t,Xt)dt+σ(t,Xt)dWt. Grimmett and Stirzaker state that the existence of such a process is beyond the scope of their account; see Stochastic differential equations.
Theorem(The Ito isometry, and the integral it defines)AdaptedAdapted from Grimmett & Stirzaker (2020), §13.8, Lemma (13.8.4) and Theorem (13.8.6), p. 606, the uniqueness argument on p. 608, Theorem (13.8.11), p. 608, Theorem (13.8.13), p. 609, and Exercises (13.8.4) and (13.8.6), p. 610
Let W be a standard Wiener process and F the filtration above. Write ∥U∥2=E∣U∣2 for a random variable U.
Isometry on step functions. If φ is a predictable step function then ∥I(φ)∥2=∥φ∥.
Approximation. If ψ∈A has continuous sample paths, there is a sequence {φ(n)} of predictable step functions with ∥φ(n)−ψ∥→0. (Grimmett and Stirzaker prove this under the continuity assumption and record that continuity is not necessary.)
The limit exists and is unique. For such a ψ the random variables I(φ(n)) converge in mean square, and any two approximating sequences give limits that agree with probability one. This limit is I(ψ).
The isometry persists. For an adapted ψ with ∥ψ∥<∞, ∥I(ψ)∥2=∥ψ∥. The source also records the bilinear version
E(I(ψ1)I(ψ2))=E(∫0∞ψ1(t)ψ2(t)dt),
whenever both sides exist, for ψ1 and ψ2predictable step functions. It extends to ψ1,ψ2∈A at once by polarisation, since I is linear and the norm identity holds throughout A: four applications of ∥I(⋅)∥2=∥⋅∥ to ψ1±ψ2 give it. (That last deduction is ours, not the source's.) In particular EI(ψ)=0.
Martingale property. If ψ∈A has continuous sample paths then Jt=∫0tψudWu is a martingale with respect to F.
Converse (martingale representation). If M is a martingale with respect to F, there is an adapted process ψ with Mt=M0+∫0tψudWu for all t≥0. Grimmett and Stirzaker state this without proof.
Item 4 is the identity usually called the Itô isometry. It is the workhorse: it turns a question about the second moment of a stochastic integral into an ordinary deterministic integral of a second moment.
Theorem(Ito's formula)AdaptedAdapted from Grimmett & Stirzaker (2020), §13.9, Theorem (13.9.2), p. 610, and Theorem (13.9.4), p. 611
General form. Let dX=μ(t,X)dt+σ(t,X)dW and let Yt=f(t,Xt), where f:[0,∞)×R→R is twice continuously differentiable. Then Y is again an Itô process, and
Grimmett and Stirzaker state (5) with the proof omitted, and record that it extends to multivariate diffusions.
The case X=W. Under the name Itô's simple formula, the source states the following. Suppose W is standard, and let the function f(s,w) possess continuous partial derivatives up to the third order on the half-plane s≥0. Then Yt=f(t,Wt) is an Itô process, and
dY=[ft(t,W)+21fww(t,W)]dt+fw(t,W)dW.(6)
This is (5) with μ=0 and σ=1. Grimmett and Stirzaker give a sketch proof of (6) and note that its differentiability assumption may be weakened; the extra derivative is an artefact of that proof, not of the result. The term 21fww, absent from the Newton-Leibniz chain rule, is the Itô correction.
Idea of the proofOriginal
The integral is built by isometric extension, exactly as the Fourier transform is extended from nice functions to L2. For a step function whose heights are fixed before each interval starts, the sum (3) is a sum of uncorrelated terms, so its second moment is the sum of the individual second moments — which is precisely ∥φ∥2. So I is an isometry from a subset of the Hilbert space A into L2(Ω). An isometry is uniformly continuous, so it carries Cauchy sequences to Cauchy sequences; since both spaces are complete, I extends uniquely and isometrically to the closure of the step functions, and the remaining work is to show that this closure is all of A.
The martingale property is the same cancellation looked at from one step back. Condition on Fs. Each term Cj(Waj+1−Waj) of the sum that lies wholly after s contributes nothing, because the increment is independent of Faj and has zero mean while Cj is known at time aj; each term wholly before s is already known. This is exactly where adaptedness is spent.
Itô's formula is a second-order Taylor expansion in which the second-order term does not die. Expand f over a short interval: the first-order terms reproduce the ordinary chain rule, and the second-order term in w carries (δW)2. That quantity has mean δt and fluctuations small enough that, summed over t/δt intervals, it may be replaced by δt in the limit. Third-order terms carry (δW)3∼(δt)3/2 and do vanish. So one extra term survives, and only one.
Proof(The isometry, and the construction of the integral)AdaptedAdapted from Grimmett & Stirzaker (2020), §13.8, proof of Lemma (13.8.4), p. 606, and the argument following Theorem (13.8.6), pp. 607-608
We prove items 1 and 3, and then item 5; items 2, 4 and 6 are quoted from the source as stated.
Step 1: square the elementary sum. For a predictable step function φ as in (2), write Δj=Waj+1−Waj. Multiplying out,
Step 2: the diagonal terms. Fix j and condition on Faj. Because φ is a predictable step function, Cj is Faj-measurable, so it comes out of the conditional expectation; and Δj is independent of Faj with E(Δj2)=aj+1−aj, since W has independent increments of variance equal to their length. By the tower property,
Step 3: the off-diagonal terms vanish. Let j<k. Condition on Fak. The three factors Cj, Ck and Δj are all Fak-measurable, again by predictability and because aj<aj+1≤ak; the remaining factor Δk is independent of Fak and has mean zero. Hence each such term has expectation 0.
Step 4: recognise the answer. Combining Steps 1-3,
the middle equality because φ takes the constant value Cj on an interval of length aj+1−aj. This is item 1.
Step 5: the approximating integrals form a Cauchy sequence. Let ψ∈A have continuous sample paths and take {φ(n)} from item 2. The difference φ(m)−φ(n) is again a predictable step function, so by linearity of I, by Step 4 applied to that difference, and by the triangle inequality in A,
as m,n→∞. The space L2(Ω) is complete, so I(φ(n)) converges in mean square to some random variable.
Step 6: the limit does not depend on the sequence. Let {ρ(n)} be a second sequence of predictable step functions with ∥ρ(n)−ψ∥→0, with limit I(ρ). Insert both approximations:
The outer two terms tend to 0 by the definition of the limits. The middle term equals ∥φ(n)−ρ(n)∥ by Step 4 and linearity, and that is at most ∥φ(n)−ψ∥+∥ρ(n)−ψ∥→0. So the left-hand side, which does not depend on n, is 0: the two limits agree with probability one. This is item 3.
Step 7: the martingale property. Let 0<s<t and choose a partition 0=a0<⋯<am=s<⋯<an=t that includes s as a partition point. Let ψu(n)=∑jψaj1(aj,aj+1](u) and Jv(n)=∫0vψu(n)dWu. Conditioning the sum (3) on Fs term by term:
a term with j<m lies entirely in [0,s], and is Fs-measurable, so it is unchanged;
a term with j≥m has aj≥s; conditioning first on Faj and using that ψaj is Faj-measurable while Waj+1−Waj is independent of Faj with zero mean, the term has conditional expectation 0.
Hence E(Jt(n)∣Fs)=Js(n). Now refine the partition so that its mesh tends to 0. By the approximation of item 2, ∥ψ(n)1(0,v]−ψ1(0,v]∥→0 for v=s and v=t, so Js(n)→Js and Jt(n)→Jt in mean square. Mean-square convergence is preserved by conditional expectation, so the identity passes to the limit: E(Jt∣Fs)=Js almost surely. □
Proof(Sketch of Ito's simple formula)AdaptedAdapted from Grimmett & Stirzaker (2020), §13.9, sketch proof of Theorem (13.9.4), p. 611, with Exercise (13.7.2), p. 604
This is a sketch; the source itself gives only a summary, and the estimates behind items (i), (ii) and (iv) below are omitted there and here. Fix t>0, put δ=t/n, write tj=jδ and Δj=Wtj+1−Wtj.
Step 1: telescope, splitting each step into a time move and a space move. Write the total increment as
Each bracket in the first sum moves only in time, each in the second only in space; the two sums are handled by different arguments, which is the reason for the split.
Step 2: Taylor-expand each sum. By the mean value theorem there are random θj∈[tj,tj+1] making the first sum ∑jft(θj,Wtj+1)δ. By Taylor's theorem with Lagrange remainder in the space variable there are random Ξj between Wtj and Wtj+1 making the second sum
The third derivative appears here, and nowhere else; this is why the hypothesis of (6) asks for three derivatives.
Step 3: the time sum becomes an ordinary integral.W has continuous sample paths and ft is continuous, so along almost every path the integrand is continuous and the sum is a Riemann sum: ∑jft(θj,Wtj+1)δ→∫0tft(s,Ws)ds almost surely.
Step 4: the first space sum becomes the Itô integral. The summand is the elementary integral (3) of the predictable step function with heights fw(tj,Wtj), which is exactly the approximating step function of the construction above; it converges in mean square to ∫0tfw(s,Ws)dWs.
Step 5: the quadratic sum becomes a time integral — the Itô correction. Here is the one genuinely stochastic step. Each Δj2 has mean δ, and Δj2 and Δk2 are independent for j=k because W has independent increments. Replacing Δj2 by its mean δ therefore introduces an error that is a sum of n independent-in-the-relevant-sense, mean-zero terms each of size δ; the variance of that error is of order nδ2=tδ→0. (For fww≡1 this is the statement that ∑jΔj2→t in mean square, the quadratic variation of the Wiener process.) Hence
in mean square, and the second of these is a Riemann sum converging to 21∫0tfww(s,Ws)ds.
Step 6: the cubic sum is negligible.E∣Δj∣3 is of order δ3/2, so n such terms total order nδ3/2=tδ1/2→0; the source states that this term converges to zero in mean square.
Step 7: assemble. Adding the limits of Steps 3-6 gives the integrated form of (6), and t was arbitrary. □
Remark(Where the hypotheses are used)Original
Predictability (the height Cj is Faj-measurable) is what makes the isometry true, and it is used twice, in Steps 2 and 3. Drop it and both fail. Take the single-interval step function with a0=0, a1=t and height C0=Wt, which is Ft-measurable but not F0-measurable. Then I(φ)=Wt(Wt−W0)=Wt2, whose mean is t, not 0, and whose second moment is E(Wt4)=3t2, while ∥φ∥2=E∫0tWt2ds=t2. Both item 4 and the martingale property are destroyed.
This is the same thing as the dependence of the limit on the evaluation point. Partition (0,t] into n equal pieces and form ∑jWθj(Wtj+1−Wtj). Taking θj=tj, the left endpoint, gives 21(Wt2−t) in mean square — the Itô value; taking θj=tj+1 gives 21(Wt2+t). A Riemann integral would not notice the difference. The Wiener process has unbounded total variation on every non-trivial interval, and that is exactly why it does. The midpoint rule gives a third answer and a different calculus, discussed in Stochastic differential equations.
The square-integrability condition (1) is what makes Jt a martingale rather than only a local martingale. Grimmett and Stirzaker record in a footnote that without it one obtains a local martingale, and that in a more general theory the correct condition on the integrand is not adaptedness but local boundedness and predictability, with W replaced by a semimartingale.
Continuity of the sample paths of ψ is used only in item 2 and in item 5, and only because the proof given for the density of step functions uses uniform continuity on compact time intervals. The source states that this condition is not necessary.
Two derivatives in Itô's formula (5) are not a technicality. The term 21fxxσ2 is a leading term, not a remainder. With f only once continuously differentiable there is no formula of this shape at all, because the first-order Taylor remainder is of size (δW)2∼δt and does not vanish on summation. The simplest instance is f(w)=w2: the Newton-Leibniz rule would give d(W2)=2WdW and hence ∫0tWsdWs=21Wt2, whereas the correct value is 21(Wt2−t). The two differ by a quantity, t/2, that is not small.
The simple formula (6) asks for three derivatives where the general formula (5) asks for two. This looks backwards, since (6) is a special case of (5). It is an honest feature of the source: (5) is quoted there without proof, while (6) is given the sketch proof above, and that sketch spends a third derivative on the Taylor remainder in Step 2. Grimmett and Stirzaker say explicitly that the assumption may be weakened. Use (5) when two derivatives are all you have.
Method(Evaluating a stochastic integral by reversing Ito's formula)Original
There is no fundamental theorem of calculus for ∫ψdW with a general adapted ψ. When the integrand has the form g(s,Ws), however, the following works.
Find an antiderivative in the space variable. Choose f(s,w) with fw(s,w)=g(s,w), three times continuously differentiable (or use (5) if you can justify two).
Compute the driftft+21fww.
Write down Itô's formula in integrated form.f(t,Wt)−f(0,W0)=∫0t[ft(s,Ws)+21fww(s,Ws)]ds+∫0tg(s,Ws)dWs.
Rearrange to isolate the stochastic integral. If the drift in Step 2 vanishes identically, the stochastic integral is simply f(t,Wt)−f(0,W0) and is a martingale.
Check admissibility before using the isometry or the martingale property. Both need g(⋅,W⋅)1(0,t]∈A, that is ∫0tE(g(s,Ws)2)ds<∞. For a polynomial g this is automatic, because all moments of a normal variable are finite.
Take second moments with the isometry.E[(∫0tgdW)2]=∫0tE(g(s,Ws)2)ds, an ordinary one-dimensional integral. Step 1 has already given the first moment, which is 0.
Check the answer independently whenever the left-hand side has a known distribution; see the example.
Example(A cubic martingale, and its second moment two ways)Original
Take f(t,w)=w3−3tw, a polynomial and so as smooth as required, and let W be a standard Wiener process with W0=0.
Step 1: the derivatives.
ft=−3w,fw=3w2−3t,fww=6w.
Step 2: the drift cancels.
ft+21fww=−3w+21(6w)=0.
So by Itô's simple formula (6), Yt=Wt3−3tWt has dY=(3Wt2−3t)dWt: a pure stochastic integral, with no drift.
Step 3: read off the integral. Since f(0,0)=0,
∫0t(3Ws2−3s)dWs=Wt3−3tWt.
Step 4: check admissibility. With g(s,w)=3w2−3s and Ws∼N(0,s),
using E(Ws2)=s and E(Ws4)=3s2. This is integrable on [0,t], so g(⋅,W⋅)1(0,t]∈A and items 4 and 5 of the statement apply. In particular Y is a martingale, so E(Wt3−3tWt)=0 for every t.
Step 5: second moment by the isometry.
E[(Wt3−3tWt)2]=∫0t18s2ds=6t3.
Step 6: the same number directly. Expand the square and use the moments of Wt∼N(0,t), namely E(Wt2)=t, E(Wt4)=3t2 and E(Wt6)=15t3:
E(Wt6−6tWt4+9t2Wt2)=15t3−18t3+9t3=6t3.
The two routes agree, which checks the isometry on this example.
A number. At t=2, E[(W23−6W2)2]=6⋅23=48.
Note how badly the naive chain rule would do here. Treating dW as an ordinary differential would give ∫0t3Ws2dWs=Wt3, and Wt3 is not a martingale in any useful sense — it has the right mean, 0, but the compensating −3tWt is exactly what the Itô correction supplies.
Application(Multiplicative noise and the gap between the mean and the typical path)AdaptedAdapted from Grimmett & Stirzaker (2020), §13.9, Example (13.9.9), p. 612, and §13.10, equation (13.10.2), p. 614
A quantity that grows by random proportional increments — a population, the size of an aerosol droplet under random coagulation, the price of a traded asset — is modelled by geometric Brownian motion,
dSt=St(μdt+σdWt).
Apply Itô's simple formula to Yt=exp(μt+σWt): with f(t,w)=eμt+σw one has ft=μf, fw=σf, fww=σ2f, so dY=(μ+21σ2)Ydt+σYdW. Therefore the solution of the displayed equation with S0=1 is not eμt+σWt but
St=exp[(μ−21σ2)t+σWt],
which is the form in which the Black-Scholes model of §13.10 of the source states it, σ there being called the volatility.
The consequence is one that is easy to get wrong in practice. Since EeσWt=eσ2t/2, the mean grows at rate μ:
E(St)=eμt.
But logSt is normal with mean (μ−21σ2)t, so the median, and the almost sure long-run growth rate of a single realisation, is governed by μ−21σ2. The two differ by exactly the Itô correction.
Take μ=0.05 and σ=0.4 per year. Then μ−21σ2=0.05−0.08=−0.03. After ten years the mean of the ensemble is e0.5=1.6487 times its initial value, while the median path is at e−0.3=0.7408 of it, and
P(S10<S0)=Φ(0.4100.03×10)=Φ(0.2372)=0.5937,
so a clear majority of realisations have gone down while the ensemble mean has gone up by nearly two thirds. The mean is carried by a thin tail of very large values. Whenever a quantity fluctuates multiplicatively, the average over an ensemble and the behaviour of a typical member are different questions, and Itô's formula is what converts between them. The same term reappears as the drift of the Fokker-Planck equation for logS; see The Fokker-Planck equation.
Gives a recipe for solving a linear stochastic differential equation with additive noise by an integrating factor, proves that its solution is Gaussian with an explicit mean, variance and autocovariance, applies this to the Ornstein-Uhlenbeck process and the Langevin equation, and explains how the Ito and Stratonovich readings of the same equation differ and how to convert between them.
Note(What an equation driven by noise can mean)Original
A physicist writes mv˙=−γv+ξ(t) and calls ξ "white noise". Nothing in ordinary analysis supports this: the intended ξ is the derivative of a Wiener process, and Brownian motion has nowhere differentiable paths. The repair is to multiply through by dt and read the equation as an integral equation in which the noise appears only inside a stochastic integral, never differentiated.
Once that is done, a second question appears that has no analogue in ordinary differential equations: which stochastic integral? The Itô integral and the Stratonovich integral give different answers to the same written equation. This unit gives the method for the case that is solvable in closed form — linear drift and additive noise — and then says what the two readings are and how to pass between them.
Definition(Stochastic differential equations, and the two readings)AdaptedAdapted from Grimmett & Stirzaker (2020), §13.8, equations (13.8.9)-(13.8.10), p. 608, §13.9, equation (13.9.1), p. 610, and §13.7, pp. 603-604
Let W be a standard Wiener process, with the filtration F built in The Ito integral and Ito's formula. A stochastic differential equation (an SDE) is written
dXt=μ(t,Xt)dt+σ(t,Xt)dWt(1)
and is shorthand for the integral equation
Xt=X0+∫0tμ(s,Xs)ds+∫0tσ(s,Xs)dWs,t≥0,(2)
in which the second integral is an Itô integral. A continuous adapted process X satisfying (2) is called a solution, and is an Itô process, or diffusion, with infinitesimal mean μ and infinitesimal variance σ2.
The last integral in (2) is the mean-square limit of sums ∑jσ(tj,Xtj)(Wtj+1−Wtj) in which the integrand is evaluated at the left endpoint of each partition interval. Evaluating instead at the midpoint produces a different limit, and that limit defines the Stratonovich integral, written ∫0tσ∘dW. An equation written with ∘dW in place of dW,
dXt=μ~(t,Xt)dt+σ(t,Xt)∘dWt,(3)
is a Stratonovich SDE: a different equation, with different solutions, not a different notation for (1). The remark below gives the dictionary between them.
Throughout, linear with additive noise means the special case
dXt=(α(t)Xt+ν(t))dt+δ(t)dWt,(4)
in which α,ν,δ are deterministic functions of t alone. (The letter ν is used for the affine forcing because β is reserved below for the relaxation rate of the Ornstein-Uhlenbeck process.) Because σ=δ does not depend on x, equations (1) and (3) coincide for (4) — the two readings differ only when the noise is multiplicative.
Method(Existence and uniqueness: what the source asserts)AdaptedAdapted from Grimmett & Stirzaker (2020), §13.9, opening paragraph, p. 610, and §13.12, Problem (13.12.21), p. 629
Grimmett and Stirzaker assert the following about equation (1). Provided μ, σ and X0 satisfy conditions that they describe as rather weak, there is exactly one solution; it is a Markov process, and its sample paths are continuous. They do not list the conditions, and they declare the proof beyond the scope of their book. We report the assertion in that form rather than substituting the Lipschitz and linear-growth hypotheses familiar from other accounts, which are not the hypotheses the source gives.
That some condition is needed is shown by their Problem (13.12.21), Tanaka's example. Let sign(x)=1 for x>0 and −1 otherwise, and consider
dXt=sign(Xt)dWt,X0=0.
The process Vt=∫0tsign(Ws)dWs is itself a standard Wiener process and solves this equation; so does −V. (The source writes the same letter W for the process used to build V and for the process driving the equation; the driving process is the one reconstructed from V itself, namely ∫0⋅sign(Vs)dVs, which is again a standard Wiener process. Nothing below depends on which reading is taken.) Solutions are therefore not unique path by path, although every solution has the finite-dimensional distributions of a Wiener process. The coefficient sign is bounded, measurable and non-zero, but discontinuous at the one point the process keeps returning to.
For the linear equation (4) there is no such difficulty, and the next statement settles existence, uniqueness and the law of the solution outright.
Method(Linear SDE with additive noise: closed form and law)Original
Let ν and δ be continuous real functions on [0,∞), let α be a continuously differentiable real function on [0,∞), let X0 be a constant, and set
Φt=exp(−∫0tα(r)dr).
Then equation (4) has one and only one continuous adapted solution, namely
Xt=Φt−1(X0+∫0tΦsν(s)ds+∫0tΦsδ(s)dWs).(5)
The solution is a Gaussian process, with the usual convention that a constant counts as normal with variance 0 (if δ≡0 the equation has no noise and the solution is the deterministic function (5)), and
A random initial value. Suppose instead that X0 is a random variable with finite variance, independent of W. Then no process starting at X0 is adapted to F at all, because F0 contains only the null events and their complements and so cannot carry a non-degenerate X0. The filtration must be enlarged to G, where Gt=Ft∨σ(X0). This costs nothing: since X0 is independent of W, the increments of W after time t are still independent of Gt, so W is a Wiener process for G too and the Itô integral in (5) is the same object as before. With "adapted" read as adapted to G, formula (5) is again the unique continuous adapted solution, and (6) and (7) hold conditionally on X0. Unconditionally, var(X0) is added to the integral in each of (6) and (7):
for s≤t (both follow from (5) and the independence of X0 and W, which makes the two contributions uncorrelated), and X is a Gaussian process whenever X0 is normal.
This is our own statement; it is proved in full below, and the special case α≡−β, ν≡0, δ≡σ recovers the law of the Ornstein-Uhlenbeck process that the source states in Problem (13.12.4).
Method(The Ornstein-Uhlenbeck process)AdaptedAdapted from Grimmett & Stirzaker (2020), §13.3, Example (13.3.11), p. 584, the representations on p. 603, Exercise (13.8.5), p. 610, and §13.12, Problems (13.12.3) and (13.12.4), p. 627
Fix β>0 and σ>0. The Ornstein-Uhlenbeck process is the diffusion V with instantaneous mean and variance
a(t,x)=−βx,b(t,x)=σ2,
that is, the solution of
dVt=−βVtdt+σdWt.(8)
Items 1-4, and the two representations displayed in item 5, are what the source records; the covariance computation in item 5 and the identifications of σ that it yields are ours.
Law at time t. If V0=u then
Vt∼N(ue−βt,2βσ2(1−e−2βt)).
Limit law.Vt is asymptotically N(0,21σ2/β) as t→∞.
Characterisation. Among processes that are at once Gaussian, Markov, stationary and equipped with an autocovariance function that is continuous, this one is the only one.
Representations. With W standard, each of
U1(t)=e−βtW(e2βt−1),U3(t)=Wt−β∫0te−β(t−s)Wsds
is an Ornstein-Uhlenbeck process, but not with the same σ; the source does not say which σ, and the rest of this item, which settles it, is ours. For s≤t, cov(U1(s),U1(t))=e−β(s+t)min(e2βs−1,e2βt−1)=e−β(t−s)−e−β(t+s), which by (7) is the case σ2=2β, V0=0; U3 is the case σ=1, V0=0, and is what integrating (8) with those values produces. Because these are built from W by continuous operations, the sample paths of V inherit continuity from those of W.
Idea of the proofOriginal
The integrating factor of ordinary linear theory still works, and Itô's formula says why nothing extra appears. For the deterministic equation x˙=αx+ν one multiplies by the integrating factorΦt=e−∫α, which turns the terms in x and x˙ into the derivative of the single product Φx. Here we must apply Itô's formula, not the product rule, to ΦtXt. Itô's formula adds a term 21fxxσ2; but f(t,x)=Φtx is linear in x, so fxx=0 and the correction is absent. That is the whole reason the classical recipe survives for additive noise, and the whole reason it fails for multiplicative noise.
Once the equation is integrated, Gaussianity is free. What is left on the right is a deterministic function plus an Itô integral of a deterministic integrand. Such an integral is, at each stage of its construction, a linear combination of independent normal increments of W, hence normal; and a mean-square limit of normals is normal. Its mean is zero and its variance is given by the Itô isometry, which here is an ordinary integral.
For the Ornstein-Uhlenbeck process the two terms in (5) have a clean physical reading. The factor e−βt damps the memory of the initial value; the stochastic integral accumulates kicks, each discounted by e−β(t−s) according to how long ago it arrived. The variance is a balance between the two, and settles at σ2/(2β): more noise raises it, more friction lowers it.
Proof(Of the linear statement)Original
Write Ft for the filtration of W, and assume throughout that X0 is constant. For a random X0 independent of W, every step below is read with F replaced by the enlarged filtration G of the statement, which changes nothing in Steps 1-4 and 6, because W is a Wiener process for G as well and the Itô integral of a deterministic integrand is built from the same increments; Step 5 is the one step that needs X0 to be constant (or normal).
Step 1: the integrating factor is smooth and deterministic.α is continuously differentiable, so t↦∫0tα is twice continuously differentiable and Φ is a positive, twice continuously differentiable, deterministic function with Φt′=−α(t)Φt, Φt′′=(α(t)2−α′(t))Φt and Φ0=1. The same holds for Φ−1, with (Φ−1)′=αΦ−1 and (Φ−1)′′=(α2+α′)Φ−1. This is the only place α′ is wanted, and the remark below says why it is wanted at all.
Step 2: apply Itô's formula to the product. Suppose X solves (4) and put f(t,x)=Φtx, so that ft=Φt′x=−α(t)Φtx, fx=Φt and fxx=0. By Step 1 the function f is twice continuously differentiable on [0,∞)×R — ftt=Φt′′x exists and is continuous because α∈C1, and ftx=Φt′, fxx=0 — so Itô's formula in the general form of The Ito integral and Ito's formula, which asks for exactly that, applies, with μ(t,x)=α(t)x+ν(t) and σ(t,x)=δ(t):
The two terms containing Xt cancel, which is what the factor was chosen for, and the Itô correction is zero because f is linear in x.
Step 3: integrate and solve. In integrated form, ΦtXt−X0=∫0tΦsνds+∫0tΦsδdWs. The Itô integral on the right is defined because Φδ is continuous and deterministic, so ∫0tΦs2δ(s)2ds<∞ and the admissibility condition of The Ito integral and Ito's formula holds. Dividing by Φt>0 gives (5). Conversely — and this is a second application of Itô's formula, not Step 2 run backwards — set
Yt=X0+∫0tΦsν(s)ds+∫0tΦsδ(s)dWs,
a continuous adapted Itô process with dYt=Φtν(t)dt+Φtδ(t)dWt, and apply Itô's formula to g(t,y)=Φt−1y along Y, which is again twice continuously differentiable by Step 1. Again gyy=0, so there is no Itô correction, and since ∂t(Φt−1)=α(t)Φt−1,
where Xt=Φt−1Yt is exactly the process (5). So (5) is continuous, adapted and solves (4): a solution exists.
Step 4: uniqueness. Let X and X~ both solve (4) with the same X0 and put Z=X−X~. Subtracting the two copies of (2), the stochastic integrals cancel — they do not involve the solution, since δ does not depend on x — and Zt=∫0tα(s)Zsds for every t, pathwise, with Z0=0. Fix T and let A=max[0,T]∣α∣, finite by continuity. Then ∣Zt∣≤A∫0t∣Zs∣ds on [0,T]. Both solutions are continuous, so K=max[0,T]∣Z∣ is finite, and feeding the bound into itself gives
∣Zt∣≤Kn!(At)n,0≤t≤T,
for every n≥0: the case n=0 is the definition of K, and if the bound holds for n then ∣Zt∣≤A∫0tK(As)n/n!ds=K(At)n+1/(n+1)!. Letting n→∞ gives Zt=0. (This is the zero-constant case of Gronwall's inequality, which Dependence on initial data and Gronwall's inequality proves in general; the iteration is spelled out here so that the step needs nothing beyond calculus.) Hence Z≡0 on [0,T] for every T.
Step 5: the solution is Gaussian. Only the term Mt=∫0tΦsδ(s)dWs is random. By construction it is the mean-square limit of the elementary sums ∑jcj(Wtj+1−Wtj) with deterministic coefficients cj, and such a sum is a linear combination of independent normal variables, hence normal with mean 0. Mean-square convergence implies convergence in distribution, and the normal family is closed under limits in distribution (move 3 of the method of Brownian motion), so Mt is normal; its variance is limnvn, since mean-square convergence carries second moments with it. The same argument applied to a finite linear combination ∑kλkMtk shows every such combination is normal, which is to say that (Mt1,…,Mtn) is multivariate normal; so M, and therefore X, is a Gaussian process.
Step 6: mean, variance and autocovariance. Taking expectations in (5), and using E(Mt)=0, gives the first half of (6). The Itô isometry gives
var(Mt)=E(Mt2)=∫0tΦs2δ(s)2ds,
and dividing by Φt2 gives the second half of (6). For s≤t, apply the bilinear form of the isometry to ψ1=Φδ1(0,s] and ψ2=Φδ1(0,t], whose pointwise product is Φ2δ21(0,s]:
cov(Ms,Mt)=E(MsMt)=∫0sΦu2δ(u)2du.
Dividing by ΦsΦt gives (7). □
Deduction of the Ornstein-Uhlenbeck law. Put α≡−β, ν≡0, δ≡σ and V0=u. Then Φt=eβt, and (5) reads Vt=ue−βt+σ∫0te−β(t−s)dWs. By (6) the mean is ue−βt and the variance is
as s,t→∞ with t−s fixed, which is the stationary autocovariance of item 3.
Remark(Where the hypotheses are used)Original
Additive noise (δ independent of x) is used in Steps 2 and 4, and nowhere can it be dropped. In Step 2 it is what makes fxxσ2=0; in Step 4 it is what makes the stochastic integrals cancel between two solutions. For multiplicative noise the integrating factor must itself be corrected: the equation dS=S(μdt+σdW) is solved not by exp(μt+σWt) but by exp((μ−21σ2)t+σWt), as Grimmett and Stirzaker (Grimmett & Stirzaker (2020), §13.9, Example (13.9.9), p. 612) show. The missing −21σ2t is exactly the term that Step 2 was able to discard.
Deterministic coefficients are used only in Step 5 — but Steps 1-2 have to be reread. If α,ν,δ are adapted processes with suitable integrability then (5) is still the solution and uniqueness still holds, but the route to it changes: Φt=exp(−∫0tα) is then a continuous adapted process, not a deterministic function, so Step 1 gives it as a process of finite variation rather than a differentiable function, and Step 2 cannot apply Itô's formula to a deterministic f. One uses instead the product rule for a finite-variation factor times an Itô process, d(ΦtXt)=ΦtdXt+XtdΦt, in which there is no cross-variation term precisely because Φ has finite variation; the two terms in Xt cancel as before. Step 4 goes through with A=max[0,T]∣α∣ replaced by a pathwise bound. What does not survive is Step 5: Mt is no longer a combination of normals, and nothing in (6)-(7) holds.
Continuous differentiability of α is used in Step 1 only, and only to make f(t,x)=Φtx twice continuously differentiable in both variables, which is what the general Itô formula of The Ito integral and Ito's formula demands. Continuity of α alone makes Φ merely C1, so ftt need not exist; the conclusions are unchanged under that weaker hypothesis, but reaching them needs the C1,2 form of Itô's formula (once continuously differentiable in t, twice in x), which this handbook does not state, so we assume α∈C1 rather than appeal to it. Mere continuity of α is all that Step 4 wants, where it bounds ∣α∣ on a compact interval and so makes the iterated bound converge; local integrability and a local bound would suffice there.
∫0tΦs2δ(s)2ds<∞ is the admissibility condition of the Itô integral, used in Step 3. Without it the stochastic integral in (5) need not exist and (6) is meaningless.
β>0 in the Ornstein-Uhlenbeck process is what makes items 1-4 of that statement true. At β=0 equation (8) is dV=σdW, whose solution is a Wiener process: the variance σ2t grows without bound, there is no limit law and no stationary version. For β<0 the mean ue−βt and the variance σ2(e2∣β∣t−1)/(2∣β∣) both diverge. Only a restoring drift produces a balance.
The existence and uniqueness assertion for the general equation (1) rests on conditions the source does not state. Tanaka's example above shows that the pathwise part of the conclusion genuinely fails for a discontinuous σ, so "rather weak" cannot be read as "none".
Remark(Ito versus Stratonovich)Original
Partition (0,t] into n equal intervals with endpoints tj=jt/n and write Δj=Wtj+1−Wtj. For the integrand W itself, the three natural evaluation points give three different limits in mean square:
evaluation point
limit of ∑jWθjΔj
left, θj=tj (Itô)
21(Wt2−t)
midpoint, θj=21(tj+tj+1) (Stratonovich)
21Wt2
right, θj=tj+1
21(Wt2+t)
The midpoint row is checked as follows. Split Δj=Aj+Bj into its halves, which are independent and each N(0,t/(2n)). Then Wθj=21(Wtj+Wtj+1)+21(Aj−Bj). The first part telescopes exactly:
The second part contributes ∑j21(Aj−Bj)(Aj+Bj)=21∑j(Aj2−Bj2), a sum of n independent mean-zero terms each of variance of order (t/n)2, so of total variance of order t2/n→0. Hence the midpoint limit is 21Wt2 in mean square. A note on the source: the prose of Grimmett and Stirzaker (Grimmett & Stirzaker (2020), §13.7, p. 604, the paragraph after (13.7.5)) names the midpoint value as Wt2; their own Exercise (13.7.3)(c) on the same page, which asks for the limit of ∑j21(Vj+1+Vj)Δj, telescopes to 21Wt2, and that is the value consistent with the two flanking entries of the table. We use 21Wt2.
The dictionary. Let X solve dX=μdt+σ(t,X)dW with σ continuously differentiable in x, and let the integrand be h(Xs) for a continuously differentiable h. Repeat the splitting above with h(Xθj) in place of Wθj. Expanding
to leading order, the midpoint sum exceeds the left-endpoint sum by
j∑(h′σ)(tj,Xtj)Aj(Aj+Bj),
whose conditional mean is ∑j(h′σ)(tj,Xtj)⋅21(t/n), since E(Aj2)=t/(2n) and E(AjBj)=0. In the limit this is 21∫0t(h′σ)(s,Xs)ds. So, formally,
The table above is the case h(x)=x and X=W, so that h′≡1 and σ≡1. The correction is 21∫0tds=t/2, which is exactly the gap 21Wt2−21(Wt2−t) between the midpoint and left rows.
If the integrand does not depend on the state then h′≡0 and there is no correction at all. Additive noise is unambiguous.
In the Stratonovich equation (3) the integrand is the diffusion coefficient, h=σ. So the correction is 21∫0t(σxσ)(s,Xs)ds, and (3) describes the same process as the Itô equation (1) with
μ=μ~+21σxσ.
This derivation is formal — it discards remainders without estimating them — and is ours, not the source's; the source gives only the case tabulated above. Where it matters, verify the conversion against a case with a known answer.
Which to use. Three practical points. The Itô integral is a martingale and obeys the isometry, which is why every result in this unit is stated for it. The Stratonovich integral obeys the ordinary chain rule, which makes it the convenient one when a change of variables or a coordinate-free formulation is wanted. And when the noise is additive, σx=0 and the two coincide — so the Langevin equation below, and every linear equation (4), is unambiguous.
Method(Solving a linear SDE with additive noise)Original
To solve dXt=(α(t)Xt+ν(t))dt+δ(t)dWt with X0 given:
Check the form. The diffusion coefficient δ must not depend on X. If it does, the method does not apply; try instead a substitution such as Y=logX, remembering that Itô's formula, not the chain rule, governs the substitution.
Build the integrating factorΦt=exp(−∫0tα(r)dr), exactly as for the deterministic equation x˙=αx+ν.
Differentiate the product with Itô's formula, not with the product rule. Confirm that f(t,x)=Φtx has fxx=0, so that the Itô correction vanishes and d(ΦtXt)=Φtνdt+ΦtδdWt.
Integrate and divide to get (5).
Read off the mean by deleting the stochastic integral, whose expectation is zero.
Read off the variance with the Itô isometry: square the coefficient of dW inside the integral, integrate it against ds, and divide by Φt2. For the autocovariance at two times, integrate only up to the earlier one, as in (7).
Declare the law. With deterministic coefficients and constant X0, the answer is Gaussian, so steps 5 and 6 determine it completely; any probability is then a normal-tail evaluation.
Sanity-check the limits. As t→0 the mean should return X0 and the variance should vanish; as t→∞ a stable equation (∫tα→−∞) should give a finite limiting variance.
Example(An Ornstein-Uhlenbeck process, with numbers)Original
Take β=21, σ2=3 and X0=2 in equation (8); the time unit is whatever makes the relaxation rate 21.
Step 1: integrating factor.α≡−21, so Φt=et/2.
Step 2: solve. By (5) with ν≡0,
Xt=2e−t/2+3∫0te−(t−s)/2dWs.
Step 3: mean.
E(Xt)=2e−t/2,so at t=4:2e−2=0.270671.
Step 4: variance by the isometry.
var(Xt)=3∫0te−(t−s)ds=3(1−e−t),so at t=4:3(1−e−4)=2.945053,
using e−4=0.0183156. The standard deviation is 2.945053=1.716116.
Step 5: check against the general formula. Item 1 of the Ornstein-Uhlenbeck statement predicts variance σ2(1−e−2βt)/(2β)=3(1−e−4)/1, which agrees. The limiting variance is σ2/(2β)=3, and 2.945053 is already short of it by only 1.83%: after four time units, that is two relaxation times, the process has nearly forgotten where it started.
In the stationary state the same question would give 1−Φ(3/3)=1−Φ(1.732051)=0.041632; the process is still slightly displaced towards its starting value 2, which raises the probability.
Step 7: autocovariance. By (7), cov(Xs,Xt)=13(e−(t−s)/2−e−(t+s)/2) for s≤t. At s=4, t=6 this is 3(e−1−e−5)=1.083424, and with var(X6)=3(1−e−6)=2.992564 the correlation is
2.945053×2.9925641.083424=0.364947,
against the stationary value e−β⋅2=e−1=0.367879. The correlation time, the lag at which the stationary correlation falls to e−1, is 1/β=2 time units.
Application(The Langevin equation, equipartition and a colloidal bead)AdaptedAdapted from Grimmett & Stirzaker (2020), §13.2, pp. 579-580
Grimmett and Stirzaker introduce the Ornstein-Uhlenbeck process as a repair of the Wiener process as a model of Brownian motion. The objection is physical: a Wiener path is nowhere differentiable, so it is not a path a particle of non-zero mass can follow. Their remedy is to let the velocity, not the position, be the fluctuating variable, with the increment of velocity made of two parts, a frictional drift proportional to the current velocity and a random kick from collisions; this is what they call the Langevin equation. The position is then recovered by integrating the velocity, and so is continuously differentiable.
In continuum form, for a particle of mass m in a fluid of drag coefficient γ,
which is (8) with β=γ/m and σ=c/m. The noise is additive, so there is no Itô-Stratonovich ambiguity, and the method above applies unchanged. By item 2 of the Ornstein-Uhlenbeck statement the velocity settles to
V∞∼N(0,2βσ2)=N(0,2γmc2).
Fixing the noise strength (this part is ours, not the source's). Thermal equilibrium at temperature T requires equipartition, 21m⟨V2⟩=21kBT, so ⟨V2⟩=kBT/m. Comparing with the display above forces
c2=2γkBT.
The strength of the random kicks is tied to the strength of the friction: the same molecular collisions cause both. This is the simplest fluctuation-dissipation relation, and it is a consequence of the Ornstein-Uhlenbeck variance formula plus one thermodynamic input.
A number. Take a silica sphere of radius a=0.5μm and density 2200kgm−3 in water at T=300K, with viscosity η=8.9×10−4Pas. Then
Both numbers explain why the Wiener process is nevertheless an excellent model of what a microscope sees. The velocity decorrelates in a tenth of a microsecond; on any observable timescale the position is the integral of an already-decorrelated velocity, which by The central limit theorem is approximately a Wiener process, with diffusion constant kBT/γ. Resolving the ballistic regime below 0.14μs, and so seeing the Ornstein-Uhlenbeck velocity directly, needs an optical trap and nanosecond timing.
States the forward (Fokker-Planck) and backward Kolmogorov equations satisfied by the transition density of a one-dimensional diffusion, derives each from Ito's formula - the forward equation by testing against compactly supported functions and integrating by parts, the backward equation from the martingale property of the conditional expectation - and gives the zero-flux recipe for the stationary density, which produces the Boltzmann distribution.
Then a is the instantaneous mean, or drift, and b the instantaneous variance of D. For s≤t the transition density is
f(t,y∣s,x)=∂y∂P(D(t)≤yD(s)=x),
a function of four variables: the pair (t,y) is forward data, the pair (s,x)backward data.
A solution of the stochastic differential equation dXt=μ(t,Xt)dt+σ(t,Xt)dWt of Stochastic differential equations is such a process, with a=μ and b=σ2; conditions (2) and (3) are the infinitesimal mean and variance of the increment read off from that equation. Condition (1) is not read off in the same way, but it holds as soon as μ and σ are bounded, say by K. Split the increment as Xt+h−Xt=Ah+Mh with Ah=∫tt+hμdr and Mh=∫tt+hσdW. The drift part is bounded pathwise, ∣Ah∣≤Kh, so Ah4≤K4h4. For the Itô part, apply Itô's formula to y↦y4 along M: Mh4=4∫0hMr3σdW+6∫0hMr2σ2dr. The stochastic integral has mean zero, and the isometry gives E(Mr2)≤K2r, so E(Mh4)≤6K2∫0hK2rdr=3K4h2. Since (u+v)4≤8(u4+v4), this gives E([Xt+h−Xt]4∣Xt=x)=O(h2), and Markov's inequality applied to the fourth power gives P(∣Xt+h−Xt∣>ϵ∣Xt=x)≤O(h2)/ϵ4=o(h). Associated with (a,b) is the second-order differential operator
"Formal adjoint" means that ∫(Ltg)pdy=∫g(Lt∗p)dy for g smooth and compactly supported — integration by parts, twice, with no boundary terms. The discrete-state analogue of this pair is the generator of a continuous-time chain and its transpose; see Continuous-time Markov chains.
Theorem(The forward and backward Kolmogorov equations)AdaptedAdapted from Grimmett & Stirzaker (2020), §13.3, equations (13.3.7) and (13.3.8) with the surrounding hypotheses, pp. 583-584
Let D be a diffusion satisfying (1), (2) and (3), and suppose in addition that D has the Markov property, so that the conditional law of D(t) given the whole history up to time s depends on that history only through D(s); let f(t,y∣s,x) be its transition density. The Markov property is a hypothesis of both equations, not a by-product: it is what makes f a function of the four variables alone, it is the route by which the source reaches (6) and (7), and it is used in Step B2 of the proof below. Subject to further technical conditions, for which Grimmett and Stirzaker refer to Feller (1971), pp. 332-335, and which they do not reproduce, f satisfies the following two partial differential equations.
Forward equation (Fokker-Planck). A note on the name before the equation: Grimmett and Stirzaker call (6) simply the forward equation, and the words "Fokker-Planck" do not appear in their chapter on diffusions. The name is the one used in physics and in the partial-differential-equations literature, where Evans (Evans (2010), §1.2, items 7 and 8, p. 4) lists it as a standard linear equation alongside Kolmogorov's equation; both names are used here, for the same equation. In the forward variables, with (s,x) held fixed,
The pair (a,b) determines everything. The source emphasises that no further information about the distribution of an increment is needed: once the instantaneous mean and variance are known, the transition density is pinned down (subject to boundary conditions, which the source does not specify at this point). This is convenient in applications, where a and b are usually what the physical description supplies.
The density of the process itself obeys (6). Fix the initial state, write p(t,y)=f(t,y∣0,x0); then p solves (6) with initial condition the point mass at x0. More generally, if D(0) has density p0, then p(t,⋅)=∫f(t,⋅∣0,x)p0(x)dx solves (6) with p(0,⋅)=p0, by linearity.
Three instances, as the source gives them. With a≡0, b≡σ2 the process is a Wiener process and (6) is the heat equation ∂tf=21σ2∂y2f, whose solution is the Gaussian transition density; see The heat kernel. With a≡m, b≡σ2 it is a Wiener process with drift, D(t)=σW(t)+mt. With a(t,x)=−βx, b≡σ2 it is the Ornstein-Uhlenbeck process and the forward equation reads ∂tf=β∂y(yf)+21σ2∂y2f.
Idea of the proofOriginal
The two equations answer two different questions, and that is the whole content of the words "forward" and "backward". The forward equation propagates a distribution: given where the process is now, where will the probability be later? The backward equation propagates an observable: given a payoff to be collected at a fixed future time, what is its expected value as a function of where and when you start? The first acts on densities and uses L∗; the second acts on functions of the state and uses L. They are adjoint because total probability is conserved.
Both come out of Itô's formula, in one line each, once the right thing is differentiated. Apply Itô's formula to g(Xt) for a smooth compactly supported g and take expectations. The stochastic integral drops out, because it is a martingale with mean zero, and what remains says
dtdEg(Xt)=E[(Ltg)(Xt)].
Written out with the density this is ∫g∂tp=∫(Ltg)p. Move the derivatives off g and onto p by integrating by parts — which is legitimate precisely because g has compact support — and the identity becomes ∫g∂tp=∫gLt∗p for every test function g. Two functions with the same integral against every test function are equal, and that is the forward equation.
The backward equation is the statement that a conditional expectation of a fixed future quantity is a martingale. Freeze a terminal time t and a payoff g, and let u(s,x) be the expected payoff started from x at time s. By the tower property and the Markov property, u(s,Xs) is a martingale in s. Apply Itô's formula to it: a martingale that is also an Itô process can have no dt term, so its drift us+Lsu must vanish identically. That is (7), tested against g.
Proof(Of the forward and backward equations)AdaptedAdapted from Grimmett & Stirzaker (2020), §13.9, Theorem (13.9.2), p. 610, §13.8, Theorem (13.8.11), p. 608, and §13.12, Problem (13.12.24), p. 629
We derive (6) and (7) for a diffusion given as the solution of dXt=a(t,Xt)dt+σ(t,Xt)dWt with b=σ2. Grimmett and Stirzaker instead obtain the equations directly from the Markov property and the moment conditions (1)-(3), deferring the technical conditions to Feller; the route below buys concreteness by assuming more, and we list the assumptions rather than hide them.
Standing assumptions.a and σ are continuous and bounded, with bounded continuous first and second derivatives in x; Xt has a density p(t,⋅) for each t>0, jointly continuously differentiable once in t and twice in y, together with ap and bp; differentiation under the integral sign in t is permitted; the equation has, for every starting point (s0,x0), a solution on [s0,∞) with Xs0=x0, and that solution is a Markov process with transition density f; and the function u of Part B below is twice continuously differentiable on [0,t]×R, with bounded derivatives — twice in s as well as twice in x, because Step B3 applies the general Itô formula of The Ito integral and Ito's formula to u(v,Xv) and that formula asks for a C2 function of the pair (s,x); the weaker C1,2 regularity that the equation itself delivers would do, but this handbook states Itô's formula only in the C2 form. These are the sort of conditions Feller's list supplies; none of them is proved here.
Part A: the forward equation.
Step A1: apply Itô's formula to a test function. Let g∈Cc∞(R), that is, smooth with compact support, and set Yt=g(Xt). Since g does not depend on t, Itô's formula gives
Step A2: the stochastic integral has mean zero. Its integrand σ(r,Xr)g′(Xr) is adapted and bounded, because σ is bounded by assumption and g′ is continuous with compact support; so ∫0tE[(σg′)2]dr<∞ and the integrand is admissible. By the martingale property of the Itô integral, the last term has expectation 0 for every t.
Step A3: take expectations and differentiate in t. Take expectations in Step A1 and interchange E with ∫0tdr. The interchange is Fubini's theorem, licensed because the integrand (Lrg)(Xr) is bounded — a and b are bounded by the standing assumptions and g′, g′′ are continuous with compact support — and the region [0,t]×Ω has finite measure. Writing Lr for the generator (4),
Eg(Xt)=Eg(X0)+∫0tE[(Lrg)(Xr)]dr.
The integrand is bounded and continuous in r, so the right-hand side is differentiable in t and
dtdEg(Xt)=E[(Ltg)(Xt)].(8)
Step A4: write (8) with the density. By the definition of p, equation (8) says
differentiation under the integral sign on the left being permitted by assumption.
Step A5: move the derivatives onto p. Integrate the first term on the right by parts once and the second twice. The support of g is compact, so every boundary term vanishes, and we get
Step A6: remove the test function. The difference h=∂tp−Lt∗p is continuous and satisfies ∫gh=0 for every g∈Cc∞. If h(y0)=0 for some y0, then by continuity h keeps its sign on an interval around y0, and choosing a non-negative g supported in that interval and positive at y0 makes ∫gh=0, a contradiction. Hence h≡0, which is (6). Applying the argument to the conditional law given Xs=x gives the same equation for f(t,y∣s,x).
Part B: the backward equation.
Step B1: the conditional expectation of a terminal payoff. Fix t and g∈Cc∞(R) and put, for 0≤s≤t and everyx∈R,
u(s,x)=E(g(Xt)∣Xs=x)=∫g(y)f(t,y∣s,x)dy.
This defines u at every point of [0,t]×R, because the transition density is given for every starting point, not only for the points one particular trajectory visits.
Step B2: fix a starting point and run the diffusion from it. Fix (s0,x0) with 0≤s0<t, and let X be the solution on [s0,t] with the deterministic initial value Xs0=x0, which the standing assumptions provide. Put Mv=u(v,Xv) for v∈[s0,t]. By the Markov property, E(g(Xt)∣Fv)=u(v,Xv), so for r≤v the tower property gives
Write αv for the bracket; by the standing assumptions it is bounded, and it is continuous in v along almost every path.
Step B4: the drift vanishes at the starting point. The stochastic integral in Step B3 has a bounded adapted integrand, so it is a martingale and has mean zero. Taking expectations in Step B3 between s0 and s, and using EMs=EMs0=u(s0,x0),
E(s−s01∫s0sαvdv)=0,s0<s≤t.
Now let s↓s0. Along almost every path the average inside tends to αs0, because v↦αv is continuous at v=s0; and Xs0=x0 is deterministic, so that limit is the number(us+aux+21buxx)(s0,x0). The averages are bounded uniformly, so bounded convergence lets the limit pass through the expectation. (The source states this zero-drift criterion for the special case of the heat equation in its Problem (13.12.24).)
The point (s0,x0) was fixed before the process was started and was arbitrary in [0,t)×R — this is why Step B2 starts the diffusion afresh from each point instead of following one trajectory, which would deliver the identity only where the law of Xs0 lives. So the identity holds at every (s,x) with s<t, and at s=t by continuity in s.
Step B6: remove the test function. Substituting u(s,x)=∫g(y)f(t,y∣s,x)dy and differentiating under the integral sign,
∫g(y)[∂s∂f+a(s,x)∂x∂f+21b(s,x)∂x2∂2f]dy=0
for every g∈Cc∞. By the argument of Step A6 the bracket vanishes, which is (7). □
Remark(Where the hypotheses are used, and two conventions)Original
Compact support of the test function is used three times, in Steps A2, A3 and A5. In A2 and A3 it is one half of the boundedness of g and its first two derivatives, which is what makes the integrand of the stochastic integral admissible and what licenses the Fubini interchange. In A5 it is what kills the boundary terms of the integration by parts, and that is the use that shapes the statement. On a bounded interval with reflecting or absorbing barriers the boundary terms do not vanish, and (6) must be supplemented by boundary conditions; Grimmett and Stirzaker note at this point that they avoid specifying boundary conditions. Zero flux at the ends is the reflecting condition, f=0 at the ends the absorbing one.
Boundedness of σ and the compact support of g′ are used in Step A2 to put the integrand in the admissible class, so that the stochastic integral is a martingale and not merely a local martingale. Without this the mean of the last term in Step A1 need not be zero, and (8) fails.
The continuity of the drift α at the starting time is used in Step B4. The zero-drift conclusion is read off from a limit of averages over [s0,s] as s↓s0, so v↦αv must be right-continuous at v=s0; it is, because u has continuous derivatives, a and b are continuous, and X has continuous paths. For merely measurable coefficients one gets the identity only for almost every s, which is enough for (7) to hold almost everywhere but not pointwise.
Starting the process afresh from each fixed point is what lets Step B5 reach every x. Had we followed one trajectory, Step B4 would give the identity at the random points (r,Xr), that is, on the support of the law of Xr, and continuity would extend it only to the closure of that support. Nothing here makes that support the whole line — the next paragraph points out that b may vanish, and then the law of Xr can be confined to part of the line. Fixing (s0,x0) in advance avoids the question.
Non-degeneracy of b is not needed for the derivation but is needed for the equations to behave like diffusion equations. In the language of The heat kernel and of Evans (Evans (2010), §7.1, equations (2)-(4), pp. 349-350), the operator ∂t−L is uniformly parabolic exactly when 21b(t,x)≥θ>0 for some constant θ. Where b vanishes the equation degenerates: for Feller's branching diffusion, a(t,x)=αx and b(t,x)=βx (Example (13.3.12) of the source, which calls the two constants a and b; those letters are already spoken for here, as the instantaneous mean and variance, so we write α and β for them), b vanishes at the origin, which is exactly the absorbing state of the branching process, and no smoothing occurs there.
Two sign conventions, and the dictionary between them. The partial-differential-equations literature writes the same pair of equations with different letters and a different sign. Evans (Evans (2010), §1.2, items 7 and 8, p. 4) lists, among the standard linear equations, Kolmogorov's equation
ut−i,j∑aijuxixj+i∑biuxi=0
and the Fokker-Planck equation
ut−i,j∑(aiju)xixj−i∑(biu)xi=0.
In one dimension these become ut=auxx−bux and ut=(au)xx+(bu)x. Comparing the second with (6) gives the dictionary
aEvans=21bG&S,bEvans=−aG&S,
so that Evans's first-order coefficient is the negative of the probabilists' drift, and his second-order coefficient is half the instantaneous variance. With that substitution his Kolmogorov equation is (7) run forward in time, τ=−s. The two sources agree; only the normalisation differs. Everything in this unit follows the probabilists' convention of (6) and (7). The n-dimensional forms are Evans's as displayed, with aij=21σσT built from the diffusion matrix.
Method(Finding a stationary density)Original
Suppose a and b do not depend on t and b>0. A stationary density is a probability density p∞ with L∗p∞=0; a process started from it has the same law at every time.
Write the forward equation as a conservation law. Define the probability flux
J(y)=a(y)p(y)−21dyd[b(y)p(y)],
so that (6) reads ∂tp=−∂yJ. A stationary density therefore has J constant in y.
Set the constant to zero when the state space is the whole line, p is to be integrable, and
a(y)p(y)→0andb(y)p(y)→0(y→±∞).
Suppose J≡c. Then 21(bp)′=ap−J=ap−c→−c, so a non-zero c would force bp to grow linearly, contradicting bp→0. Hence c=0. (On a finite interval with a current driven through it, keep the constant; that is the non-equilibrium case.)
The decay hypothesis, not boundedness of a and b, is what the argument needs. When a and bare bounded it is automatic, since an integrable p with integrable derivative has limits at ±∞ and those limits must be 0. When they are not — the Ornstein-Uhlenbeck drift a(y)=−βy of step 5 is unbounded — it has to be checked on the candidate density produced by step 3, and there it holds, because that density decays like a Gaussian and so kills any polynomial factor.
Solve the first-order equation J≡0. Put q=bp. Then 21q′=(a/b)q, a separable linear equation, so
Normalise. Fix C by ∫p∞=1. If the integral diverges there is no stationary density, and that is the correct answer: a Wiener process with or without drift has none.
Check against a known case. For the Ornstein-Uhlenbeck process, a=−βy and b=σ2, the exponent is −βy2/σ2 and (9) gives a centred normal density of variance σ2/(2β), matching the limit law quoted in the statement above.
Recognise the Boltzmann form. Whenever b≡2κ is constant and a=−V′ is a gradient, (9) collapses to
p∞(y)=Z1e−V(y)/κ,
with Z the normalising constant. The letter κ is used here, not D, because D already names the process. The next block identifies κ with kBT divided by the friction, and the exponent with the Boltzmann factor.
Example(A Wiener process with drift, from the equation)Original
Take a≡m and b≡σ2, both constant, with D(0)=0. We solve the forward equation rather than quoting the process.
Step 1: the equation. By (6), with a and b constant the derivatives pass through:
∂t∂p=−m∂y∂p+2σ2∂y2∂2p,p(0,⋅)=δ0.
Step 2: remove the drift by moving with it. Put z=y−mt and p~(t,z)=p(t,z+mt). By the chain rule ∂tp~=∂tp+m∂yp, so the first-order term cancels and
∂t∂p~=2σ2∂z2∂2p~,
the heat equation with diffusivity κ=σ2/2.
Step 3: quote the heat kernel. By The heat kernel, the solution from a point mass at the origin is p~(t,z)=(4πκt)−1/2e−z2/(4κt). With κ=σ2/2 this is (2πσ2t)−1/2exp(−z2/(2σ2t)).
Step 4: transform back.
p(t,y)=2πσ2t1exp(−2σ2t(y−mt)2),
that is, D(t)∼N(mt,σ2t) — which is the law of σW(t)+mt, as the source records in its Example (13.3.10). Substituting this p back into the equation of Step 1 confirms it. Writing u=y−mt, one finds ∂yp=−pu/(σ2t) and ∂y2p=p[u2/(σ4t2)−1/(σ2t)], while differentiating in t (remembering that u depends on t) gives the same thing: both sides equal
p[σ2tmu+2σ2t2u2−2t1].
Step 5: no stationary density. Run the method: a/b=m/σ2 is constant, so (9) gives p∞∝e2my/σ2, which is not integrable for any m, including m=0. There is no stationary density, as there should not be: the process wanders off.
Step 6: a number. Take m=0.8 and σ2=0.25 in their respective units. At t=4 the mean is mt=3.2 and the variance is σ2t=1.0, so
P(D(4)>5)=1−Φ(1.05−3.2)=1−Φ(1.8)=0.035930.Application(Boltzmann statistics and a bistable colloid)Original
An overdamped particle at position X in a potential U, in a fluid of friction coefficient γ at temperature T, obeys the Smoluchowski equation: the drift is the force divided by the friction and the diffusion constant is κ=kBT/γ, so
dXt=−γU′(Xt)dt+2κdWt,a(y)=−γU′(y),b(y)=2κ.
The forward equation is ∂tp=∂y[(U′/γ)p+κ∂yp], whose bracket is minus the flux: the first term is the drift down the potential, the second is Fick diffusion. Step 6 of the method applies with V=U/γ and gives
p∞(y)=Z1exp(−γκU(y))=Z1exp(−kBTU(y)),
the Boltzmann distribution, once κ=kBT/γ is inserted. This is worth pausing on: equilibrium statistical mechanics is here a one-line consequence of the stationary Fokker-Planck equation plus the Einstein relation between diffusion and mobility. It is also a consistency check on that relation — any other choice of κ would give the wrong equilibrium.
A bistable trap, with numbers. A colloidal bead in a double optical trap sees, in reduced units, the potential
U(y)=y4−2y2,
with minima at y=±1 of depth U(±1)=−1 and a barrier at y=0 of height 1. Take γκ=kBT=0.25, that is a barrier of four times the thermal energy. Then p∞(y)=Z−1e−4(y4−2y2) with, by numerical quadrature, Z=51.7502.
Barrier suppression.p∞(1)p∞(0)=e−[U(0)−U(1)]/(γκ)=e−4=0.018316: the bead is found at the barrier about one fiftieth as often as at a minimum.
Occupation of the barrier region.P(∣X∣<0.5)=Z−1∫−0.50.5e−4(y4−2y2)dy=0.041655, so the bead spends about 4% of its time within half a unit of the barrier top, and the rest in the two wells.
Mean square displacement.E(X2)=0.917671 in the stationary state, close to the 1 it would be if the bead sat exactly in one well or the other.
Symmetry.P(X>0)=21 exactly, by the evenness of U.
The last point is the one to be careful with in practice. The stationary density is symmetric, but a single experimental run of finite length is not: the bead sits in one well for a long time, hops, and sits in the other. The time between hops grows like eΔU/(kBT) — Kramers' law, which the stationary density alone does not give, since it is a statement about first passage rather than about the equilibrium distribution. With the barrier only four thermal units high the hops are frequent enough to see; at ten thermal units they are not, and the measured histogram over any feasible run would be one well, not two.