17. Stochastic Processes¶
PDF pages 697–832
Chapter 17 Stochastic Processes Before tackling quantum measurements head-on, we will first examine some of the basic mathematics for handling measurements, and in particular continuous quantum measurements. We will thus need to look at some basics in the area of stochastic processes—that is, the mathematics for modeling systems as having underlying randomness influencing the dynamics.1 Historically, the term ‘‘stochastic’’ has also been used to refer to low-dimensional, chaotic dynamics of Hamiltonian systems, which is not what we mean here. By stochastic we are referring to a truly random element that is not predictable even in principle. This is sensible for modeling quantum measurements, which are considered to be sources of true randomness. However, despite the inherent unpredictability, we can fruitfully model stochastic systems by building on the basic formalism introduced here. 17.1 Finite Random Walks, Diffusion, and the Central Limit Theorem One central problem in statistical mechanics that is useful in quantum optics—and indeed underlies much of the formalism of quantum measurement that we will develop—is the random-walk process. Suppose a random walker takes a random step of size X with probability density f(x) between periodic intervals of duration ∆t. Let’s assume that all the steps are statistically independent, and the probability distribution is characterized by
(17.1) After N steps (N large), where has the walker ended up? The central limit theorem says that the probability density of the accumulated displacement SN := N X j=1 Xj (17.2) for N steps is Gaussian with zero mean and variance N\sigma2. That is, the width (standard deviation) is \sigma \sqrt N.2 The probability distribution thus becomes asymptotically Gaussian with a time-dependent width of
r t ∆t. (17.3) This random-walk behavior is characteristic of a diffusion process, which is a transport process by which the distribution grows as t1/2, ∆x ∼D t1/2, (17.4) 1Note that we will be giving just an introductory overview, and will sacrifice rigor in favor of intuition; a good rigorous introduction is W. Horsthemke and R. Lefever, Noise-Induced Transitions: Theory and Applications in Physics, Chemistry, and Biology (Springer, 1984). Another good introduction is the classic C. W. Gardiner, Handbook of Stochastic Methods 1st ed. (Springer, 1983). 2Recall that the variance of X is defined by Var[X] := D
=
X2 −\langle X\rangle 2, and the standard deviation is the square root of the variance.
17.1.1 Two-Step Distribution¶
Chapter 17. Stochastic Processes where for the random walker the diffusion coefficient is D = \sigma/ \sqrt ∆t. Note that within certain restrictions, the final distribution is Gaussian, independent of the one-step distribution. 17.1.1 Two-Step Distribution Before proving the full central limit theorem, we will examine the probability density after exactly two steps. The mathematical problem is as follows: let X1 and X2 be independent random variables with probability density functions f1(x) and f2(x), respectively. That is, the probability that X1,2 is between x and x + dx is f1,2(x) dx. Then we can ask, what is the probability density of X1 + X2? To answer this, we can note that X1 + X2 = x for any pair of values of X1 and X2 that happen to add up to x. But then we must sum over all such pairs. The probability that both X1 and X2 will both have particular probabilities is the product of the individual probabilities since the variables are independent. Thus, expressing what we said in equation form,
X x′,x′′
(17.5) We can translate this statement in terms of the probability densities and implement the constraint as a \delta-function (with a factor of dx, so that the \delta-function registers unity when the condition is met). Letting f+(x) denote the probability density of X1 + X2,
Z \infty −\infty dx′ Z \infty −\infty dx′′ f1(x′) f2(x′′) \delta(x′ + x′′ −x) dx. (17.6) Evaluating the x′′ integral, we see that the probability density of the sum is the convolution of the individual densities,
Z \infty −\infty
(17.7) where we use the ∗symbol to denote convolution of two functions. Note that this result is general in that it doesn’t assume any particular form for f1(x) or f2(x).
the probability density for two steps is
(17.8) (two-step probability density) i.e., the convolution of the one-step distribution with itself. Recall that the convolution ‘‘smears’’ one function with another, and so as the effect of the second step is to smooth the one-step distribution. The idea behind the central limit theorem is that this smoothing continues until the distribution is Gaussian after many steps. 17.1.1.1 Example 1: Convolution with a Delta Function As an example of the general idea of the convolution of two functions f and g,
Z \infty −\infty dx′ f(x′) g(x −x′), (17.9)
then
Z \infty −\infty
(17.10) The effect of convolution with a delta function is thus simply to do nothing: convolution with a delta function is just the identity operation.
17.1.2 Convolution Theorem¶
17.1 Finite Random Walks, Diffusion, and the Central Limit Theorem In terms of the random walk, \delta(x) as a one-step probability function simply corresponds to a step of zero length, or just taking no step at all. Thus, it makes intuitive sense that the distribution isn’t changed by convolution with \delta(x). In general, when g(x) is some other function, the convolution ‘‘smears’’ f(x) with the convolution kernel g(x). Typically, we will use centered kernels; the effect of a displaced kernel is simply to displace the convolution by the same amount. For example, if
(17.11) then
Z \infty −\infty
(17.12) which is just the displaced version of the original. 17.1.1.2 Example 2: Convolution of Box Functions As a slightly more complicated example, consider the convolution of box functions, both given by
1, |x| \le 1/2 elsewhere, (17.13) which here are properly normalized to correspond to probability distributions. The convolution consists of displacing g(x′) by x, multiplying the functions together, and integrating. For this simple case (box functions of unit height), the convolution (product) just turns out to be the area where the two functions overlap. xoo' xo'o=o0 1/2 -1/2 x fo(xoo') go(xo-oxoo') (foog)o(x) When the displacement is large, |x| > 1, the boxes don’t overlap at all, so the convolution is zero. Otherwise, the overlap area varies linearly with the displacement, so the convolution is a triangle function. x xo=o0 1/2 -1/2 -1 fo(x) (foofo)o(x)
or ‘‘blurring’’ effect of the convolution. The original functions were discontinuous, but the convolution is continuous. The convolution is also wider than the original functions. As we will see, continued, successive convolutions will make the distribution look Gaussian. 17.1.2 Convolution Theorem Now that we brought up the convolution, we may as well discuss how to compute it. The convolution theorem gives an easy way to evaluate the convolution integral in Eq. (17.7), both in an intuitive and a
Chapter 17. Stochastic Processes computational sense. The convolution theorem states that the Fourier transform of the convolution is the product of the Fourier transforms of the individual functions: F[f ∗g] = F[f]F[g]. (17.14) (convolution theorem) To prove this, we’ll just compute the explicit form of F[f ∗g]. This will be very much a physicist’s proof, not a mathematician’s proof, in that we’ll just assume the functions are nice enough that all the integrals simply exist. First of all, in our notation here, the Fourier and inverse transforms have the form
2\pi Z \infty −\infty dk ˜f(k) eikx,
Z \infty −\infty dx f(x) e−ikx, (17.15) where ˜f(k) ≡F[f(x)]. It’s important to make this explicit, since the result depends on the normalization convention we choose for the Fourier transform. Then computing the Fourier transform of f ∗g, F[f ∗g] = F Z \infty −\infty dx′ f(x′) g(x −x′) = Z \infty −\infty dx Z \infty −\infty dx′ f(x′) g(x −x′) e−ikx = Z \infty −\infty dx Z \infty −\infty dx′ f(x′) e−ikx′g(x −x′) e−ik(x−x′). (17.16) Letting x −\rightarrow x + x′, F[f ∗g] = Z \infty −\infty dx Z \infty −\infty dx′f(x′) e−ikx′g(x) e−ikx = Z \infty −\infty dx′ f(x′) e−ikx′ Z \infty −\infty dx g(x) e−ikx = F[f]F[g]. (17.17) Thus, to convolve two functions, just follow this recipe: Fourier transform both functions, multiply them together, then compute the inverse Fourier transform. Mathematically, we can write f ∗g = F −1 {F[f]F[g]} . (17.18) Since Fourier transforms of common function are usually already known, the convolution theorem provides a shortcut for evaluating the full convolution integral. 17.1.2.1 Example: Convolution of Two Gaussians Since it’s easy to compute the Fourier transform of Gaussian distributions, let’s use the convolution theorem to convolve two Gaussians. Let’s write the two functions as
(17.19) The Fourier transform of a Gaussian is also a Gaussian, and in particular
(17.20) Then the product of the Fourier transforms is
(17.21)
17.1.3 Proof of the Central Limit Theorem¶
17.1 Finite Random Walks, Diffusion, and the Central Limit Theorem Finally, we invert the Fourier transform to obtain the convolution:
p
− x2
. (17.22) Recall that the standard (normalized) form of the Gaussian is \sqrt
−(x −µ)2 2\sigma2 , (17.23) where the µ is the mean and \sigma is the standard deviation (\sigma2 is the variance). The standard deviation is a common measure of the width of a Gaussian function. Note that f(x) has standard deviation \alpha/ \sqrt 2, g(x) has standard deviation \beta/ \sqrt 2, and (f ∗g)(x) has standard deviation p (\alpha2 + \beta2)/2, so that the standard deviations add in quadrature as a result of the convolution. Thus, the convolution of Gaussians is still Gaussian, but the blurring effect of the convolution makes the convolved Gaussian wider than the original functions. 17.1.3 Proof of the Central Limit Theorem Now we extend the two-step analysis above analysis to N steps. Let X1, . . . , XN be independent, identically distributed random variables. Let f(x) be the probability density function of each of the Xj. Defining the sum by SN := N X j=1 Xj, (17.24) we will now ask, what is the probability density fSN (x) of SN? Evidently, we can iterate Eq. (17.8) to obtain
(17.25) where the result is the successive convolution of N copies of f (for N −1 total convolution operations). However, it turns out that this distribution becomes simple for large enough N. The central limit theorem states that, provided that the mean and variance of the Xj exist, with
for large N with
(17.26) (central limit theorem) (The mean and variance are in fact exact, whereas the form of the distribution is valid for large N.) This is a rough statement, since ‘‘becomes asymptotically Gaussian’’ is an imprecise statement. So let’s clean this up a bit. The central limit theorem states that the probability density function fZN (x) of the centered, scaled statistic ZN := SN −Nµ \sigma \sqrt N (17.27) converges to the ‘‘standard normal’’ (Gaussian) distribution fZN (x) −\rightarrow \sqrt
(17.28) which is the special Gaussian with mean 0 and unit variance. Let’s prove this now3. To evaluate the convolutions in Eq. (17.25), we need to employ the convolution 3This is the physicist’s proof; the rigorous version is in T. W. Körner, Fourier Analysis (Cambridge, 1988), starting on p. 349.
Chapter 17. Stochastic Processes theorem. Taking the Fourier transform of f(x),
Z \infty −\infty dx f(x) e−ikx = \infty X j=0 Z \infty −\infty dx f(x)(−ikx)j j!
- O(k3). (17.29) Here, we Taylor-expanded e−ikx and then used the fact that the terms of the expansion were proportional to expectation values \langle Xj\rangle . In particular, note that in probability theory the characteristic function of a probability density, given by
(17.30) is an important tool for manipulating probabilities. This is more cumbersome than necessary, so let’s recompute the expansion in Eq. (17.29) for the centered, scaled variable Zj = Xj −µ \sigma \sqrt N , (17.31) with corresponding probability density fZ(x). The centering effectively zeroes the mean, and the rescaling changes the factor in front of the variance, with the result
2N + O " k \sqrt N 3# . (17.32) The convolution theorem says that to calculate the transform of the N-fold convolution, we just compute ˜fZ(k) to the Nth power:
h ˜fZ(k) iN =
1 −k2 2N + O " k \sqrt N 3#!N . (17.33) As N becomes large, we can neglect the higher order terms beyond the first, and then use the formula lim n−\rightarrow \infty 1 + x n n = ex (17.34) to see that for large N, the transform becomes
−k2 . (17.35) But now the inverse Fourier transform of exp(−k2/2) is exp(−x2/2)/ \sqrt 2\pi, so fZN converges to a standard normal distribution as N −\rightarrow \infty. 17.1.3.1 Example: Square Distribution As a simple example of the central limit theorem, let’s try out the unit box function as the one-step distri- bution, as we tried out in Section 17.1.1.2:
1, |x| \le 1/2 elsewhere. (17.36) First note that the this function is normalized, so it represents a proper probability distribution. Thus, so do all of its self-convolutions. Let f ∗N(x) denote the convolution of f(x) with itself N −1 times. This is the
17.1 Finite Random Walks, Diffusion, and the Central Limit Theorem same as fSN (x) for the random-walk interpretation of this distribution. The central limit theorem says that asymptotically, the self-convolution becomes Gaussian,
\sqrt
N , (17.37) with zero mean, since f(x) is centered. The variance of f(x) is Z \infty −\infty
Z 1/2 −1/2 x2 dx = 1 12, (17.38) so that the width of the asymptotic Gaussian is
r N 12. (17.39) Here, f(x) is plotted with several self-convolutions f ∗N(x), along with the asymptotic form, the Gaussian of width \sigmaN = p N/12. Noáo" No=o1 No=o2 No=o3 No=o4 No=o10
-1.5 1.5 No1/o2(probability density) 1.5 0.5 As N increases, the widths of the distributions increase and the peak values decrease, but we have rescaled the axes by appropriate factors of \sqrt N to keep the distributions comparable at each step. The box function is very different from the asymptotic Gaussian. However, even the first self-convolution (a triangle function) is already pretty close to the Gaussian, and the successive self-convolutions converge fairly rapidly to the asymptotic form. 17.1.3.2 Application: Standard Deviation of the Mean Returning again to error analysis, suppose we make independent measurements X1, . . . , XN of some quantity in the laboratory. The sample mean is µN := 1 N N X j=1 Xj. (17.40) We can rewrite this as µN = SN
\sqrt N , (17.41) where the first term represents the true mean, and the second is the experimental error (statistical fluctuation in the sample mean). Applying the central limit theorem, ZN is approximately standard normal for large N, so µN is Gaussian with mean µ and standard deviation \sigma/ \sqrt N, where \sigma is the standard deviation of a single measurement. Thus, the standard deviation of the mean (also called the standard error) is \sigma/ \sqrt N. This is why, by making many measurements, it is possible to increase the accuracy of a measured quantity.
17.1.5 A Walk on the Cauchy Side¶
Chapter 17. Stochastic Processes 17.1.4 Variances Add in Quadrature In discussing random walks so far, we have been discussing the asymptotic, N-step probability distribution, and its scaling with time. However, we can make a simpler statement that does not explicitly refer to the distribution. Let X1, . . . , XN be independent random variables, but now we won’t even require them to be identically distributed. For the moment, let’s also assume\langle Xn\rangle = 0. Now consider the sum X1 +X2. Clearly the mean vanishes, and thus the variance becomes
(X1 + X2)2 =
X 2
+
X 2
=
X 2
+
X 2
= Var[X1] + Var[X2], (17.42) where we used the fact that X1 and X2 are independent, and thus their correlation function\langle X1X2\rangle factorizes into\langle X1\rangle \langle X2\rangle [the joint probability density f(x1, x2) for independent processes must have the factored form f1(x1)f2(x2)]. Thus, the variances of independent random variables, and regarding the variance as the square of the ‘‘width’’ of the corresponding probability distributions, we see that the widths add in quadrature when we add together the random variables. By subtracting \langle X1 + X2\rangle 2 =\langle X1\rangle 2 +\langle X2\rangle 2 + 2\langle X1\rangle \langle X2\rangle from each intermediate expression, it isn’t hard to see that the same result holds when \langle Xn\rangle ̸ = 0. Iterating this process, we see that the variance of the sum defined as before, SN := N X j=1 Xj, (17.43) is simply Var[SN] = N X j=1 Var[Xj]. (17.44) Again, if we take each Xn to be identical as for the random walk, and we take the variance as the square of
(17.45) (variance of the sum) or
p Var[SN] = \sqrt N \sigma. (17.46) (standard deviation of the sum) This is the same as one of the results of the central limit theorem (17.26), but this is not an asymptotic statement, it is exact. Thus, we expect the width of the sum to be precisely \sqrt N \sigma. Nevertheless, we often expect this scaling of the distribution width to hold only asymptotically, since in general the ensemble of walkers will have an initial distribution that does not match the one-step distribution (or any N-step distribution), and thus we also need to include the convolution with this initial state. 17.1.5 A Walk on the Cauchy Side Consider independent, identically distributed random variables X1, . . . , XN with Cauchy (Lorentzian) prob- ability density functions
(17.47) (Cauchy distribution) The Fourier transform is given by
(17.48)
17.1 Finite Random Walks, Diffusion, and the Central Limit Theorem as we can see by computing the inverse Fourier transform of ˜f(k):
2\pi Z \infty −\infty dk e−|k|eikx = 1 2\pi Z \infty dk e−k(1−ix) + Z \infty dk e−k(1+ix) = 1 2\pi Z \infty dk e−k(1−ix) + c.c. =
= 1 + ix
=
(17.49)
Now let’s compute the probability density of the mean µN := 1 N N X j=1 Xj. (17.50) The probability density function of the sum SN := N X j=1 Xj (17.51) is (f ∗f ∗. . . ∗f)(x) (N copies of f or N −1 convolutions), so using the convolution theorem,
h ˜f(k) iN = h e−|k|iN
(17.52) In general, if f(x) and ˜f(k) are a Fourier transform pair, then so are \alphaf(\alphax) and ˜f(k/\alpha). Thus, the inverse transform of ˜f(Nk) is f(x/N)/N. The variable µN is the same as SN except for a scaling factor of 1/N, so the probability density must be the same, but N times wider. So to get the probability density of µN, we
SN, which gives f(x) dx. Thus, f(x) is also the probability density of µN. This is different than what we expect from the central limit theorem: there, we expect the mean to have a width that is smaller than that of the one-step distribution by a factor of \sqrt N. Stated otherwise,
of the one-step distribution, which says that the widths add for the Cauchy random walk. But the central limit theorem said that variances add, or the widths should add in quadrature. Is there a contradiction here? Obviously there should be some simple resolution. The problem is that variance of Xj does not exist for a Cauchy distribution. This is because the Cauchy distribution only falls off as 1/x2 for large |x|, and so the variance integral Z \infty −\infty dx f(x) x2 (17.53) diverges. The central limit theorem implicitly assumes that the variance exists; thus, the central limit theorem does not apply to this case. This is one case of anomalous diffusion, where the diffusion coefficient diverges, because the width of the N-step distribution does not scale diffusively (i.e, it scales as t rather than \sqrt t).
17.1.6 Arbitrary Combinations of Random Variables¶
Chapter 17. Stochastic Processes One important application of such ‘‘fat-tailed’’ distributions is in the area of financial modeling,4 where Gaussian random walks do not adequately model the large jumps observed, for example, in the histories of stock prices. 17.1.6 Arbitrary Combinations of Random Variables In considering random walks, we have been considering the probability distribution corresponding to the sum (17.2) of independent random variables, in terms of the probability distributions of the separate variables—it turned out to be just the convolution of the individual distributions. But how far can we push this? Here we will develop some concepts seemingly unrelated to probability, and use them to deduce the probability density for an arbitrary (scalar) function of a set of independent random variables. 17.1.6.1 Divergence Theorem The divergence theorem is fundamental in the study of electrostatics, and the standard form states that for a vector field A, Z V
I S A \cdot ˆn dS, (17.54) where V is the volume of integration, S is the surface of the volume, ˆn is the (outward-pointing) normal vector to the surface, and dS is the surface-area element for integration over the surface of S. Let’s briefly derive this in d dimensions. Consider a box of infinitesimal volume, given by
(17.55) The flux of A(x) through the surface of this volume is given by summing over the fluxes of the two sides bounding each dimension; the ‘‘area’’ of the jth side is dx1 \cdot \cdot \cdot dxj−1dxj+1 \cdot \cdot \cdot dxd = dV /dxj, and thus we have a total flux
X j A(x + dxj ˆxj) \cdot ˆxj dV dxj −A(x) \cdot ˆxj dV dxj = X j \partial A(x) \partial xj \cdot ˆxj dV
(17.56) where the divergence is interpreted in d dimensions. Now we integrate over the volume, adding up the volume and surface contributions due to all the infinitesimal elements in the integration volume. Whenever two elements contact, their fluxes on their common surface cancel, so the only contribution in the surface integral is due to the flux on the outer surface of the integration volume, and thus we have Z V
I S A \cdot ˆn dS (17.57) (divergence theorem, n dimensions) as the generalized divergence theorem, converting a volume integral in n dimensions to an integral over the bounding hypersurface, a manifold in n −1 dimensions. 17.1.6.2 Transformation of Surface Delta Functions This divergence theorem is useful in establishing a chain-rule formula for a delta function. Recall that given a function f(x), the delta function obeys the chain rule
X x\in f −1(a) \delta(x −a) |f ′(a)| , (17.58) 4Rama Cont and Peter Tankov, Financial Modelling with Jump Processes (Chapman & Hall/CRC, 2004), Chapter 1.
17.1 Finite Random Walks, Diffusion, and the Central Limit Theorem¶
x. For a coordinate change between coordinates x and y (both in d dimensions), this formula generalizes to
\deltad(x)
(17.59) so that the scaling factor is now the Jacobian determinant of the coordinate transformation. But what if
determinant turns out to be the Euclidean norm of the vector of partial derivatives of h, as we will now show. We will start by considering arbitrary scalar functions h(x) and A(x), respectively, on Rd. Now consider the step function \Theta(h), which defines a volume (or possibly multiple, unconnected volumes). Then consider the integral Z
(17.60) which follows from changing to a surface integral via the divergence theorem, and assuming that either the surface at infinity is outside the volume defined by h (any physically relevant volume should be finite), or that A vanishes when necessary at infinity. Then using
(17.61) we have Z
Z
= − Z
(17.62) Now put
(17.63) and consider Z
Z ddx \Theta[h(x)] \nabla \cdot A(x) = − Z V
= − I S A \cdot ˆn dS = I h−1(0)
|\nabla h| dS = I h−1(0) f(x) |\nabla h| dS, (17.64) where we have identified V as the volume defined by \Theta(h); S as the surface, defined by the locus of points where h vanishes; and the normal vector ˆn = −\nabla h/|\nabla h|, since the gradient is normal to the surface, but points towards the interior of the volume (where h is increasing away from the surface). To summarize, we have identified Z
I h−1(0) f(x) |\nabla h| dS, (\delta-function constraint in integration) (17.65) as the effect of a \delta-function constraint in a volume integral, where |\nabla h| is the Euclidean norm of the gradient vector:
sX j (\partial xjh)2. (17.66)
Chapter 17. Stochastic Processes Again, the constraint here changes the volume integral to an integral over the bounding hypersurface. In this derivation, we have assumed that any function f(x) can be represented as A \cdot \nabla h, but since A is arbitrary,
any surface manifold we like). 17.1.6.3 Direct Derivation
Consider a neighborhood of a point x0 on the surface, in which one of the partial derivatives is nonvanishing, say \partial h/\partial x1̸ = 0. Then we can consider the action of the \delta function along this coordinate, as Z
Z dx1 \delta(x1 −x01) f(x)
\partial h \partial x1
x=x0
\partial h \partial x1
x=x0 , (17.67)
(17.58). Then using
v u u t d X j=1 \partial h \partial xj 2 =
\partial h \partial x1
v u u t1 + d X j=2 \partial x1 \partial xj 2 , (17.68) and then identifying dS = v u u t1 + d X j=2 \partial x1 \partial xj 2
v u u t d X j=1 \partial x1 \partial xj 2
(17.69) as the local surface element (this bears more explanation), and finally integrating over the remaining coor-
Now to return to the business of identifying the surface element (17.69). We can define the surface element in terms of the volume element as dV = dS dℓ, (17.70) where dℓis the line element normal to the surface. Then dS = dV
dℓ . (17.71) Now suppose that we express dℓin terms of x1:
(17.72) where ˆn is normal to the surface, so that dℓrepresents only the component of dx1 normal to the surface. Then dividing the volume by dℓ, we only remove the dimension normal to the surface. Using ˆn = \nabla h/|\nabla h|, we have dℓ= |\nabla h| \partial h \partial x1 dx1 = dx1 v u u t1 + d X j=2 \partial x1 \partial xj 2 , (17.73) after using Eq. (17.68), and then putting this expression into Eq. (17.71), we obtain the surface element (17.69). 5This derivation and the connection to the divergence theorem in the previous section are adapted from (the much more terse treatment in) Lars Hörmander, The Analysis of Linear Partial Differential Operators I: Distribution Theory and Fourier Analysis (Springer-Verlag, 1983), p. 136.
17.1 Finite Random Walks, Diffusion, and the Central Limit Theorem 17.1.6.4 Chain Rule for Coordinate Transformations Using Eq. (17.65), we can also derive a chain rule for integrals of the form of the left-hand side, when we compare equivalent constraints specified by different functions. First, by rescaling the function f(x), we can write Z
I h−1(0) f(x) dS(x). (17.74) Now consider the constraint h(x) and the alternate but equivalent constraint k(x):
⇐\Rightarrow
(17.75) Then the translation of Eq. (17.74) to the alternate constraint function is Z
I k−1(0) f(x) dS(x). (17.76) The right-hand side here is equivalent to that of Eq. (17.74), so eliminating the surface integral, we have Z
Z ddx \delta[k(x)] |\nabla k| f(x) (17.77) This hold for any test function f(x), so
|\nabla h|. (\delta-function chain rule for scalar constraints) (17.78) Therefore, the transformation just involves the ratio of Euclidean vector-gradient lengths for the two con- straint functions, evaluated in each case at the boundaries. These should both exist for a sensible surface- constraint function, since these act to define the normal vector. 17.1.6.5 Probability Density for Combinations of Random Variables Now we can return to the main point. Suppose we have random variables X1, . . . , XN, with probability
(17.79) Then in the same way as in the result (17.6) for the two-step random walk, we can compute the probability density fy(y) for Y by integrating over all possible values of x, using a delta function to enforce the relation between the variables:
Z dNx f(x) \delta[y −h(x)]. (17.80) Then using Eq. (17.65), we have
I
f(x)
I
f(x1, . . . , xN) v u u t N X j=1 \partial h(x) \partial xj 2 dS, (transformation law, arbitrary function of random variables) (17.81)
desired value Y ], a set that forms a hypersurface S of dimension N −1.
Chapter 17. Stochastic Processes 17.1.6.6 Example: Convolution For example, taking
(17.82) we have |\nabla h| = \sqrt 2, taking x2 = y −x1, and parameterizing the ‘‘surface’’ integral with s = x1 + x2, so that ds = q dx 2 1 + dx 2 2 = \sqrt 2 dx1, (17.83) Eq. (17.81) becomes
Z ds f(x1, y −x1) \sqrt = Z dx1 f1(x1) f2(y −x1), (17.84) where in the last equality we have assumed independence of X1 and X2. This result recovers the convolution (17.7). 17.1.6.7 Example: Quotient of Normal Deviates As a slightly more complicated example, consider two standard-normal deviates X1 and X2, with (separable) joint distribution
2\pi e−(x 2 1 +x 2 2 )/2. (17.85) Now consider the quotient Y = X1/X2 of the two variables, such that the transformation function is
x2 . (17.86) Then the gradient norm is
s x 2 + x 2 x 4 = p x 2 1 + x 2 x 2 = p 1 + y2 |x2| , (17.87) where we have set x1 = yx2, taking x2 as the independent variable. The line element for the ‘‘surface’’ integration is ds = q dx 2 1 + dx 2 2 = p 1 + y2 dx2, (17.88) and thus
Z ds |x2| f(yx2, x2) p 1 + y2 = 1 2\pi Z dx2 |x2| e−x 2
= 1 2\pi 1 + y2 . (17.89) Then we see that the distribution function for the quotient is
(17.90) which is a standard Cauchy distribution. 17.2 Continuous Random Walks: Wiener Process Let’s first define the Wiener process W(t) as a sort of ‘‘ideal’’ random walk with arbitrarily small, inde- pendent steps taken arbitrarily often. That is, the usual random walk is usually taken to be a sequence of random steps of finite average (rms) size, taken after every finite time interval ∆t. Recall from Section 17.1
17.2 Continuous Random Walks: Wiener Process that under fairly reasonable assumptions (such as the existence and finiteness of the one-step variance, and independence of the individual steps), that the central limit theorem guarantees that for long times, the prob- ability density for the walker’s location is Gaussian, independent of the one-step distribution, and the width (standard deviation) increases as \sqrt t. The Wiener process is essentially the idealized limit where ∆t −\rightarrow 0, but where the size of each step decreases as appropriate to maintain the same asymptotic distribution. In this sense, the Wiener process is scale-free, since it has random steps on arbitrarily small time scales, and in fact is a fractal object: a Wiener process with appropriate but arbitrary rescaling (magnification) is still a Wiener process. We choose the Wiener process to correspond to a symmetric random walk, and so W(t) is a normally distributed random variable with zero mean. To fix the scale of the random walk, we choose the variance of W(t) to be simply t. That is, the (rms) width of the distribution is \sqrt t, as is characteristic of a diffusive process. In particular, W(t) has the dimensions of \sqrt t. We can thus write the probability density for W(t) as
\sqrt
(17.91)
emphasize that in view of the central-limit theorem, any simple random walk gives rise to a Wiener process in the continuous limit, independent of the one-step probability distribution (so long as the one-step variance is finite). To get an idea what these look like, 5 and 200 Wiener processes are respectively shown in the two
t Wo(t) -5 t Wo(t) -5
Chapter 17. Stochastic Processes Intuitively, W(t) is a function that is continuous but everywhere nondifferentiable. (Of course, any such statement necessarily includes the proviso that the statement is true except for possibly a set of realizations of zero measure.) Naturally, the first thing we will want to do is to develop the analogue of the derivative for the Wiener process. We can start by defining the Wiener increment
(17.92) corresponding to a time interval ∆t. Again, ∆W is a normally distributed random variable with zero mean and variance ∆t. Note again that this implies that the rms amplitude of ∆W scales as \sqrt ∆t. We can understand this intuitively since it is the variances that add for successive steps in a random walk, not the standard deviations. Mathematically, we can write the variance as
(∆W)2 = ∆t, (17.93) where the double angle brackets \langle \langle \rangle \rangle denote an ensemble average over all possible realizations of the Wiener process. This relation suggests the notion that second-order terms in ∆W contribute at the same level as first-order terms in ∆t, thinking about both of these variables in a small-time expansion of the evolution. In the infinitesimal limit of ∆t −\rightarrow 0, we will write ∆t −\rightarrow dt and ∆W −\rightarrow dW. Then dW(t) is the Wiener differential, which is a fundamental object underlying stochastic calculus, and our analogue of the derivative of the Wiener process. Thought of as a ‘‘signal,’’ it is everywhere discontinuous. Notice that it is somewhat unusual: one sometimes writes a noise process as
dt , (17.94) but this object is singular (i.e., has unbounded variance at any given time t), because as ∆t −\rightarrow 0, ∆W ∆t ∼ \sqrt ∆t ∆t = \sqrt ∆t −\rightarrow \infty. (17.95) It is possible to work with this singular fraction so long as you are careful with it, in the same sense that you can work with the singular delta function. We will tend to stick to the notation of differentials dt and dW(t), but note that while dW is ‘‘zero,’’ it is not ‘‘quite as small as’’ dt. There is, in fact, a deeper connection of dW with the delta function. If we think of dW as a temporal ‘‘noisy’’ signal, the reason dW/dt is singular is that it contains contributions from all frequencies with equal weights—it’s white noise—and that’s the reason why the Wiener process contains random steps on all time scales. However, the total power for such a system, if the power in any band is finite, must be infinite. This is consistent with the fact that dW/dt diverges on average. On the other hand, if we do anything that limits the bandwidth of this signal, such as convolution with a finite function, or using a bandpass filter, it makes the resulting signal finite and well-behaved. Of course, any physical calculation or physical system involves just such a procedure, say, via dissipation through friction. This is exactly analogous to the delta function. The difference is that in the delta function, the frequency components have well-defined relative phases, while for the Wiener process they have effectively random phases. 17.3 It¯o Calculus Now that we have introduced a white-noise process, we will explore the formalism for handling this, par- ticularly for handling the evolution of systems that are driven by white noise. It turns out that adding a white-noise stochastic process changes the basic structure of the calculus for treating the evolution equa- tions. In particular, the usual Riemann integral is undefined for stochastic processes. There is more than one formulation to treat stochastic processes, but we will start out with It¯o calculus,6 which is the one most commonly used in treating quantum systems. We will start by showing how to use this calculus, since the rules are a bit different than what you’re probably used to, and then we will justify the rules of usage. 6The name It¯o is also commonly transliterated as Ito or Itô.
17.3.2 It¯o Rule: Justification¶
17.3 It¯o Calculus 17.3.1 Usage First, let’s review the usual calculus in a slightly different way. A differential equation dy
(17.96) can be instead written in terms of differentials as
(17.97)
try calculating the differential dz for the variable z = ey in terms of the differential for dy as follows:
e\alpha dt −1 . (17.98)
(17.99) This is, of course, the same result as that obtained by using the chain rule to calculate dz/dy and multiplying through by dy. The point here is that calculus breaks up functions and considers their values within short intervals ∆t. In the infinitesimal limit, the quadratic and higher order terms in ∆t end up being too small to contribute. In It¯o calculus, we have an additional differential element dW representing white noise. The basic rule of It¯o calculus is that dW 2 = dt, while dt2 = dt dW = 0. We will justify this later, but to use this calculus, we simply note that we ‘‘count’’ the increment dW as if it were equivalent to \sqrt dt in deciding what orders to keep in series expansions of functions of dt and dW. As an example, consider the stochastic differential equation (SDE)
(17.100) We obtain the corresponding differential equation for z = ey by expanding to second order in dy: dz = ey edy −1 = z dy + (dy)2 . (17.101) Only the dW component contributes to the quadratic term; the result is dz = z
(17.102) The extra \beta2 term is crucial in understanding many phenomena that arise in continuous-measurement processes. 17.3.2 It¯o Rule: Justification We now want to show that the Wiener differential dW satisfies the It¯o rule dW 2 = dt. We already noted above that by definition, the ensemble average of (∆W)2 is equal to ∆t. However, in the infinitesimal limit, we will show that dW 2 = dt holds without the ensemble average. This is surprising, since dW is a stochastic quantity, while dt obviously is not. To show this, consider the probability density function for (∆W)2, which we can obtain by a transforming the Gaussian probability density for ∆W:
\sqrt 2\pi∆t e−(∆W )2/2∆t. (17.103)
X x\in f −1(y) Px(x) dx, (17.104)
17.3.3 Ensemble Averages¶
Chapter 17. Stochastic Processes which is also equivalent to the Frobenius–Peron equation
Z dx′ Px(x′)\delta[y −f(x′)]. (17.105) Then we may write P (∆W)2 = e−(∆W )2/2∆t p 2\pi ∆t (∆W)2 . (17.106) In particular, the mean and variance of this distribution for (∆W)2 are
(∆W)2 = ∆t (17.107) and Var (∆W)2
(17.108) respectively. To examine the continuum limit, we will sum the Wiener increments over N intervals of duration ∆tN = t/N between 0 and t. The corresponding Wiener increments are
(17.109) Now consider the sum of the squared increments N−1 X n=0 (∆Wn)2, (17.110) which corresponds to a random walk of N steps, where a single step has average value t/N and variance 2t2/N 2. According to the central limit theorem, for large N the sum (17.110) is a Gaussian random variable with mean t and variance 2t2/N. In the limit N −\rightarrow \infty, the variance of the sum vanishes, and the sum becomes t with certainty. Symbolically, we can write Z t
N\rightarrow \infty N−1 X n=0
Z t dt′. (17.111) For this to hold over any interval (0, t), we must make the formal identification dt = dW 2. (17.112) (It¯o rule) This means that even though dW is a random variable, dW 2 is not, since it has no variance when integrated over any finite interval. Incidentally, we can also write down a similar expression for dW with itself, but at
∆t < |t −t′|, since the Wiener increments are independent. By a similar argument to the one above, the variance vanishes in the continuum limit—the variance of ∆W(t) ∆W(t′) is bounded above by the variance
have an explicit ensemble average when replacing this product by zero. 17.3.3 Ensemble Averages Finally, we need to justify a relation useful for averaging over noise realizations, namely that
y dW
= 0 (17.113) (It¯o ensemble average) for a solution y(t) of Eq. (17.100). This makes it particularly easy to compute averages of functions of y(t) over all possible realizations of a Wiener process, since we can simply set dW = 0, even when it is multiplied
17.3.4 Correlation Function¶
17.3 It¯o Calculus¶
discrete relation
(17.114) This discrete form here turns out to be the defining feature of It¯o calculus, as we will see. Thus, y(t) depends on ∆W(t −∆t), but is independent of dW(t), which gives the desired result, Eq. (17.113). This gives the important feature of It¯o calculus that makes it useful for computing ensemble averages: at a given time, the state of the noise process and the state of the system are independent. In particular, it is simple to write down an equation for the ensemble average of Eq. (17.100), d
y(t)
=
\alpha(y, t)
dt, (17.115) which we obtain simply by setting dW −\rightarrow 0 in the SDE. This leads us to some common terminology. Any process y(t) satisfying an SDE of the form of
d
y(t)
=
y(t + dt) −y(t)
= 0. (17.116) Any process satisfying this average condition is called a martingale, and is special in that each step in time is unbiased as a random walk. 17.3.4 Correlation Function Now we are in a position to justify the ‘‘whiteness’’ of the noise. Recalling the singular noise signal
dt , (17.117) let’s compute the ensemble average
, (17.118) which is just the correlation function of the noise signal. Note that the ensemble average here can just as well be replaced by a time average, to get the correlation function in the time-averaged sense. If t̸ = t′, we can simply write
= ** dW(t) dW(t′) dt dt′ ++ =
dW(t)
dW(t′)
dt dt′ = 0
(17.119)
=
(dW)2 (dt)2 = 1 dt −\rightarrow \infty. (17.120) Thus we see the divergent behavior. In fact, we can get the normalization from Z dt
= 1, (17.121)
any other point in the integration range. Thus, we can infer that \xi(t) is delta-correlated:
(17.122) (white-noise correlation function) This justifies the notion of dW [equivalently, \xi(t)] as representing white noise, since the power spectrum is the Fourier transform of the correlation function according to the Wiener–Khinchin theorem, which in this case turns out to be a constant function over all frequencies. Note also the peculiarity that everything in
17.3.5 Diffusion¶
Chapter 17. Stochastic Processes 17.3.5 Diffusion The noise term in the It¯o stochastic differential equation (SDE)
(17.123) causes, as you might expect, diffusion of the trajectories y(t). To see this, we need the evolution of the width of the ensemble. Using the It¯o rule
(17.124) we find the mean-square trajectory d
y2 =
=
dt. (17.125) Then defining the ensemble variance by Vy :=
(17.126) we can use
d
(17.127) to write the variance evolution as
= h
i dt. (17.128) (SDE variance evolution) Thus, the variance is affected by gradients of \alpha with y, or ‘‘spatial’’ dependence of the drift coefficient that can stretch or compact the distribution. This is the deterministic component of the variance. The noise part of the equation also contributes the \beta2 term, so that the noise always tends to increase the ensemble variance, thus causing diffusion. 17.3.5.1 Fokker–Planck Equation The evolution of the mean (17.127) and variance (17.128) are equivalent to the mean and variance according to the deterministic Fokker–Planck equation for the probability density f(y, t) (Problem 5.18)
2\partial 2 y \beta2(y, t)f(y, t). (equivalent Fokker–Planck equation) (17.129) In fact, this Fokker–Planck equation turns out to be the correct one to evolve the ensemble density. Recall that the standard form for the Fokker–Planck equation in one dimension is [from Section (5.8.6.1)]
2\partial 2 y D(y, t)P(y, t), (general Fokker–Planck equation) (17.130) where A(y, t) is the drift coefficient, and D(y, t) is the diffusion coefficient. Thus, we identify the stochastic drift coefficient \alpha(y, t) with the Fokker–Planck drift A(y, t), while we identify the squared stochastic coefficient \beta2(y, t) with the diffusion coefficient D(y, t). (For an alternate connection between stochastic trajectories and diffusion-type equations, see Section 17.11.) To prove this, let’s review a couple of concepts regarding probability theory. The conditional prob- ability density P(y, t|y0, t0), is a probability density in y, with P(y, t|y0, t0) dy representing the probability density for finding the particle between y and y + dy at time t, given the particle was at y0 at time t0.
17.3 It¯o Calculus This is distinct from the joint density P(y, t; y0, t0), which is a probability density in both y and y0, where P(y, t; y0, t0) dy dy0 is the probability for finding the particle between y and y + dy at time t and between y0 and y0 + dy0 at time t0. The individual probability densities are given by integrating out the other variable,
Z dy0 P(y, t; y0, t0),
Z dy P(y, t; y0, t0). (17.131) The joint and conditional densities are related by the conditional probability relation, which states that the probability for A and B to occur is the product of the probability for A given that B occured and the probability for B to occur:
(17.132) For Markovian evolution, that is, evolution where the entire state of the system for all future times is determined by the state of the system at the present time, this means that P(y, t) is determined by P(y0, t0). In this case, the conditional density satisfies the Chapman–Kolmogorov equation,
Z dy′ P(y, t|y′, t′)P(y′, t′|y0, t0) (Chapman–Kolmogorov equation) (17.133) which certainly seems a reasonable property of the conditional density: two steps of the evolution of the density may be composed into a single step by integrating over all possible intermediate values. Now to derive the Fokker–Planck equation.7 To do this, consider the evolution of the ensemble average of an arbitrary function g[y(t)], where y(t) is a solution to the SDE (17.123): d
g(y)
=
g′(y) dy + 1 2g′′(y) (dy)2
=
g′′(y) dt
. (17.134) We can obviously rewrite this as \partial t
g(y)
=
\partial 2 y g(y)
. (17.135) The operator acting on g(y) on the right-hand side,
\partial 2 y , (17.136) (generator of the SDE) is often called the generator or infinitesimal generator corresponding to the SDE (17.129), because it ‘‘generates’’ the average change in a function of y over an infinitesimal time step dt. Now let us write out the ensemble average explicitly, using the conditional density P(y, t|y0, t0) for y(t): Z
Z dy P(y, t|y0, t0)
\partial 2 y g(y) . (17.137) Integrating by parts and discarding boundary terms, Z
Z dy g(y)
2\partial 2 y \beta2(y, t)P(y, t|y0, t0) . (17.138) Since g(y) is arbitrary, we may equate the integrands, and thus P(y, t|y0, t0) satisfies an equation with the form of the Fokker–Planck equation:
2\partial 2 y \beta2(y, t)P(y, t|y0, t0). (Kolmogorov forward equation) (17.139) 7As in C. Gardiner, Stochastic Methods: A Handbook for the Natural and Social Sciences, 4th ed. (Springer, 2004) (ISBN: 9783540707127), Section 4.3.5, p, 93.
Chapter 17. Stochastic Processes This equation is called the Kolmogorov forward equation, from which the Fokker–Planck equation (17.130) follows by multiplying through by P(y0, t0) and integrating over y0. The above argument may also be adapted to give an evolution equation in terms of the initial time t0.8 Then writing
\deltat\rightarrow 0 P(y, t|y0, t0 + \deltat) −P(y, t|y0, t0) \deltat = lim \deltat\rightarrow 0 \deltat Z dy′ P(y′, t0 + \deltat|y0, t0)
, (17.140) where the first term takes advantage of P(y′, t0 + \deltat|y0, t0) acting as a normalized distribution in y′, and the second term uses the Chapman–Kolmogorov equation (17.133). Under the assumption that the distributions P(y, t|y0, t0) are continuous functions, the two conditional densities in the difference above can be expanded to lowest order in \deltat, with first- and higher-order terms not contributing in the limit \deltat −\rightarrow 0:
\deltat\rightarrow 0 \deltat Z dy′ P(y′, t0 + \deltat|y0, t0) P(y, t|y0, t0) −P(y, t|y′, t0) . (17.141) Now using the forward equation (17.139) to replace the first conditional distribution, \partial t0P(y, t|y0, t0) = Z dy′
2\partial 2 y′\beta2(y′, t0)P(y′, t0|y0, t0) P(y, t|y0, t0) −P(y, t|y′, t0) . (17.142) Integrating by parts, we have (after discarding surface terms)
Z dy′ P(y′, t0|y0, t0)
y′P(y, t|y′, t0) . (17.143)
\partial 2 y0P(y, t|y0, t0). (Kolmogorov backward equation) (17.144) This peculiar partial differential equation for the initial values y0 and t0 is called the Kolmogorov backward equation. It has a form similar to the Fokker–Planck equation, except for the order of the derivatives and the coefficients, and the minus sign on the time derivative. 17.3.5.2 Multidimensional Fokker–Planck Equation In multiple dimensions, it is fairly straightforward to generalize the equivalent Fokker–Planck equation (17.129). We can start with a multidimensional generalization of the SDE (17.123),
(17.145) where repeated indices are summed. In this case, the equivalent Fokker–Planck equation becomes (Prob- lem 17.2)
2\partial i\partial jDij(x, t)f(x, t), (equivalent Fokker–Planck equation) (17.146) where the diffusion tensor is
(equivalent Fokker–Planck equation) (17.147) in terms of the noise-coupling matrix \beta. 8The beginning of this argument is as in C. Gardiner, op. cit., Section 3.6, p. 55.
17.3.6 Ornstein-Uhlenbeck Process¶
17.3 It¯o Calculus 17.3.6 Ornstein–Uhlenbeck Process As an example of using It¯o calculus, we will consider the Ornstein–Uhlenbeck process, which we can define as the solution of the damped equation driven by a Wiener process:
(17.148) (Ornstein–Uhlenbeck process) This equation is the Langevin equation. To solve this, we write the equation in the form d ye\gammat
(17.149) which we can integrate to obtain
Z t e\gammat′dW(t′), (17.150) or simply
Z t e−\gamma(t−t′)dW(t′). (17.151) The first term is clearly a decaying transient due to the initial condition, while the second is a convolution of the Wiener process with an exponential kernel, effectively smoothing the white noise. Note that since the Wiener differentials are Gaussian random variables, we see from this that the Ornstein–Uhlenbeck process is the sum over Gaussian random variables and is thus itself Gaussian. It is thus sufficient to completely characterize it by computing the mean and autocorrelation function. The mean is simply given by a decaying transient induced by the initial condition,
y(t)
(17.152) since in It¯o calculus we compute ensemble averages by setting dW = 0. The correlation function is given by (taking t′ > t)
y(t) y(t′)
= y 2
Z t ds Z t′ ds′ e−\gamma(t−s)e−\gamma(t′−s′) dW(s) ds dW(s′) ds′
= y 2
Z t ds Z t′ ds′ e−\gamma(t−s)e−\gamma(t′−s′)\delta(s −s′) = y 2
Z t ds e−\gamma(t−s)e−\gamma(t′−s) = y 2
Z t ds e−2\gammas = y 2 0 −1 2\gamma
(17.153) If we regard y0 as an belonging to an ensemble of variance 1/2\gamma, or if we consider the limit t −\rightarrow \infty(with t −t′ fixed), we can ignore the transient part and take the correlation function to be
y(t) y(t′)
= 1
(17.154) (Ornstein–Uhlenbeck correlation)
Uhlenbeck process, although corresponding to Gaussian noise, does thus not have independent increments at different times. We see that the damping introduces a ‘‘memory’’ in the dynamics. Further, since the correlation function is exponential, we immediately see that the power spectral density for the Ornstein– Uhlenbeck process is Lorentzian, and thus decays asymptotically as \omega−2, and thus corresponds to ‘‘1/f 2’’ noise. But also notice that the correlation function is finite for all times, and thus the damping, which has cut off the high frequencies, has made the noise bounded and well-behaved.
Chapter 17. Stochastic Processes 17.3.6.1 Brownian Motion The Ornstein–Uhlenbeck process is a model for Brownian motion, if we use it to model the velocity of a particle subject to friction and to frequent ‘‘kicks,’’ as from collisions with many background-gas atoms:
(17.155) We are assuming that we are dealing with a scaled velocity such that the units come out right. Again, we can equivalently write this in the (possibly) more familiar form
(17.156) so that we have the usual equation of motion for the velocity, but driven by a white noise ‘‘force’’ \xi(t). The Ornstein–Uhlenbeck process also corresponds to a white-noise voltage signal \alpha\xi(t) passing through a low-pass filter, with no connection at the output V (t). I
R C V (t) To see this, note that the current I flowing through the resistor is
R , (17.157) and the voltage across the capacitor is related to the current by
RC , (17.158) which we can write as
(17.159) where \gamma = 1/RC, as we expect for the low-pass filter. Thus, the output of the low-pass filter corresponds to a scaled Ornstein–Uhlenbeck process, and thus a physical signal despite the idealized input. Note that from Eq. (17.154), an Ornstein–Uhlenbeck process of the form
(17.160) has \langle \langle y2\rangle \rangle = 1/2\gamma, and thus an rms fluctuation yrms =
(17.161) If we write Eq. (17.159) in the form
(17.162) and then we let t −\rightarrow t/(\alpha\gamma)2 and W −\rightarrow W/\alpha\gamma, we have dV = −1
(17.163) which is in the form of (17.160). Thus, the rms output voltage of the low-pass filter is
r\gamma 2 , (17.164) (amplitude of filtered white noise) for an input signal of \alpha\xi(t). Note that the rms output voltage increases as the square root of the filter
transmitted power is proportional to the width of frequency band passed by the filter. Evidently, \alpha has the dimensions of V/ \sqrt Hz.
17.4 Stratonovich Calculus 17.4 Stratonovich Calculus The main alternative to It¯o calculus is Stratonovich calculus, which we will introduce primarily to gain more insight into It¯o calculus. Although Stratonovich calculus has some aesthetically nice features, it is often easier to perform calculations in It¯o form, and we will mainly stick to It¯o equations in our discussion of quantum measurement. Consider the deterministic ODE
(17.165) This is the continuous limit of the discrete relation
(17.166) where \tau is an arbitrary time in the range [t, t+∆t]. This is because the formal (implicit) solution of (17.165) is given by the Riemann integral
Z t \alpha(y(t′), t′) dt′. (17.167) The Riemann integral is approximated by successively finer refinements of discrete ‘‘rectangle’’ areas of width ∆t, where the area of each rectangle is determined by the value of the integrand at any point within the interval ∆t. Returning to the SDE
(17.168) we must also interpret the solution of this equation in terms of the implicit integral
Z t
Z t \beta(y(t′), t′) dW(t′). (17.169) The first integral is an ordinary Riemann integral, but the second is of the form Z t
Z t \beta(y(t′), t′) dW dt′ dt′. (17.170)
to save this is that in the successive finite approximations, if you consistently pick the same point within each interval of equal length ∆t, the integral is defined. However, the result that you get by evaluating the integral will depend on your choice. Of course, if \beta is constant or even a smooth function of time—in the case of additive noise—then this won’t be a problem, since the result amounts to integrating dW to get W(t). The problem arises in the case of multiplicative noise when \beta is a function of y. Thus, when we regard the SDE (17.168) as the continuum limit of the finite-difference equation
(17.171) where \tau \in [t, t + ∆t], we obtain a different limit depending on where in the interval we choose \tau. In view of
Since we expect different results according to what calculus we intend to use to solve the equation, we must use a notation to distinguish the calculus that goes with the SDE. Thus, we will write as usual
(17.172) (notation: It¯o SDE) to denote an It¯o SDE, while we will use the special notation
(17.173) (notation: Stratonovich SDE) to refer to a Stratonovich equation.
17.4.1 Example: Stochastic Integration¶
Chapter 17. Stochastic Processes 17.4.1 Example: Stochastic Integration To illustrate the consequences of this choice, we will compute the sample It¯o integral I = Z t t0 W(t′) dW(t′), (17.174) and compare it to the Stratonovich integral J = Z t t0 W(t′) ◦dW(t′), (17.175) of the same form, to see that they indeed give different results. Note that if these were ordinary Riemann integrals, we would simply have Z t t0
f 2(t) −f 2(t0) (17.176) for a sufficiently well-behaved function f(t). The It¯o integral follows from the continuum limit of the N-step
I = lim N\rightarrow \infty N−1 X j=0 W(tj)∆W(tj) = lim N\rightarrow \infty N−1 X j=0 W(tj) [W(tj+1) −W(tj)] = lim N\rightarrow \infty N−1 X j=0 [W(tj) + W(tj+1) + W(tj) −W(tj+1)] [W(tj+1) −W(tj)] = lim N\rightarrow \infty N−1 X j=0 W 2(tj+1) −W 2(tj) −lim N\rightarrow \infty N−1 X j=0 [W(tj+1) −W(tj)]2 = 1 W 2(t) −W 2(t0) −lim N\rightarrow \infty N−1 X j=0 [∆W(tj)]2 = 1 W 2(t) −W 2(t0) −1 Z t t0 [dW(t′)]2 = 1 W 2(t) −W 2(t0) −1 2(t −t0). (17.177)
so that I = Z t t0
Z t t0 d(W 2) −dt = 1 W 2(t) −W 2(t0) −1 2(t −t0). (17.178) We thus see how the It¯o rule enforces the choice of approximating integration intervals by the beginning point of the interval.
17.4.2 It¯o-Stratonovich Conversion¶
17.4 Stratonovich Calculus to s = 0 and Stratonovich to s = 1/2: Js := lim N\rightarrow \infty N−1 X j=0 W(tj+s)∆W(tj) = lim N\rightarrow \infty N−1 X j=0 W(tj+s) [W(tj+1) −W(tj)] = lim N\rightarrow \infty N−1 X j=0 W(tj) [W(tj+1) −W(tj)] + lim N\rightarrow \infty N−1 X j=0 [W(tj+s) −W(tj)] [W(tj+1) −W(tj)]
N\rightarrow \infty N−1 X j=0 [W(tj+s) −W(tj)] [W(tj+1) −W(tj)]
N\rightarrow \infty N−1 X j=0 [W(tj+s) −W(tj)] [W(tj+1) −W(tj+s)] + lim N\rightarrow \infty N−1 X j=0 [W(tj+s) −W(tj)]2 (17.179) The second term corresponds to the continuum limit of a product of independent Wiener increments, which vanishes according to the last argument of Section (17.3.2). The last term is the sum of squared, independent Wiener increments corresponding to time intervals s∆t, and is thus given by s(t −t0). Thus,
W 2(t) −W 2(t0) + s −1 (t −t0). (17.180) In particular, the Stratonovich integral is
W 2(t) −W 2(t0) . (17.181) Note that this is exactly the same result as if we had just used ordinary calculus, so that in Stratonovich
We will prove this after we see how to convert between It¯o calculus and Stratonovich calculus. 17.4.2 It¯o–Stratonovich Conversion It¯o and Stratonovich SDEs in general give different results for the same coefficients, but in what sense are they equivalent? That is, how do we convert between It¯o and Stratonovich SDEs? Suppose we have the It¯o SDE
(17.182) and the Stratonovich SDE
(17.183) Then what we will show is that these two SDEs are equivalent if and only if
(17.184) (It¯o–Stratonovich conversion) Clearly, this distinction only matters in the case of multiplicative noise. To show this, recall that the It¯o SDE is the continuum limit of the discrete relation
(17.185) while the Stratonovich SDE is the continuum limit of
(17.186)
17.4.3 Stratonovich Calculus and the Chain Rule¶
Chapter 17. Stochastic Processes Now noting that
h
i + O(∆t), (17.187)
and we are dropping terms of order ∆t and higher since this expression will be multiplied by ∆W. We can thus use this result to write
(17.188) In the continuous limit, we can write ∆W(t) ∆W (1/2)(t) −\rightarrow dt 2 , (17.189) since [∆W (1/2)]2 −\rightarrow dt/2 with certainty, and the product of ∆W (1/2) from two different time intervals will converge to zero with certainty. Thus, the continuum limit of (17.186) is the It¯o-form SDE dy =
(17.190) Since this is the continuum limit of the same equation as the Stratonovich form (17.183), we can thus conclude that the It¯o and Stratonovich forms are equivalent in the case
(17.191) Thus, when writing down an SDE, we again see that it is crucial to specify which calculus it assumes, since the solution would otherwise be ambiguous. 17.4.3 Stratonovich Calculus and the Chain Rule Recall that the It¯o equation
(17.192)
2f ′′(y)(dy)2 =
2f ′′(y) \beta2
(17.193)
We will now show that in Stratonovich calculus, the ‘‘extra’’ \beta2 term does not appear, and thus the usual
Thus, consider the Stratonovich equation
(17.194)
dy = \alpha + 1
(17.195) and now we can use It¯o rules to accomplish the transformation:
2f ′′(y)(dy)2 = f ′(y) \alpha + 1
+ 1 2f ′′(y)\beta2
(17.196)
17.4.4 Comparison¶
17.4 Stratonovich Calculus Now transform this back into Stratonovich form: dz = f ′(y) \alpha + 1
+ 1 2f ′′(y)\beta2 −1
(17.197) Noting that
f ′(y)\partial y, (17.198) we can write the last dt term as −1
= −1 2f ′′(y)\beta2 −1
(17.199) Thus, this term, which we obtained from switching from It¯o to Stratonovich form, cancels the other two dt terms that involve \beta, leaving
(17.200)
Thus, Stratonovich calculus obeys the usual chain rule, and we have no need for terms of order dW 2. 17.4.4 Comparison We will now summarize the differences between It¯o and Stratonovich SDEs, and then explain why we tend to favor It¯o calculus. For It¯o calculus: • The rules of stochastic integration are slightly more complicated than for ordinary integration, since
• The solution y(t) of an It¯o SDE and the driving Wiener process dW(t) are statistically independent at equal times, so that ensemble averages are simply computed by setting dW = 0. • An It¯o SDE is ‘‘natural’’ as the continuum limit of an evolution constructed by a discrete-step process, since
(17.201) is the continuous limit of
(17.202) On the other hand, for Stratonovich calculus: • The rules of stochastic integration are those of ordinary Riemann integration, since the usual chain
• The solution y(t) of a Stratonovich SDE and the driving Wiener process dW(t) are not statistically independent at equal times. This is clear from the above It¯o–Stratonovich conversion, since setting dW = 0 in a Stratonovich SDE does not give the same result as setting dW = 0 in the equivalent It¯o SDE. In fact, the easiest rule for computing an ensemble average is to convert the SDE to It¯o form and then set dW = 0. • A Stratonovich SDE is ‘‘natural’’ as the idealization of a physical noise process in the following sense. If one models a stochastic system as being driven by a physical noise of finite bandwidth and bounded variance, then the normal rules of calculus apply. For example if dO represents an Ornstein–Uhlenbeck process, then we could model a system by the SDE
(17.203)
Chapter 17. Stochastic Processes Then if you take the white-noise limit for the driving process (\gamma −\rightarrow 0), then dO goes over to dW, but because the ODE always obeyed the rules of ordinary calculus, the white-noise limit
(17.204) should be interpreted as a Stratonovich SDE. Also note that most proofs, including the construction of Stratonovich calculus, are usually proved in It¯o calculus, so its advantages tend to outweigh its peculiarities. For handling quantum measurements, we will often want to compute ensemble averages to obtain unconditioned master equations, and we will also in general construct continuous measurements as limits of discrete processes of the form (17.202). Thus, we will virtually always use It¯o-form SDEs to handle continuous quantum measurements. 17.5 Poisson Process Recall that the Poisson probability distribution of mean \lambda is
n! , (17.205) (Poisson distribution) where n is a nonnegative integer. The variance of the Poisson distribution is equal to the mean \lambda. The Poisson distribution models the number of independent random events that occur in a given interval of time, such as the number of cars that arrive at an intersection or the number of atoms that decay in a large ensemble. Poisson random variation is also responsible for shot noise, which occurs as noise in electrical current due to random fluctuations in the rate at which electrons flow through a device, or as noise in the detected intensity of classical light due to the random detection times of individual photons. In general, we can speak of a rate at which events occur by setting \lambda = \Gamma ∆t for finite time interval ∆t, where \Gamma is the mean rate of occurence (also called the intensity of the Poisson process). Then
n! . (17.206) Note that the Poisson distribution implies an exponential waiting time for the first event, because the probability for the event to occur after waiting a time ∆t is given by setting n = 0 in the above probability function:
(17.207) Then according to our interpretation, this probability is related to the probability density P ′(t) for the time of first occurence by e−\Gamma ∆t = Z \infty ∆t P ′(t) dt, (17.208) so that
(17.209) Thus, Poisson random variables are intimately connected with exponential-decay processes, such as sponta- neous emission from an atom prepared in the excited state.
higher). The probability for a single event occuring during an interval of duration dt thus becomes \Gamma dt, with no events occuring otherwise. We can denote this by the infinitesimal random variable dN(t)—the Poisson process—which has an ensemble mean
dN(t)
= \Gamma dt. (17.210) (Poisson process: ensemble mean) In the standard Poisson process N(t), the intensity \Gamma is a constant that characterizes the process. Thus, in general, when writing down a Poisson process, you must always also specify the intensity, which is
17.5.1 The Poisson Process Implies the Poisson Distribution¶
17.5 Poisson Process not specified in the notation dN in contrast to the Wiener process dW (where there is no freedom to specify the moments). In generalizations of the Poisson process to time- or state-dependent intensities (see Section 17.5.2), an explicit specification of the intensity is even more critical. Again, as an integer-valued differential random variable, dN can take on only the values 0 and 1, where the value of 1 occurs with probability equal to the mean. Because dN(t) \in {0, 1}, it immediately follows that dN 2 = dN, (17.211) (Poisson-process property) so that the ensemble-averaged variance is equal to the ensemble mean,
dN 2 =
dN
= \Gamma dt, (17.212) as we expect for a Poisson-distributed random variable. Note that the infinitesimal variance here is just the second moment, since the square of the mean is O(dt2). In another view, note that dN(t)/dt is zero except in isolated intervals of length dt, where the value is 1/dt. Thus, we can write this form of the Poisson process as the sum of delta functions, dN(t) dt = X j \delta(t −tj), (17.213) if the events occur at times tj. In view of our discussion above, ∆tj := tj+1 −tj is a random variable with probability density
(17.214) since the waiting time until the next event is given by the exponential distribution. 17.5.1 The Poisson Process Implies the Poisson Distribution If we take the Poisson process dN(t) of intensity \Gamma (i.e., of mean \Gamma dt) as the fundamental object, we should also be able to derive the Poisson distribution for the frequency of events in finite time intervals. That is, we can show that the Poisson distribution arises if there is a constant probability per unit time of a single event occuring during an arbitrarily short time interval. Specifically, defining the time integral of the Poisson process, ∆N := Z t+∆t t dN(t′) (17.215) we can ask, what is the probability distribution for ∆N? For a given value n of ∆N, this means that during exactly n infinitesimal intervals, dN(t) took on the value unity, while it took on the value of zero during the remaining intervals. The probability of doing so is the product of three factors: 1. The probability for having exactly n such events, one in each of n particular time intervals: (\Gamma dt)n. 2. The number of ways to distribute the n events among all such intervals. In the time interval [t, t+∆t), there are ∆t/dt such time intervals, and so there are (∆t/dt)n ways to distribute n events among all possible time intervals. But we divide by n! since we take the n events to be indistinguishable, so we don’t overcount, so the total factor is (∆t/dt)n/n!. 3. The probability of having zero events in all other time intervals. Again, there are (∆t/dt) total time intervals, and the probability of zero events in a given time interval is (1−\Gamma dt), so the total probability is
M\rightarrow \infty 1 −\Gamma ∆t M M = e−\Gamma ∆t. (17.216) We are being cavalier in slinging around factors of dt, but our manipulations here are equivalent to the ‘‘correct’’ approach of using finite, small subintervals \deltat, where we neglect n compared to ∆t/\deltat, and we
17.5.3 White-Noise Limit¶
Chapter 17. Stochastic Processes neglect the probability that two events end up in the same subinterval. Both of these approximations are appropriate (and exact) in the continuum limit. The total probability is thus
∆t dt n 1 n!e−\Gamma ∆t
n!
n! , (17.217) which is the Poisson distribution, where again the mean is \lambda = \Gamma ∆t. Thus, the Poisson distribution results when there are many ‘‘trials’’ (short time intervals), where there is a vanishingly small probability of ‘‘success’’ (an event occurrence) in each trial.
is piecewise constant, taking on only constant, nonnegative integer values in intervals of the form [tj, tj+1), and with ‘‘jumps’’ of unit size occurring at the times tj, so that N(t) is also nondecreasing. Further, N(t) must follow a Poisson distribution (with mean \Gammat in the homogeneous case, where \Gamma is the intensity). A more general process that relaxes the requirement of the Poisson distribution is called a counting process, and a jump process additionally relaxes the requirement of unit jumps (and therefore of monotonicity, if negative jumps are allowed). 17.5.2 Inhomogeneous Poisson Process and State Dependence In general, the rate \Gamma of event occurence may depend on time, either explicitly or via dependence on the state y(t) of the system. We can handle this by noting that according to our construction for the Poisson distribution above, then if X1 and X2 are Poisson-distributed random variables with means \lambda1 and \lambda2, then X1 + X2 is a Poisson random variable of mean \lambda1 + \lambda2. This statement amounts to agglomerating two adjacent time intervals of different duration in the above derivation of the Poisson distribution. Then if \Gamma is time-dependent, we can subdivide the interval [t, t + ∆t) into sufficiently fine increments such that \Gamma is constant over each increment, and sum them to find that the number of events occuring in the interval [t, t + ∆t) is still a Poisson variable with mean
Z t+∆t t \Gamma(t′) dt′. (17.218) (inhomogeneous Poisson process: mean) The variance is also of course just ¯\lambda. However, while we can define the inhomogeneous Poisson process as above, a generalization to a process with state-dependent intensity \Gamma(y), where y(t) is some process driven by dN(t), is not a Poisson process: the argument above for the inhomogeneous process does not apply, because dN(t) is no longer statistically independent at different times.10 Since in this case N(t) is no longer Poisson-distributed, it is more proper to refer to it as a counting process, as we defined in Section 17.5.1. 17.5.3 White-Noise Limit Consider the scaled Poisson process dy = dN \sqrt \Gamma , (17.219) where
(17.220) 9see, e.g., Rama Cont and Peter Tankov, Financial Modelling with Jump Processes (Chapman & Hall/CRC, 2004), Section 2.5.3, p. 48. 10Note, however, that this is still sometimes referred to as a ‘‘state-dependent Poisson process.’’ See, e.g., Edoardo Daly and Amilcare Porporato, ‘‘Intertime jump statistics of state-dependent Poisson processes,’’ Physical Review E 75, 011119 (2007) (doi: 10.1103/PhysRevE.75.011119).
17.5 Poisson Process It may be that in a given system, the rate \Gamma of Poisson events is much faster than the processes of physical interest. In such a case, we can ignore the discreteness of the events, and coarse-grain the dynamics to approximate the Poisson events by white noise. Note in particular that the mean of dy is
\sqrt \Gamma = \sqrt \Gamma dt, (17.221) while the variance is
\Gamma
\Gamma = dt. (17.222) Thus, if events occur rapidly on time scales of interest—that is, we only measure ∆y over time intervals ∆t ≫1/\Gamma, by the central limit theorem, we may effectively regard dy as a Gaussian random variable of mean \sqrt \Gamma dt and variance dt. In particular, the Poisson process corresponds to a random walk of steps of length \sqrt \Gamma, in one direction only, taken at random times, as we can see by writing
\sqrt \Gamma Z t
\sqrt \Gamma Z t dN(t′) dt′ dt′, (17.223) and recalling from Eq. (17.213) that dN/dt is a sum of delta functions. After many events, the result is the same as a biased random walk, and the central limit theorem again guarantees an asymptotic Gaussian probability density. In this limit, it thus is a good approximation to write dy = \sqrt \Gamma dt + dW. (17.224) Thus, in the limit where Poisson events occur at a very large rate \Gamma, we can make the formal replacement dN −\rightarrow \Gamma dt + \sqrt \Gamma dW, (17.225) (white-noise limit of Poisson process) to approximate the Poisson process with a mean drift plus white noise. 17.5.3.1 Shot Noise As an example, let’s consider shot noise of an electrical current. Let Q denote the total charge that has crossed a certain point along a wire. The current is given by I = dQ dt , (17.226) so that dQ = I dt. (17.227) We will model the current as a stream of independent electrons of charge −e with Poisson arrival times, so that dQ = −e dN, (17.228) with \langle \langle dN\rangle \rangle = \Gammadt as usual. Then equating the two expression for dQ, we find I dt = −e dN. (17.229) The mean current is then given by taking the ensemble average of this relation, so that
(17.230) Frequencies of interest for measuring electrical currents generally range from ns to s, whereas \Gamma ∼1019 s−1. Thus, the white-noise approximation is quite appropriate, and thus dQ \approx −e\Gamma dt −e \sqrt \Gamma dW =
I
\sqrt \Gamma dW =
I
dt − p
(17.231)
Chapter 17. Stochastic Processes Thus, we see that due simply to the discreteness of charge, a mean current \langle \langle I\rangle \rangle is accompanied by white noise of amplitude
\sqrt \Gamma = p
(17.232) Note that the SI units of the current noise amplitude are in A/ \sqrt Hz, since when multiplied by dW/dt, which has dimensions 1/\sqrts, the noise amplitude takes the dimensions of current. The alternate way to view this is that physically, the noise is always bandwidth-limited (as in the RC model of the Ornstein–Uhlenbeck process), and so the filtered noise amplitude is given by multiplying the above noise amplitude by the square root of the circuit bandwidth. More explicitly, the above white noise corresponds to a uniform spectral density of signal power. According to Eq. (17.164), an input signal \alpha\xi(t) corresponds to rms fluctutations of \alpha p \gamma/2 at the output of a low-pass filter. Thus, the rms current fluctuation through a low-pass filter due to shot noise is given by
p
r\gamma 2 , (17.233) where \gamma = 1/RC is the angular cutoff frequency for the low-pass filter. The equivalent-power bandwidth ∆\nu is defined as the bandwidth of the ‘‘brick wall’’ filter (with flat response up to a sudden cutoff at frequency ∆\nu, where ∆\nu is a frequency in Hz, not an angular frequency). The low-pass filter function for power transmission is (up to a constant factor, set by requiring the dc transmission to unity) the Fourier transform of the Ornstein–Uhlenbeck correlation function e−\gamma\tau, so that we may write the transmission function as
\gamma2
(17.234) Obviously, \gamma is the (angular) ‘‘3 dB’’ frequency, or the frequency where the transmission drops to 1/2 the dc value:
2\pi . (17.235) Since Z \infty
(17.236) where the second result applies to the brick-wall filter, we can write
2 f3 dB, (17.237) and thus the shot-noise magnitude is
p
(17.238) (shot-noise amplitude) This expression applies to filters beyond the low-pass filter, so long as the appropriate equivalent-power bandwidth is used in this relation. Thus, a 1 A average current detected in a 1 MHz bandwidth has an rms noise current of 0.57 µA, a fluctuation at under the ppm level. Shot noise clearly gets much worse for smaller currents: for the same bandwidth, an average current of 1 µA has fluctuations of 0.57 nA rms, or 0.057% relative noise, and an average current of 1 pA has fluctuations of 0.57 pA rms, or 57% relative noise. Note that this model assumes the independence of electrons, and gives an appropriate result, e.g., for semiconductor junctions, but not in metallic-wire circuits, where long-range correlations between electrons suppress shot noise.11 Essentially, this is just because of Coulomb interactions between electrons, which causes them to antibunch: in a conducting, crystalline lattice, it is energetically favorable to have two electrons in different lattice sites, as compared to having them occupy the same lattice site. Probabilities of seeing more than one electron pass in a metallic wire in a short time interval are thus suppressed compared to the Poissonian expectation. Of course, we can adapt this result to the case of optical shot noise. Instead of an electrical current, we have a detected optical power P. We can treat the photon arrival times as independent in the case of coherent 11Paul Horowitz and Winfield Hill, The Art of Electronics, 2nd ed. (Cambridge, 1989), pp. 431-2.
17.6.1 Free-Particle Limit¶
17.6 Particle Subject to a Stochastic Force light (recall that in a coherent state the photon-number occupation probabilities are Poisson-distributed). The rms fluctuations for a mean optical power \langle \langle P\rangle \rangle detected in an equivalent-power bandwidth ∆\nu are then given by
p
(17.239) (optical shot-noise amplitude) where the photon energy ¯h\omega plays the role of the electron charge. For 780 nm light in a 1 MHz detection bandwidth, a 1 W power has fluctutations of 0.71 µW, less than the ppm level. For the same bandwidth, a 1 µW mean power has fluctuations of 0.71 nW (0.071%), and a 1 pW mean power has fluctuations of 0.71 pW (71%). Of course, the relative fluctuations will be even larger for thermal (bunched) light, and smaller for antibunched light (nonclassical light, as for resonance fluorenscence of a two-level atom). 17.6 Particle Subject to a Stochastic Force As an important example to tie together the stochastic-process concepts that we have developed so far, we will treat in some depth the dynamics of a particle subject to a stochastic force, together with damping and other external forces. In particular consider the equations of motion dx = p m dt dp = F(x, p; t) −\gammap
(damped, stochastically forced particle) (17.240) with F defining an external force, \gamma the damping coefficient, and \sigma defining the magnitude of the stochastic force. This is for a particle in one spatial dimension, but the generalization to multiple degrees of freedom is straightforward. Also, note that \gamma could depend on x and p (and even t), as could \sigma. Note that if \sigma depends on p, then the particle is coupled to multiplicative noise (Section 17.4), and it is important to specify which calculus to use in treating the SDE system. As discussed in Section 17.4.4, if the noise arises as the continuum limit of a physical noise process of finite bandwidth, then Eqs. (17.240) should be interpreted in the Stratonovich sense. On the other hand, if the noise arises due to a process like spontaneous emission, where the emission probability depends on the present state (or arises in terms of a discrete Poisson process, which goes over to white noise in the limit of frequent emission events, as in Section 17.5.3), then the equations of motion should be interpreted in the It¯o sense. For concreteness we will treat this system in the It¯o sense, although it is useful to remember that as we make assumptions on \sigma in the discussion below, the importance of this distinction drops away. From Eqs. (17.145) and (17.146), we can write out the equivalent Fokker–Planck equation for the probability density f(x, p; t) of the particle state as
F(x, p; t) −\gammap f + 1 2\partial 2 p \sigma2f. (equivalent Fokker–Planck equation) (17.241) Note that any dependence of F, \gamma, and \sigma on p enforces the ordering of the derivative operators as shown; in the case where there is no such momentum dependence, the derivative operators \partial p can commute with these functions and operate only on the probability density. But again notice that dependence on x causes no problem; for example, the diffusion term can have the form \sigma2(x) \partial 2 p f in the case that \sigma depends only on position. 17.6.1 Free-Particle Limit One of the simplest limits of the above equations of motion is to assume vanishingly small damping (\gamma −\rightarrow 0) and external force (F −\rightarrow 0). In this case we can integrate the momentum equation to obtain
Z t dt′ \sigma(x, p; t′) dW(t′). (17.242)
Chapter 17. Stochastic Processes Then putting this into the position equation and integrating, we find
m + 1 m Z t dt′ Z t′ dt′′ \sigma(x, p; t′′) dW(t′′). (17.243) It is difficult to carry this further in the general case. However, in the special case where \sigma is constant, we can characterize the transport of x quite precisely. In this case, we obtain
m Z t dt′ W(t′), (17.244)
for obvious reasons it goes by the names integrated Brownian motion and Brownian area (this is also the iterated integral I10 that we discuss in Section 27.3), and here we will write it as
Z t dt′ W(t′). (17.245) (integrated Brownian motion) Since I(t) is a linear combination of Gaussian increments, it is itself Gaussian, so we can characterize it fully by computing its mean and correlation. The mean is simple, since in It¯o calculus we just set dW = 0:
I(t)
= 0. (17.246) To compute the variance (see also Problem 27.1), we can write DD I2(t) EE =
dt′ Z t′ dW(t′′)
. (17.247) Then using the integral identity R t 0 dt′ R t′
R t 0 dt′′ R t t′′ dt′g(t′, t′′) to interchange the order of integration, we have DD I2(t) EE =
dW(t′′) Z t t′′dt′
=
(t −t′′) dW(t′′)
. (17.248)
DD I2(t) EE =
(t −t′′)2 dt′′ ++ = t3 3 . (17.249) This implies a covariance
I(t)I(t′)
= 1 min{t, t′} 3, (covariance of integrated Brownian motion) (17.250) because for example if t′ > t, then the Wiener increments associated with the time interval (t, t′) are uncorrelated with any increments associated with the interval (0, t), so the time interval (t, t′) does not contribute to the correlation function. Thus, the stochastically driven, undamped particle (17.244) tends to a Gaussian distribution in x, whose width grows as (\sigma/m)t3/2, which is one power of t higher than the t1/2 dependence we expect for diffusive growth in momentum.
17.6.2 Thermal Equilibrium¶
17.6 Particle Subject to a Stochastic Force 17.6.2 Thermal Equilibrium Now consider the particle equations (17.240) with damping at rate \gamma but no external force. This could correspond, for example, to a particle moving (in one dimension) in a static fluid at some temperature T. For simplicity we will assume constant temperature, noise coefficient \sigma, and damping rate \gamma. First, note that the momentum equation in this case is basically a scaled Ornstein–Uhlenbeck process:
(17.251)
(17.252) which has the form of the Ornstein–Uhlenbeck process in Eq. (17.148). Then the results following from that analysis carry through here with the replacements \gamma −\rightarrow \gamma/\sigma2 and t −\rightarrow \sigma2t. For example, the mean from Eq. (17.152) is only transient,
p(t)
(17.253) (momentum mean) and the correlation function (17.153) becomes
p(t) p(t′)
= p 2 0 −\sigma2 2\gamma
(17.254) (momentum correlation) where the first term is again transient and the second persists in steady state, and t, t′ \ge 0. In thermal equilibrium, we have 2m
v2 =
p2 2m = 1 2kBT. (17.255)
p 2m\gammakBT, (17.256) (temperature related to diffusion) which fixes the noise coefficient in terms of the temperature and damping. Thus, for example, Eq. (17.254) becomes
p(t) p(t′)
= mkBT e−\gamma|t′−t| (17.257) in terms of temperature, after ignoring any transients. The position equation (17.240) then says that the position is the (scaled, shifted) integrated Ornstein– Uhlenbeck process:
Z t dt′ p(t′). (17.258) For the ‘‘pure’’ Ornstein–Uhlenbeck process (17.148),
(17.259) we can define the integrated process by
Z t dt′ y(t′). (17.260) Then clearly G(t) is again a Gaussian process, with
G(t)
= 0. (17.261)
17.6.3 Strong-Damping Limit: Brownian Motion Revisited¶
Chapter 17. Stochastic Processes Computing the correlation function, we start with
G(t) G(t′)
= Z t ds Z t′ ds′
y(s) y(s′)
= 1 2\gamma Z t ds Z t′ ds′
. (17.262)
but keeping the transient term, which is important when integrating the correlation function. Then if t′ \ge t
G(t) G(t′)
= 1 2\gamma Z t ds "Z s
Z t′ s ds′ e\gamma(s−s′) # −1 2\gamma Z t ds Z t′
(17.263) Carrying out the integration and simplifying gives
G(t) G(t′)
= 2\gamma3 h
- t \gamma2 . (17.264) Again, using the requirement that the correlation function should be symmetric in t and t′, we find
G(t) G(t′)
= 2\gamma3 h
- min{t, t′} \gamma2 . (integrated-Ornstein–Uhlenbeck correlation) (17.265)
Returning to the transport of the particle in the fluid, we are interested in the asymptotic growth of the variance of G(t). For large times, Eq. (17.265) gives an asymptotic growth DD G2(t) EE ∼t \gamma2 (17.266)
R t 0 dt′ p(t′)/m and Eq. (17.254) for the momentum correlation function gives
x2(t)
∼ \sigma m\gamma 2 t = 2kBT m\gamma t. (17.267) (position diffusion) That is, x(t) grows diffusively with diffusion coefficient (\sigma/m\gamma)2. Thus, the damping ‘‘transfers’’ the diffusive behavior from momentum (which is driven by the Wiener process) to position. We will see this more directly in the following section. 17.6.3 Strong-Damping Limit: Brownian Motion Revisited Back in Section 17.3.6.1, we commented that the Ornstein–Uhlenbeck process is a model for Brownian motion—as we just showed, a damped and stochastically driven momentum leads to diffusion in position. At the same time, the Wiener process W(t) itself is often referred to as ‘‘Brownian motion,’’ although in the physical model it enters via the force, not directly in the particle’s position. Here we will make the connection more directly of the stochastic force leading to diffusion in position, in the case of strong damping. Starting again with the equations of motion (17.240), dx = p m dt dp = F(x, p; t) −\gammap
(17.268) with the same comments about the (x, p)-dependence of \gamma and \sigma, we will begin by assuming a large damping coefficient \gamma. In the limit of large damping, we will assume that the momentum is always approximately in
17.6 Particle Subject to a Stochastic Force equilibrium with respect to the position state, so that we may perform an adiabatic elimination of momentum
adiabatic approximation p \approx 1 \gamma
dt . (17.269) In assuming large damping, we are essentially assuming that the damping term ‘‘absorbs’’ the entire effect of the external and fluctuating forces, provided we coarse-grain on time scales much longer than \gamma−1. The difficulty here is that dW/dt is divergent and fluctuates on all time scales. We will proceed for now, but we will soon see that we need to be somewhat more careful with this approximation.
dx \approx F
m\gamma dW. (particle motion, adiabatic approximation) (17.270)
Brownian motion for W(t). And again appealing to Eqs. (17.145) and (17.146), the equivalent Fokker–Planck equation for this approximate SDE is
F
2\partial 2 x \sigma2 m2\gamma2 f. (Smoluchowski equation) (17.271) Here, F, \gamma, and \sigma can still depend on position and momentum; however, the dependence on momentum
damping ‘‘concentrates’’ the momentum around the mean value. This diffusion equation is the counterpart of the original Fokker–Planck equation (17.241) after momentum is adiabatically eliminated. In the case of coupling to a uniform bath in thermal equilibrium, \sigma is constant and given by Eq. (17.256), such that
m\gamma \partial 2 x f, (17.272) (Smoluchowski equation) after also assuming constant damping \gamma. This is called the Smoluchowski equation.12 Note that in the form (17.272), at any given temperature, the steady-state solution is independent of
write
x f = 0 (17.273) in steady state. Removing one derivative gives
kBT f, (17.274) where we can see that any constant of integration vanishes by comparing this equation with its x −\rightarrow −x counterpart. Integrating this equation from 0 to x gives
(17.275) (Boltzmass distribution) which is just the Boltzmann distribution. 12M. San Miguel and J. M. Sancho, ‘‘A Colored-Noise Approach to Brownian Motion in Position Space. Corrections to the Smoluchowski Equation,’’ Journal of Statistical Physics 22, 605 (1980) (doi: 10.1007/BF01011341), Eq. (2.39).
Chapter 17. Stochastic Processes 17.6.3.1 Momentum Distribution Given the reduced distribution f(x, t) to the Smoluchowski equation (17.272), it is useful to see how to reconstruct the information about momentum that is implicit in the solution. First, as we already mentioned, the adiabatic relation (17.269) yields directly the mean momentum at each position:
p(x)
\gamma . (17.276) Here, to be precise, the ensemble average is over all dW, but at fixed position x. The mean position corresponding to the entire ensemble of particles must be averaged over the distribution f(x, t):
p(t)
x \approx Z \infty −\infty
\gamma . (17.277) (ensemble-mean momentum) Here, the subscript on the ensemble average emphasizes the average over the spatial distribution. However, the same approach applied to the second moment \langle \langle p2\rangle \rangle is problematic, because the square of
associated with the white-noise force dW/dt. This is a sign that we need to treat the adiabatic approximation with more care. To do this we return to the momentum equation of motion in Eqs. (17.240), dp = F(x, p; t) −\gammap
(17.278) and note that the same procedure that we used for the Ornstein–Uhlenbeck process in Section 17.3.6 applies here. Thus, we may rewrite the solution to Eq. (17.278) as
Z t dt′ e−\gamma(t−t′)
dt′ . (17.279) The essential point here is that the exponential convolution kernel here smooths the white-noise force, endowing it with finite power and removing the divergence in the second momentum moment. What we did before was to assume that the bracketed factor in the integrand changed slowly over the time scale \gamma−1, and that t ≫\gamma−1, so that we could make the replacement e−\gamma(t−t′)\Theta(t −t′) −\rightarrow 1
(17.280) Making this replacement and carrying out the integral leads directly to the previous adiabatic relation (17.269). However, we are only justified in doing this for the external-force term; the white-noise term is not constant over any time scale. Thus, we may write Eq. (17.279) as
Z t e−\gamma(t−t′) dW(t′). (17.281) Under the assumption that F, \gamma, and \sigma vary slowly over the damping time \gamma−1, we may regard these as constants. Coarse graining over the damping time scale, we can take the long-time limit of the above equation to remove transient behavior:
Z t e−\gamma(t−t′) dW(t′). (17.282) Now we regard the momentum as a function of the particle’s position, which is approximately constant. Under these conditions, the first term is the same mean position that we found before in Eq. (17.276), while the second represents (small) Gaussian fluctuations about the mean. Thus, the momentum distribution is a Gaussian tightly localized about the mean, and all that remains is to characterize the second moment. For
may adapt Eq. (17.154) to give the variance
2EE = \sigma2
(17.283)
17.7 Stochastic Boundary-Value Problems: Brownian Bridges where the latter expression follows from Eq. (17.256), and in the adiabatic regime can be regarded as correct even with spatial variation in temperature. Now when measuring the momentum distribution of the entire ensemble, we should average the momen- tum variance over the spatial distribution. However, it is more appropriate to average the second moment of momentum rather than the variance, because local variations in the mean momentum are reflected as variance in the global distribution. Thus, we require
p2(x)
\gamma2 . (17.284) Then averaging this over the spatial distribution from the Smoluchowski equation (17.272), we find
p2(t)
x = Z \infty −\infty dx f(x, t) \sigma2
\gamma2 . (distribution-average of second momentum moment) (17.285) At this point one can subtract the square of the distribution-averaged mean (17.277) to obtain the momentum variance of the whole distribution. However, Eq. (17.284) is already sufficient to characterize the momentum distribution in important cases such as the steady-state distribution in a potential well, where the mean momentum vanishes. 17.7 Stochastic Boundary-Value Problems: Brownian Bridges We have already studied the formalism to handle the simple stochastic differential equation dB = dW, (17.286)
in either It¯o or Stratonovich calculus, since the noise is additive.) However, suppose we add the additional
more difficult, as W(t) in general tends to wander away from zero. However, a (vanishingly) small subset of solutions obey this final condition, so in principle we could simulate many possible realizations of W(t), and discard them until we find one that returns sufficiently close to zero at t = 1. This kind of constrained random walk, or ‘‘stochastic loop’’, comes up, for example, in quantum field theory.13 This problem is also a nice example, showing alternate approaches to solving stochastic equations, and providing more insight into regular diffusion W(t). One simple guess at a solution is simply to force a regular Wiener path W(t) back to its initial point
(17.287) (Brownian bridge) This is called a Brownian bridge, and somewhat surprisingly, this solution satisfies the conditions above for our constrained Wiener path (with some cautions). This is something like viewing a Wiener path W(t) as composed of a linear drift to a final destinations plus fluctuations about zero, and then subtracting off the drift. To see that this is the case, first we note that B(t) is still a Gaussian random variable, since it is a
proper Wiener process. Dividing the unit time interval into N increments, with time steps ∆t = 1/N, with
(17.288) 13Holger Gies, Kurt Langfeld, and Laurent Moyaerts, ‘‘Casimir effect on the worldline,’’ Journal of High Energy Physics 06,
Stecker's Quantum Optics Notes — Page 941 (page-0738)¶
[D] Energy of the field in a cavity¶
Let \(a = (\hat{b}-\tilde{c})\), let \(\omega_{\bf k}=\Omega\). Then, since there are no counterpropagating fields with frequency \(\pm \hbar\omega\) in any quantum circuit containing only single-mode systems. The Hamiltonian takes the form $$ H = -i(\hat{x}_{in},a^\dagger)\frac{\partial}{\partial t} a $$
[D] Vacuum field fluctuations and radiation pressure noise¶
The Heisenberg uncertainty relation for \([\hat{p}], [\hat{x}]\) implies that they cannot be simultaneously measured to arbitrary precision. When these observables are defined as canonical operators with mass \(M\), frequency \(\Omega = \omega_0\). Then in the limit of zero temperature, one finds $$ \langle n(\infty) - 1\rangle / 2 $$
The field energy density is dominated by vacuum fluctuations arising from spontaneous emission processes. These give rise to radiation pressure noise proportional to \(|\frac{dE}{dz}|^2\) which couples back into the mechanical resonator dynamics via the force operator \(\hat{x}(t)\) satisfying $$ \ddot{\tilde{x}} + 2\gamma_{c}\dot{\tilde{x}}/\hbar = -i(\alpha_0,\Omega)a^\dagger $$
Summary of results for page-NNN and related sections¶
The above derivation shows that the field energy is dominated by vacuum fluctuations when \(k > \omega\). In particular, we found: $$ H(-\hat{x}_r) + i k = 2|\psi|^2 e^{-iE} $$
This leads to a non-trivial relation between the input and output quadratures. The result confirms that for large detuning \(\Delta\), one has good transmission through the cavity, while near resonance conditions are met only when \(T_1 \ll T_3\). In this regime $$ \frac{d}{dt}\langle n\rangle_{\infty} \approx -i(\hat{x}_c,\alpha^\top)\dot{\psi}(t) $$
The final expression demonstrates that the field energy scales as \(\Omega^2\) times a factor proportional to \(T_3/T_1\), confirming our earlier analysis.
17.7 Stochastic Boundary-Value Problems: Brownian Bridges Finally, we can examine the fluctuations of the (loop-style) Brownian bridge,
(17.294) so that we find
(17.295) (variance of Brownian bridge) Thus, the bridge fluctuates most when it is farthest away from either fixed endpoint, which is sensible. Again, to get a better idea of what these look like, 5 and 200 Brownian bridges are respectively shown
t Bo(t) -1 t Bo(t) -1 Both the symmetry about t = 0.5 and the character of diffusing away and returning are readily apparent here.
17.7.1 Finite Bridge Generation: Homogeneous Case¶
Chapter 17. Stochastic Processes 17.7.1 Finite Bridge Generation: Homogeneous Case Here we will consider generating a finite numerical approximation15 to a closed Wiener path in a more direct way than before. If we again divide the path into N increments, then we have time steps ∆t = 1 N , (17.296)
distributed random variable of zero mean and variance ∆t = 1/N. Thus we have the multidimensional probability density P(B1, . . . , BN−1) ∝exp −N N X j=1 (Bj −Bj−1)2 , (17.297) where by construction B0 = 0 and BN = 0 are not dependent variables. We will proceed by a coordinate transformation, changing variables to obtain an standard normal Gaussian distribution in every dimension. First, consider the sum in the exponent, which we may write as N X j=1
N−1 X j=2 (Bj −Bj−1)2 + B 2 1 + B 2 N−1 = 2B 2 1 −2B1B2 + 2B 2 2 −2B2B3 + 2B 2
N−2 −2BN−2BN−1 + 2B 2 N−1. (17.298) Now separating out the B1 dependence of the exponent, N X j=1
B1 −1 2B2 2 + 3 2B 2 2 −2B2B3 + 2B 2
N−2 −2BN−2BN−1 + 2B 2 N−1 = B′2 1 + 3 2B 2 2 −2B2B3 + 2B 2
N−2 −2BN−2BN−1 + 2B 2 N−1. (17.299) where we completed the square of B1, and we defined the transformed coordinate B′ 1 := \sqrt B1 −1 2B2 , (17.300) which encompasses all the dependence on B1, and enters in the exponent to give a normally distributed random variable with zero mean and variance 1/N. We now continue to complete squares and factor out the dependence on B2, B3, and so on. At the nth step, we have N X j=1
n−1 + cnB 2
N−2 −2BN−2BN−1 + 2B 2 N−1 = B′2
n−1 + cn Bn −1 cn Bn+1 2 + 2 −1 cn B 2
N−2 −2BN−2BN−1 + 2B 2 N−1 = B′2
N−2 −2BN−2BN−1 + 2B 2 N−1 (17.301) where in the basis step above we began with c1 = 2, and we have defined
2 −1 cn (17.302) 15Gies et al., op. cit.
17.7 Stochastic Boundary-Value Problems: Brownian Bridges and B′
Bn −1 cn Bn+1 . (17.303) This has the same form as the (n + 1)th step, which inductively completes all the squares. With these variables, the probability distribution is P(B′ 1, . . . , B′ N−1) ∝exp −N N−1 X j=1 B′2 j , (17.304) so that again the B′ j are independent, Gaussian random numbers of zero mean and variance 1/N. These can be chosen independently, and Eq. (17.303) can be solved to give Bn in terms of B′ n and Bn+1: Bn = B′ n \sqrtcn
cn . (17.305) Thus, the bridge coordinates should be generated in a backwards recurrence, given the forward recurrence for the coefficients cn. This is even more conveniently given in terms of standard-normal deviates zj, which can replace the B′ j: Bn = zn \sqrtNcn
cn . (17.306) Note also that the recurrence (17.302) has the solution
n . (17.307) These relations give the bridge directly in terms of easily generated deviates zn. Note that in principle we must also involve the Jacobian determinant of the coordinate transformation from Bj to B′ j. Effectively, we have taken the exponent (17.298), which we can write as a quadratic form as N X j=1
(17.308) where (Aab) is a square, tridiagonal matrix of dimension N −1 with a twos along the diagonal and ones for ev- ery other nonzero element. Since the matrix is symmetric, it is diagonalized by an orthogonal transformation
aDabB′ b, where (Dab) is diagonal (and in fact has only one nonzero eigenvalue, whose value is N −1), and B′ a := PabBb. We have effectively performed this diagonalization in constructing the above recurrence relations. Notice then that the Jacobian determinant is just the determinant of (Pab), which is just a constant factor that only affects the normalization of the probability density. Thus, we have justified our coordinate transformation to a new Gaussian distribution. To summarize, the algorithm to generate a Brownian bridge of N steps (i.e., N + 1 points B0, . . . BN,
- Generate the positions Bn for n = N, . . . , 1, according to the backwards recurrence BN = 0 Bn = zn r n
n n + 1 Bn+1, n = N −1, . . . , 1 B0 = 0. (17.309)
17.7.2 Finite Bridge Generation: Inhomogeneous Case¶
Chapter 17. Stochastic Processes We can similarly write a forward recurrence B0 = 0 Bn = zn rcn N + cnBn−1, n = 1, . . . , N −1 BN = 0, (17.310) where to simplify the notation, we have defined the recurrence coefficient cn := N −n N −n + 1. (17.311) This forward scheme gives finitely sampled Brownian bridges with the same statistics. This algorithm gives a simulated (finite) realization of a closed, stochastic Wiener path in terms of easily generated standard-normal random deviates. Note that this algorithm is easily generalized to generate samples of a Brownian bridge Ba\rightarrow b(t) running from a to b in unit time: the Bn in Eqs. (17.310) should be thought of as Bn −b (i.e., measured in terms of their distance to the ‘‘target’’ b), so that B0 = a Bn = zn rcn N + cn(Bn−1 −b) + b, n = 1, . . . , N −1 BN = b, (17.312) will act as a recurrence for this ‘‘open’’ bridge, with the cn as before. In the case of a more general Brownian bridge Bt(a\rightarrow b)(t′) running from a to b in time t, this recurrence is further generalized by changing the first term on the right-hand side of the second equation of Eqs. (17.312) from zn p cn/N to zn \sqrtcn∆t, where
17.7.2 Finite Bridge Generation: Inhomogeneous Case A slightly more complicated variation on the above recurrences arises when we allow for time-dependent drift and diffusion rates, according to
(17.313)
of finite differences, but strictly speaking this shouldn’t matter since the noise is still additive. In finite form, we have the multivariate probability density P(y1, . . . , yN−1) ∝exp −N N X j=1
\sigma 2 j−1 , (17.314) where again by construction y0 = 0 and yN = 0 are not dependent variables. We now thus have the exponent sum N X j=1
\sigma 2 j−1
\sigma 2
\sigma 2 + N X j=3
\sigma 2 j−1 = ¯y 2 \sigma 2
\sigma 2 + N X j=3
\sigma 2 j−1 , (17.315) where we have defined
N (17.316)
17.7 Stochastic Boundary-Value Problems: Brownian Bridges in order to begin eliminating the mean drifts. Continuing in this process, we have N X j=1
\sigma 2 j−1 = ¯y1 \sigma 2 + (¯y2 −¯y1)2 \sigma 2
\sigma 2 N−1 , (17.317) where ¯yn := yn −1 N n−1 X j=0 \alphaj, n \in 1, . . . , N, (17.318) again remembering yN = 0. Completing the first square as before, N X j=1
\sigma 2 j−1 = 1 \sigma 2 + 1 \sigma 2 ¯y 2 1 −2¯y1¯y2 \sigma 2 + ¯y 2 \sigma 2 + (¯y3 −¯y2)2 \sigma 2
\sigma 2 N−1 = \sigma 2
\sigma 2 0 \sigma 2 ¯y1 − \sigma 2 \sigma 2
¯y2 2 + 1 \sigma 2 − \sigma 2 \sigma 2 1 (\sigma 2
1 ) ¯y 2 + (¯y3 −¯y2)2 \sigma 2
\sigma 2 N−1 = y′2 1 + 1 \sigma 2 − \sigma 2 \sigma 2 1 (\sigma 2
1 ) ¯y 2 2 + (¯y3 −¯y2)2 \sigma 2
\sigma 2 N−1 , (17.319) where we have defined y′ 1 := s\sigma 2
\sigma 2 0 \sigma 2 ¯y1 − \sigma 2 \sigma 2
¯y2 . (17.320) At the nth stage of completing the square, we must handle terms of the form cn¯y 2 n −2¯yn¯yn+1 \sigma 2 n + 1 \sigma 2 n + \sigma 2 n+1 ¯y 2
¯yn − cn\sigma 2 n ¯yn+1 2 + 1 \sigma 2 n + \sigma 2 n+1 − cn\sigma 4 n ¯y 2 n+1 = y′2
n+1, (17.321) where we have defined the decoupled square y′
¯yn − cn\sigma 2 n ¯yn+1 (17.322) and the recursion
\sigma 2 n + \sigma 2 n+1 − cn\sigma 4 n , (17.323) thus inductively completing all the squares. Again, the y′ n are Gaussian numbers, such that we may solve to find the shifted positions ¯yn = y′ n \sqrtcn + cn\sigma 2 n ¯yn+1, (17.324) or in terms of standard-normal deviates, ¯yn = zn \sqrtNcn + cn\sigma 2 n ¯yn+1. (17.325) Then solving Eq. (17.318),
N n−1 X j=0 \alphaj (17.326) we find the actual bridge positions. To summarize, the algorithm is, to generate a stochastic, constrained path of N points y1, . . . yN, where
17.7.3 SDE and Integral Representations of the Brownian Bridge¶
Chapter 17. Stochastic Processes 1. Begin with means \alphan and standard-deviations \sigman for n \in 0, . . . , N = 1. 2. Generate the coefficients cn for n = 1, . . . , N −1, according to the recurrence c1 = 1 \sigma 2 + 1 \sigma 2 ,
\sigma 2 n + \sigma 2 n+1 − cn\sigma 4 n . (17.327) If many paths are to be generated, these coefficients only need to be generated once.
- Generate the shifted positions ¯yn for n = N, . . . , 1, according to the backwards recurrence ¯yN = − N−1 X j=0 \alphaj ¯yn = zn \sqrtNcn
cn\sigma2n . (17.328) 5. Generate the path positions yn for n = 1, . . . , N, using
N n−1 X j=0 \alphaj. (17.329) This gives a simulated realization of a closed, stochastic path with nonuniform drift and diffusion. Note that
in the previous section, owing to the conditioning. 17.7.3 SDE and Integral Representations of the Brownian Bridge The definition (17.287) of the Brownian bridge involves the pro-rated subtraction of the global drift of a Wiener path. However, given the completed-square construction of Eq. (17.7.1), we can also derive a local representation of the Brownian bridge as the solution of an SDE. Recall the backward recurrence, Eq. (17.309), which we may solve for Bn+1:
n + 1 n Bn −zn r n + 1 nN . (17.330) We can further rewrite this as
1 n Bn −∆Wn r n + 1 n , (17.331) where ∆Wn = zn/ \sqrt N, and then
N n Bn∆t −∆Wn r n + 1 n , (17.332)
and we note that (n+1)/n is only different from unity in a vanishingly small interval of small t, a correction we will ignore in the continuum limit: dB = B t dt −dW. (17.333)
integrated backwards in time. We can fix this by letting t −\rightarrow 1 −t, so that dt −\rightarrow −dt and dW −\rightarrow −dW, and thus dB = − B 1 −t dt + dW. (17.334) (SDE form of Brownian bridge)
17.7.4 State-Dependent Diffusion in Brownian Bridges¶
17.7 Stochastic Boundary-Value Problems: Brownian Bridges Thus, we have a representation of a Brownian bridge as an SDE solution. Here, the SDE is similar to an Ornstein–Uhlenbeck process, where the damping rate increases with time, diverging at t = 1 to ensure the return of the bridge. The solution to Eq. (17.334) is (see Problem 17.7)
Z t dW(t′) 1 −t′ , (17.335) (integral form of Brownian bridge) which can be verified by differentiating this expression. According to this definition, the Brownian bridge is clearly a Gaussian process, and the correlation function computed from this definition
(correlation function of Brownian bridge) (17.336) matches the same correlation function as computed from the first definition (17.287) (see Problem 17.8). This is sufficient to characterize the Gaussian process, and thus the two definitions are equivalent (at least in the sense of ensemble averages, not in the sense of individual paths). 17.7.4 State-Dependent Diffusion in Brownian Bridges Suppose we consider an inhomogeneous bridge problem that is slightly different from the inhomogeneous equation (17.313):
(state-dependent boundary-value SDE) (17.337) This is more complicated than our previous problem, since the drift and diffusion coefficients depend on the state of the system. While this equation is straightforward to integrate numerically, it is not at all straight- forward to do this integration subject to a bridge condition. This problem could be handled iteratively. For example, start by generating a solution with zero drift and constant diffusion, use this fiducial trajectory to generate the drift and diffusion coefficients, which then generates a new path; continue this process until con- vergence is reached. That is, if convergence is reached. In principle, this is a many-dimensional root-finding procedure, but a more stable method such as Newton iteration can be numerically cumbersome. We will, however, treat in more detail the somewhat simpler problem of state-dependent Stratonovich diffusion,
(17.338) (bridge with state-dependent diffusion)
to the initial-value SDE
(17.339) (solution for state-dependent diffusion)
has the same local statistics (i.e., statistics of increments) as the Wiener process dW(t), we need only verify the closure of the path. Now we turn to the closure of the solution to the SDE (17.339). First, we rewrite this SDE as
(17.340) still emphasizing the Stratonovich nature of this SDE. Now for some function S(y), where y(t) is defined by the SDE (17.339), the Stratonovich chain rule (17.200) gives
(17.341) where we treat dB(t) equivalently to dW(t) as far as the chain rule is concerned [dB(t) being a particular realization of dW(t)]. Suppose we take
\sigma(y), (17.342)
Chapter 17. Stochastic Processes or that is, S(y) is the antiderivative of 1/\sigma(y). Then Eq. (17.340) becomes
(17.343) Integrating both sides from t = 0 to 1 gives
(17.344) where we have again used the Stratonovich chain rule, so that the left-hand side was just the integral of dS[y(t)]. Thus,
(17.345) which implies
(17.346) as desired, so long as S(y) is an invertible function. Since it is the antiderivative of 1/\sigma(y), this is guaranteed if \sigma(y) is everywhere positive and finite, which is a reasonable restriction on the form of the SDE (17.338). However, the Stratonovich nature of this SDE is crucial: an It¯o equation naïvely of the same form has a drift term when converted to a Stratonovich equation, and our solution here does not cover that case. In fact, the integration procedure applied to arbitrary time t instead of t = 1 gives
(17.347)
Z y y(0) dy′ \sigma(y′). (17.348) Then the explicit solution in terms of a Brownian bridge is
(explicit solution for state-dependent diffusion bridge) (17.349) where as we have already stated, under reasonable assumptions S−1 exists. (For an example solution where \sigma(y) is a step function, see Problem 17.12.)Note that the solutions here can be generalized to the case of a return at time t = T instead of t = 1 by the replacements B(t) −\rightarrow \sqrt T B(t/T), dB(t) −\rightarrow \sqrt T dB(t/T) (17.350) in the above solutions. 17.7.4.1 Drift Returning to the more general SDE (17.337), we see that the presence of a drift term \alpha(y, t) dt is more difficult because even when integrating with respect to a bridge, closure of the solution is no longer guaranteed: the drift will in general make the solution ‘‘miss’’ the initial point, because in the formal bridge solution (interpreting the SDE as a Stratonovich equation),
Z t
Z t \sigma(y, t′) ◦dB(t′), (17.351) the first term is not guaranteed to vanish, and in fact the second term is also no longer guaranteed to vanish due to the influence of the first term on y(t). Additionally, any explicit time dependence in \sigma(y, t) in general causes the solution (17.339) in the drift-free case to fail. As we already mentioned, an iterative procedure to generate a closed solution may work for this, but in practice iteration tends not to converge well if the SDE coefficients change rapidly with y. An alternate approach is analogous to ‘‘shooting’’ methods for ODEs with boundary conditions. The idea here is to note that while the solution (17.351) may not close when integrated with respect to a a bridge
17.7 Stochastic Boundary-Value Problems: Brownian Bridges B(t), it will close for some bridge B0\rightarrow c(t) connecting 0 to c, as defined by Eq. (17.292). The extra drift compensates for any drift introduced by the SDE coefficients, and the existence of a value c that closes the path is guaranteed provided that the SDE is ‘‘reasonable enough’’ to guarantee continuous solutions (and solutions that vary continuously as the input noise varies). In general, this closing value c must be found numerically, via a root-finding algorithm. A concern with this method is the measure with which we generate the solutions. When we perform the analogous procedure with additive noise by closing a Wiener path to form a bridge as in Eq. (17.287), there is no problem: each Wiener path is associated with a unique bridge (with many paths associated with a single bridge), and choosing Wiener paths with the usual measure results in bridges being chosen with the correct measure (all Brownian bridges being equally likely). However, the nonlinear transformation in solving Eq. (17.337), as well as the unique association of a bridge B0\rightarrow c(t) with each solution y(t), where c is different for each solution, makes things more complicated. In particular, when generating Wiener paths, the relative probability of generating a bridge from 0 to c is the usual probability density for a Wiener path to end up at c after time t = 1:
\sqrt
(17.352) Therefore, if we require a bridge B0\rightarrow c\alpha(t) to generate a particular closed solution y\alpha(t), the relative proba- bility for this trajectory to occur is given by
\sqrt 2\pi e−c 2
(17.353) Thus, to generate an ensemble of solutions of Eq. (17.337), each generated solution should only be kept with a relative probability w\alpha (e.g., by the rejection method). Alternately, when computing an ensemble average with respect to the solutions y\alpha(t), the average should be computed as a weighted average, where the weight of each member of the ensemble is w\alpha. This procedure is valid for It¯o or Stratonovich equations, provided that appropriate integration methods are used in each case. Finally, note that this measure business is not
probability. One final possibility for an algorithm in the case where the SDE coefficients are time-independent [i.e.,
then the two paths can be spliced together to realize a bridge (i.e., the solution is ya(t) for t < \tau, and yb(t) thereafter). If there is no such crossing, the paths are rejected and the process is repeated until successful. 17.7.4.2 Lamperti Transform The above considerations of rescaling stochastic processes can be elegantly viewed in the framework of the Lamperti transform,17 which is a transformation that equalizes the variances within a stochastic process. This idea applies to SDEs as follows. Consider the It¯o SDE
(17.354) and the transformation
(17.355) Then the It¯o chain rule (17.193) gives dz =
2S′′(y) \beta2
(17.356) 16Stefano M. Iacus, Simulation and Inference for Stochastic Differential Equations (Springer, 2008), Section 2.13. 17after J. W. Lamperti, ‘‘Semi-stable stochastic processes,’’ Transactions of the American Mathematical Society 104, 62 (1962) (doi: 10.1090/S0002-9947-1962-0138128-7). See Stefano M. Iacus, op. cit., Section 1.11.4.
Chapter 17. Stochastic Processes and the Lamperti transformation obtains by choosing
\beta(y), (17.357) or
Z y y0 dy′ \beta(y′). (17.358) (Lamperti transform)
dz = \alpha
dt + dW, (17.359) (Lamperti-transformed process) which now is driven by additive noise. The complexity of the multiplicative noise in the dy SDE is thus moved from the stochastic to the deterministic term, and reconstructing the dynamics in z requires inverting
same sign over the domain of y). Recall from the It¯o–Stratonovich conversion (17.184) that a Stratonovich diffusion equation
(17.360)
S—something not in general possible if the drift term remains in Eq. (17.354). 17.7.4.3 Temporal Rescaling A similar transformation, closer to the original transformation introduced by Lamperti to study self-similar processes18 is a trajectory-dependent temporal rescaling, with a transformation from time t to t′, such that
(17.361) (temporal rescaling) In the It¯o SDE (17.354), this transformation leads to the equivalent SDE
(17.362) (rescaled It¯o SDE) Note that, as in Eq. (17.359), the variable diffusion rate disappears, and we are left with additive noise. However, the drift term is different, and in particular there is not a second-order It¯o correction term. On the other hand, the Stratonovich SDE (17.360) transforms simply to
(17.363) (rescaled Stratonovich SDE) with no drift, just as in the case of the Lamperti transform. As in our discussion in Section 17.7.4, the interpretation here is simple if we consider Brownian bridges. In the It¯o case, if the noise dW(t) is a standard Brownian bridge, then the closure of the bridge is preserved under the temporal rescaling if and only if \alpha = 0. On the other hand, the Stratonovich SDE always preserves the closure of the path for any temporal rescaling. Note, however, that in the It¯o case, even an It¯o equation that is equivalent to the
18J. W. Lamperti, op. cit.; see also Krzysztof Burnecki, Makoto Maejima, and Aleksander Weron ‘‘The Lamperti transfor- mation for self-similar processes,’’ Yokohama Mathematical Journal 44, 25 (1997).
17.8.1 Wiener Process¶
17.8 Boundary Crossings 17.8 Boundary Crossings¶
what is the probability that W(t) will cross a ‘‘boundary’’ at d after some time t? This is a useful analysis, for example, in financial mathematics, and a wide range of other areas.19 17.8.1 Wiener Process
threshold d > 0 over this time interval from 0 to t? One approach to this problem is via the Reflection Principle.20 The basic idea is illustrated in the plot below of a Wiener path. Once the path hits the threshold at d, it continues to diffuse on. However, note that we can construct an equally likely path by mirroring the path after the crossing about d, as shown in the light, mirrored trace. This is true when the probability distribution for the increments is symmetric, which is of course the case for standard Brownian motion. Note that we are also implicitly assuming that there is a well-defined crossing, that is, given that
process. t Wo(t) d T Now let’s apply this to the crossing probability, which we will write
(17.364) where \taud is the time of the first crossing of W(t) through d:
(17.365) Then we can partition the probability according to whether W(t) ends up below or above the boundary:
(17.366)
(17.367) 19David Siegmund, ‘‘Boundary crossing probabilities and statistical applications,’’ The Annals of Statistics 14, 361 (1986)
20A particularly readable reference for the Reflection Principle appears in Joseph T. Chang, Stochastic Processes, available at http://www.stat.yale.edu/~jtc5/251/stochastic-processes.pdf. We follow his arguments here for single-boundary crossing probabilities for Wiener paths and bridges. For the Reflection Principle for Wiener paths, see also Kurt Jacobs, Stochastic Processes for Physicists: Understanding Noisy Systems (Cambridge, 2010).
17.8.2 Standard Brownian Bridge¶
Chapter 17. Stochastic Processes According to the Reflection Principle,
- (17.368) By construction in the figure, for each path that has crossed through t and ends up with W(t) < d, there is an equally likely (mirrored) path with W(t) > d. So the probability of ending up above or below d is the same, given that the boundary is crossed. Then we have
(17.369) and solving for Pcross(d, t),
(17.370) Now W(t) is normally distributed with variance t, so we have
Z \infty d dW \sqrt
(17.371) or finally21
d \sqrt 2t (crossing probability of W(t) through d in time t) (17.372) for the crossing probability past d in time t. Note that this probability converges to unity as t −\rightarrow \infty; even though sample paths can start off in the negative direction, given enough time, they will tend to return to the origin and go far enough into the positive direction to cross the boundary anyway. For small t, this probability reduces to
r 2t \pid2 , (17.373) so that the crossing probability is exponentially suppressed as t −\rightarrow 0. The result (17.372) can also be interpreted as a cumulative probability distribution for the crossing to occur before time t. Then the probability density for the first-passage time \taud is given by22
d \sqrt 2t t=x = d \sqrt
(probability density for first-passage time) (17.374) where we have simply differentiated the probability (17.372). Note that the probability density is suppressed exponentially at short times, as e−d2/2\taud, but decays only as \tau −3/2 d at large times. Thus, this distribution has no mean or variance, meaning that may well take a long time to cross a boundary, although it will almost surely happen. However, the most likely value of \taud (i.e., that maximizes the probability density) is d2/3. 17.8.2 Standard Brownian Bridge A similar argument works to compute boundary-crossing probabilities for the standard Brownian bridge B(t). Here, however, the Reflection Principle is slightly different. The path crossing the boundary at d must return to zero. The equivalent path that is mirrored after the crossing then must return to 2d instead of 0, as shown below. 21cf. Andrei N. Borodin and Paavo Salminen, Handbook of Brownian Motion—Facts and Formulae, 2nd ed. (Birkhäuser, 2002), p. 153, formula 1.1.4. 22cf. Andrei N. Borodin and Paavo Salminen, op. cit., p. 198, formula 2.0.2.
17.8 Boundary Crossings t Bo(t) d 2d Here we will represent the bridge B(t) as an ordinary Wiener path W(t), but subject to the condition
(17.375)
. (17.376) To compute the numerator,
(17.377) and using the Reflection Principle,
(17.378) Thus, the crossing probability becomes
(17.379)
between x and x + dx; that is, we are referring to probability densities here, so that we are taking the ratios
\sqrt 2\pi, then we have23
(17.380) (crossing probability of B(t) through d) for the crossing probability of the bridge. Note that the crossing probability is the same as the probability for the peak of the bridge to cross d:
(17.381) 23cf. Andrei N. Borodin and Paavo Salminen, op. cit., p. 154, formula 1.2.8.
17.8.4 First-Passage Time of the Brownian Bridge¶
Chapter 17. Stochastic Processes Thus, the probability density function for the maximum of B(t) is
(17.382) or24
(probability density for Brownian-bridge maximum) (17.383) We can also compute the moments of this distribution, via DD [sup{B(t)}]nEE = Z \infty
Z \infty 4xn+1e−2x2, (17.384) with the result DD [sup{B(t)}]nEE
1 + n , (moments for Brownian-bridge maximum) (17.385) which gives (1/4) p
p \pi/2, and 1/2, for n = 1, 2, 3, and 4, respectively. Note that the results
since the pinning essentially rescales the size of the bridge by \sqrt T, which is equivalent to scaling any distances by 1/ \sqrt T. Thus, for example, the crossing probability (17.380) is given by letting d2 −\rightarrow d2/T, and the peak density (17.383) is given by letting x2 −\rightarrow x2/T and x dx −\rightarrow x dx/T in fsup{B(t)}(x) dx. 17.8.3 Brownian Bridge The treatment above for the standard Brownian bridge is easy to generalize to a bridge B0\rightarrow c(t) that connects
(17.386) such that Eq. (17.380) generalizes to
(crossing probability of B0\rightarrow c(t) through d) (17.387) Then Eq. (17.383) similarly becomes
(probability density for Brownian-bridge maximum) (17.388) These results can be generalized to a bridge pinned to c at time t by scaling d −\rightarrow d/ \sqrt t and c −\rightarrow c/ \sqrt t. 17.8.4 First-Passage Time of the Brownian Bridge In analyzing the crossing probability for the Wiener path, we were able to simply derive the result (17.374) for the density of the first-passage time \taud. It is only slightly more complicated to do this for a Brownian bridge, so we will carry out the derivation here. For the Brownian bridge, we can define the first-passage time as
(17.389)
up to time \tau (before the final time t) is
(17.390) 24cf. Eq. (2) in Jim Pitman and Marc Yor, ‘‘On the distribution of ranked heights of excursions of a Brownian bridge,’’ Annals of Probability 29, 361 (2001) (doi: 10.1214/aop/1008956334).
17.8 Boundary Crossings when written in terms of a Wiener path. We don’t know this crossing probability, but from Eq. (17.387) we do know the crossing time of the bridge up to the final running time t. To rewrite Eq. (17.390) in terms of this crossing probability, we can proceed as follows:
\sqrt
= \sqrt
Z \infty −\infty
(17.391) The second factor in the integrand is just a Gaussian probability density,
p
(17.392) and the first factor is given by Eq. (17.387) as
\sqrt
= \sqrt
= \sqrt
(17.393) for the case where the intermediate point z is not past the boundary d, and
\sqrt
= \sqrt
(17.394) for the case where the path has already crossed the boundary by time \tau. Putting these expressions into
crossing probability
2 erfc
p
! + 1 2 e−2d(d−c)/t erfc
p
! . (first-passage cumulative probability, Brownian bridge) (17.395) To check the normalization of this probability note that erfc(x) −\rightarrow 0 as x −\rightarrow \inftyand erfc(x) −\rightarrow 2 as x −\rightarrow −\infty. Then as \tau −\rightarrow t, Eq. (17.395) becomes
(17.396) which gives the correct crossing probability from Eq. (17.387). The probability density for the first-passage time is then simply given by differentiating with respect to \tau:25
\sqrt t p 2\pix3(t −x) d e−(d t−cx)2/2tx(t−x). (probability density for Brownian-bridge first-passage time) (17.397) Note that in the limit of large t, this reduces to the Wiener-process result (17.374), as it should. This expression is also unnormalized if c < d: in this case the normalization is given by Eq. (17.396), since the first-passage time is undefined in the event that a path does not cross the boundary. 25L. Beghin and E. Orsingher, ‘‘On the maximum of the generalized Brownian bridge,’’ Lithuanian Mathematical Journal 39, 157 (1999) (doi: 10.1007/BF02469280), Eq. (2.15).
17.9.1 Wiener Process¶
Chapter 17. Stochastic Processes The moments of the first-passage time can be written in terms of the density (17.397) as
\tau n d
= r t 2\pi d Z t dx xn−3/2 \sqrtt −x e−(d t−cx)2/2tx(t−x). (17.398) Changing variables to y = t/x −1 gives the alternate expression
\tau n d
\sqrt 2\pi Z \infty dy (1 + y)−n \sqrty
(17.399) For n = 0, this expression gives the correct normalization DD \tau 0 d EE
(17.400) The integral is easier to evaluate for n < 0 than n > 0. For example DD \tau −1 d EE
td2 e−d(d−c+|d−c|)/t (17.401) gives the first inverse moment—the well-defined value here and of the other inverse moments is an indication of how heavily the value \taud = 0 is suppressed.
is also the statistic for the last passage time, because of the time-reversal symmetry of the bridge. Moments of the last passage time can be calculated in the same way as the regular moments. Adapting the same variable change leading to Eq. (17.399), the post-first-passage-time moments are
(t −\taud)n
\sqrt 2\pi Z \infty
(17.402) Then, for example, the first inverse moment is given by
(t −\taud)−1
|d −c|t e−d(d−c+|d−c|)/t. (17.403)
to have the same form as the boundary-crossing probability e−2d2/t from Eq. (17.380), but with an extra factor of 2 (recall that this applies to the case where paths that do not cross the boundary count as zero in this ensemble average). 17.9 Escape Probability Similar to the boundary-crossing problem is the escape problem, which is concerned whether a stochastic process leaves an interval. There are implicitly two boundaries involved, and the idea is to see whether the process touches either boundary. [This is closely related to whether the process touches both boundaries,
P(A) and P(B) are given by the appropriate single-boundary-crossing probabilities.] 17.9.1 Wiener Process We will set up the problem as follows: a Wiener process W(t) begins between two barriers separated by distance L, and the starting point is a distance a from one of the barriers. For concreteness, we take the
define the interval [−a, L −a].
17.9 Escape Probability L a Of course, we would obtain the same answers by instead using the interval [a −L, a]. The question is, what is the probability to touch either boundary, and thus to escape the interval, in time t? We will approach this problem in the same way as the boundary-crossing, making use of the boundary- crossing probability (17.372) in the process.26 Actually we will first consider a slightly different problem, which is, what is the probability that W(t) touches the upper boundary first? That is, we only count the escapes where the first escape is through the upper boundary. To compute this, we will define some events (sets of outcomes). First, we define event U1 to be set of all outcomes where the process touches the upper boundary: U1 : We are illustrating the trajectory schematically here; the trajectory may be much more complicated and touch either boundary many more times than we have indicated. The conical ‘‘spray’’ of trajectories to the right indicates that we don’t particularly care what happens to the trajectory afterwards. Now we are only interested in the cases where the process touches the upper boundary first, but we have included cases where the process touches the lower boundary and then the upper boundary, since U1 includes any outcome that touches the upper boundary. We will denote this set L1: L1 : To properly count the events we want, we should delete events in L1. But not all of them! In L1 we included paths that touch the upper boundary before touching the lower boundary, and we want to count these. We will denote this set U2: U2 : But again, in this set, we are counting paths that touch the lower boundary before the indicated touchings, and we don’t want to count these. We will denote this set L2. L2 : Continuing in this way, we should define the set Uj to be the set of all paths that touch the upper boundary j times, with j −1 touchings of the lower boundary ‘‘interleaved,’’ and Lj to be the set of all paths that touch the lower boundary and then alternating between the upper boundary and lower boundary, with j touchings of each boundary. (Once a boundary touches a boundary, it is okay for it to touch it again before touching the other boundary in these definitions.) Thus, we have argued that the set of all paths that touch the upper boundary first is Aupper first = U1 −L1 + U2 −L2 + . . . (17.404) The probabilities of the events on the right-hand side are easy to compute, using the Reflection Principle. (Indeed, notice how the ‘‘reflections’’ pop up here, in a way analogous to the infinity of reflections in a Fabry– Perot cavity.) The idea as before is to ‘‘unwrap’’ trajectories via reflections, and the resulting probability 26The method here was applied to the Wiener process in J. L. Doob, ‘‘Heuristic Approach to the Kolmogorov–Smirnov Theorems,’’ The Annals of Mathematical Statistics 20, 393 (1949) (doi: 10.1214/aoms/1177729991); and T. W. Anderson, ‘‘A Modification of the Sequential Probability Ratio Test to Reduce the Sample Size,’’ The Annals of Mathematical Statistics 31, 165 (1960) (doi: 10.1214/aoms/1177705996).
17.9.2 Standard Brownian Bridge¶
Chapter 17. Stochastic Processes is just the single-boundary crossing probability (17.372), where the distance d is the total vertical distance traversed in each diagram (noting that the distances to the boundaries are L −a and a above and below the dashed line, respectively. Thus, for example,
L −a \sqrt 2t
L + a \sqrt 2t
3L −a \sqrt 2t
3L + a \sqrt 2t , (17.405) and so on. The probability to touch the upper boundary before the lower boundary is then
L −a \sqrt 2t −erfc L + a \sqrt 2t + erfc 3L −a \sqrt 2t −erfc 3L + a \sqrt 2t + . . . = \infty X j=1 erfc (2j −1)L −a \sqrt 2t −erfc (2j −1)L + a \sqrt 2t . (17.406) The probability to touch the lower boundary before the upper boundary is simply given by the replacement a −\rightarrow L −a in the above expression:
\infty X j=1 erfc (2j −2)L + a \sqrt 2t −erfc (2j)L −a \sqrt 2t = erfc a \sqrt 2t + \infty X j=1 erfc 2jL + a \sqrt 2t −erfc 2jL −a \sqrt 2t . (17.407) The escape probability is the sum of the above two probabilities,
(17.408) since they represent two disjoint sets of outcomes:27
a \sqrt 2t + \infty X j=1 (−1)j erfc jL + a \sqrt 2t −erfc jL −a \sqrt 2t = \infty X j=−\infty (−1)j sgn(j + 0+) erfc |a + jL| \sqrt 2t . (escape probability for Wiener process) (17.409) Note that the 0+ is included in the sgn function so that the j = 0 term is positive. 17.9.2 Standard Brownian Bridge The calculation for the Brownian bridge goes in essentially the same was as for the Wiener path.28 We set
the Wiener path B(t). 27cf. Andrei N. Borodin and Paavo Salminen, op. cit., p. 167, formula 1.7.4(2), and see p. 641 for the function definition. 28Bruno Casella and Gareth O. Roberts, ‘‘Exact Monte Carlo Simulation of Killed Diffusions,’’ Advances in Applied Probability 40, 273 (2008) (doi: 10.1239/aap/1208358896). (Note the typo in the expression for qj.) Also Alexandros Beskos, Stefano Peluchetti, and Gareth Roberts, ‘‘ϵ-Strong Simulation of the Brownian Path,’’ arXiv.org preprint (arXiv: 1110.0110v1).
17.9 Escape Probability L a Again, we consider the altered problem of: what is the probability that W(t) touches the upper boundary first? Again, we define events Uj and Lj as in the Wiener-path case, except that now the final points of the paths are pinned down: U1 : L1 : U2 : L2 : and so on. Continuing in this way, we again define the set Uj to be the set of all paths that touch the upper boundary j times, with j −1 touchings of the lower boundary ‘‘interleaved,’’ and Lj to be the set of all paths that touch the lower boundary and then alternating between the upper boundary and lower boundary, with j touchings of each boundary. Thus, the set of all paths that touch the upper boundary first is again Aupper first = U1 −L1 + U2 −L2 + . . . (17.410) The probabilities of the events on the right-hand side are easy to compute, using the Reflection Principle. The probability in each diagram is just the single-boundary (bridge) crossing probability (17.380), where the distance d is half the total vertical distance traversed in each diagram. Thus, for example,
(17.411) and so on. The probability to touch the upper boundary before the lower boundary is then
= \infty X j=1 h e−2(jL−a)2 −e−2(jL)2i . (17.412) The probability to touch the lower boundary before the upper boundary is again simply given by the replacement a −\rightarrow L −a in the above expression: Plower first = \infty X j=1 h e−2[(j−1)L+a)2 −e−2(jL)2i
\infty X j=1 h e−2(jL+a)2 −e−2(jL)2i (17.413) The escape probability is once again the sum of the above two probabilities, Pescape = Pupper first + Plower first. (17.414)
17.9.3 Brownian Bridge¶
Chapter 17. Stochastic Processes The resulting expression is Pescape = e−2a2 + \infty X j=1 h e−2(jL−a)2 + e−2(jL+a)2 −2e−2(jL)2i (0 < a < L)
\infty X j=−\infty h e−2(a+jL)2 −e−2(jL)2i . (escape probability, standard Brownian bridge) (17.415) Notice that besides the obvious difference in the functions appearing here compared to the Wiener case, the structure is different: an a-independent term (the last term here) does not appear in the Wiener case. 17.9.3 Brownian Bridge It is not difficult to generalize the above escape probability to a more general Brownian bridge B0\rightarrow c(t),
various distances d. According to (17.387) we just need to change these to terms of the form exp[−2d(d−c)], with
\infty X j=1 h e−2(jL−a)(jL−a−c) + e−2(jL+a)(jL+a−c) −2e−2(jL)(jL−c)i (0 < a < L; −a < c < L −a) (escape probability, Brownian bridge) (17.416) as the result. 17.10 Dirichlet Problem and Connection to Electrostatics 17.10.1 Laplace Equation Another interesting context in which boundary crossing and escape arises is in the Laplace equation, subject to Dirichlet boundary conditions. That is, suppose we want the solution to the Laplace equation
∀r\in \partial D. (17.417) (Dirichlet problem) on some bounded domain D, where \phi\partial (r) fixes the solution on the domain boundary \partial D. (That is, find the electrostatic potential in a charge-free region, given the potential/voltage on a bounding surface.) We will now show that the solution can be written as the path average
DD
EE
(17.418) (stochastic solution) where W(t) is a vector Wiener process, and \tau\partial is the first passage time of W(t) through \partial D. That is, we start a bunch of Wiener paths from r, let them go until they hit the boundary, and then take the average of the boundary values where the paths hit. Actually, this solution requires that the boundary be sufficiently nice, which is certainly true for physical boundaries in electrostatic problems. (More technically, the surface should satisfy the ‘‘Poincaré cone condition,’’ which basically says that at each point on the surface, you can attach a code of finite angle and length that doesn’t intersect the boundary except at the attachment point. This rules out, for example, interior boundaries of arbitrarily small area, or a sufficiently severe ‘‘kink’’ in the surface.29) That the stochastic solution (17.418) has the correct boundary values is reasonably obvious, because as the point r moves close to a point on the boundary, the paths starting from that point will hit the 29Kiyosi Itô and Henry P. McKean, Jr., Diffusion Processes and their Sample Paths (Springer, 1974), pp. 257, 261-4 (doi:
17.10 Dirichlet Problem and Connection to Electrostatics nearest surface with probability approaching unity. This can be seen, for example, from the first-passage- time density for the Wiener path in Eq. (17.374), where the peak of the density shifts to zero as d −\rightarrow 0, or that the crossing probability (17.372) in any finite time converges to unity as d −\rightarrow 0. However, the boundary function \phi\partial should be continuous, as the paths from r will always average over a small region of the boundary, even as the source point approaches the boundary. Now to show that the expression (17.418) satisfies the Laplace equation. Suppose we draw a sphere of radius R, centered at r, such that the sphere lies entirely within D, as shown below. R r W(t) D \partial D \tauR \tau\partial The idea is that the Wiener path must cross the sphere before it crosses the boundary. Then suppose we rewrite Eq. (17.418) as
DD
EE
W(\tauR) , (17.419) where \tauR is the first crossing time of the sphere, and
(17.420) We haven’t really done much here, except to split (pathwise) the time interval into pre- and post-\tauR, and we are explicitly taking the ensemble average after \tauR separately from the ensemble average over all possible first crossings of the sphere represented by \tauR. However, since the path ∆W(\tau\partial ) after touching the sphere acts itself like a Wiener path, we can use the solution (17.418) to replace the inner ensemble average:
DD
EE W(\tauR). (17.421) Now since the Wiener path is equally likely to have its first touching point W(\tauR) at any point on the sphere, we can simplify the notation a bit to write
DD
EE |R|=R, (17.422) where R determines some point on the sphere, and the ensemble average is simply a uniform average over the surface of the sphere. That \phi(r) is the average value of \phi on a sphere centered on r is a necessary and sufficient condition for \phi to be a harmonic function (i.e., a solution to the Laplace equation).30 17.10.2 Laplacian as Spherical Average We can make the above argument about the averaging property of harmonic functions more precise by working it out somewhat more explicitly, and in the process arrive at a useful representation of the Laplacian. First, consider the Taylor expansion
(17.423) 30See David J. Griffiths, Introduction to Electrodynamics, 4th ed. (Prentice Hall, 2013), Section 3.1.4, p. 117. The idea is to consider a point charge outside the sphere, and show that the averaging statement is true for this case. Then the same statement must be true for any collection of charges, and conversely for any solution of the Laplace equation, since any physical solution may be regarded as being produced by some charge distribution.
Chapter 17. Stochastic Processes where repeated indices are summed. Averaging over all orientations of R at fixed distance R gives DD
EE
DD \partial 2
R 2 \alpha EE
(17.424) where all first derivatives have vanished under the symmetric average. Now we can rewrite the last term using ** d X \alpha=1 \partial 2
R 2 \alpha ++ |R|=R = d X \alpha=1 \partial 2
DD R 2 \alpha EE |R|=R = d X \alpha=1 \partial 2
R2 d = R2
(17.425) where we are writing out the sum explicitly now over d dimensions, and we used the fact that R 2 \alpha is independent of \alpha once averaged over all orientations. Then putting this result into Eq. (17.424) and taking the limit R −\rightarrow 0
R\rightarrow 0+ 2d R2 DD
EE
, (Laplacian as spherical average) (17.426) which is a representation of the Laplacian operator in terms of an average over a small sphere around r. This equation immediately implies that Eq. (17.422) is equivalent to the Laplace equation (17.417), completing our proof of the Wiener-path solution. 17.10.3 Poisson Equation The same basic approach31 works for the Poisson equation
∀r\in \partial D, (17.427) (Dirichlet Poisson problem) which is the same as the original problem (17.417), but with the addition of a source \rho(r) (a factor of 1/ϵ0 is absorbed into the source for simplicity). The solution is the same as before, but with the addition of a source-averaging term:
**
++
. (stochastic Poisson solution) (17.428) The same argument leading to Eq. (17.419) applies, but the second term splits into two parts, before and after \tauR. The result is then
DD
EE |R|=R − ** Z \tauR
++ W(\tauR) . (17.429) The portion of the integral from \tauR to \tau\partial was already absorbed into the solution in the first term here. Rearranging and multiplying by 2d/R2, 2d R2 DD
EE
= d R2
++ W(\tauR) . (17.430) 31For a rigorous version of the argument here, see Sidney C. Port and Charles J. Stone, Brownian Motion and Classical Potential Theory (Academic Press, 1978) (ISBN: 0124335942), Proposition 5.2 on p. 14 and Section 4.6 on p. 114.
17.11 Feynman–Kac Formula Taking the limit R −\rightarrow 0, we can use Eq. (17.426) on the left-hand side, and on the right-hand side, we can treat \rho as being a constant with respect to the integral. Thus
R2
\tauR
W(\tauR). (17.431) At this point, we have dropped the limit R −\rightarrow 0 on the right-hand side, as it is not necessary. The statistic involved here is the mean of the first-passage time through a sphere of radius R in d dimensions. Now to calculate the remaining expectation value. First, recall that for a vector Wiener process in d dimensions, the Euclidean norm is given on average by
∥W(t)∥2 = td, (17.432) since there are d independent directions of displacement, each contributing t to the variance. Now what we would like to show is that this is still true if t is replaced by the stopping time \tauR, in the sense that
∥W(\tauR)∥2
(17.433) or indeed any other stopping time. First, let’s take t to be some very large time, such that almost certainly \tauR < t (i.e., 0 < \tauR < t). For some paths this will not be true, such that the argument below will miss them; in these cases we can take \tauR to be equal to t, which will produce an error that vanishes in the limit t −\rightarrow \infty. Then starting with Eq. (17.432), td =
∥W(t)∥2 =
=
∥W(t) −W(\tauR)∥2 +
∥W(\tauR)∥2 −
, (17.434) where the last term is zero, because the parts of the path W(t) before and after \tauR are independent, and at least the first factor has zero mean. Then taking on the first term on the right-hand side, we can think of all paths that start at a particular stopping point W(\tauR), such that the ensemble average of ∥W(t) −W(\tauR)∥ is just (t −\tauR)d. Then continue the ensemble average over all stopping times \tauR, so that
∥W(t) −W(\tauR)∥2 =
(t −\tauR)d
(17.435) Putting this together with Eq. (17.434), we find the desired result
∥W(\tauR)∥2
(17.436)
ensemble averages involving \tauR are correct. By definition of \tauR, the left-hand side is R2, so
d . (17.437) Putting this into Eq. (17.430), we see that it reduces to the Poisson equation (17.427) as desired. 17.11 Feynman–Kac Formula To introduce a very powerful method of computing expectation values for SDEs, we will introduce the Feynman–Kac formula,32 which solves diffusion problems in terms of integrals over solutions to SDEs. 32R. P. Feynman, ‘‘Space-Time Approach to Non-Relativistic Quantum Mechanics,’’ Reviews of Modern Physics 20, 367 (1948) (doi: 10.1103/RevModPhys.20.367); M. Kac, ‘‘On Distributions of Certain Wiener Functionals,’’ Transactions of the American Mathematical Society 65, 1 (1949) (doi: 10.1090/S0002-9947-1949-0027960-X); M. Kac, ‘‘On Some Connections between Probability Theory and Differential and Integral Equations,’’ in Proceedings of the Second Berkeley Symposium on Mathematical Statistics and Probability (University of California Press, 1951), p. 189 (http://projecteuclid.org/euclid. bsmsp/1200500229).
Chapter 17. Stochastic Processes There are numerous forms of this formula, and we will develop a relatively simple form33 that considers a forced diffusion equation for the distribution f(x, t),
2\partial 2 x f −V (x, t)f + g(x, t), (17.438) (PDE for Feynman–Kac formula) subject to the initial condition
(17.439) (initial condition for Feynman–Kac formula) The Feynman–Kac formula gives the solution of (17.438) for t > 0 as
** f0[x + W(t)] exp − Z t dt′ V [x + W(t′), t −t′] + Z t dt′ g[x + W(t′), t −t′] exp
− Z t′ dt′′ V [x + W(t′′), t −t′′]
, (Feynman–Kac formula) (17.440) where the ensemble average is over all realizations of W(t). Before continuing to prove this formula, notice that it reduces correctly to f0(x) as t −\rightarrow 0. Also, for a simple, undriven diffusion,
2\partial 2 x f, (17.441) the formula reduces to
DD f0[x + W(t)] EE = Z \infty −\infty dW f0(W) \sqrt
(17.442) which is the convolution of f0(x) with a Gaussian of variance t, as we expect. The extra terms V and g introduce damping (or ‘‘killing’’ of diffusing particles) and particle sources, respectively, that clearly make the solution more complicated. 17.11.1 Proof: Simple Diffusion But now to prove the Feynman–Kac formula (17.440), which we will do in increasingly complex stages. First, take again the simple case where V = g = 0, so that we have the simple diffusion equation,
2\partial 2 x f, (17.443) with solution
DD f0[x + W(t)] EE . (17.444) We will show explicitly that this is the solution by simply differentiating it with respect to time. Technically, we will just compute the differential df(t), regarding x as a fixed parameter,
DD \partial xf0[x + W(t)] dW EE + 1 DD \partial 2 x f0[x + W(t)] dt EE , (17.445) where we have differentiated according to the It¯o rule. The dW term vanishes in the ensemble average,
DD f ′′ 0 [x + W(t)] EE dt, (17.446) and differentiating (17.444) with respect to x to evaluate the averages on the right-hand side, we see that the first derivatives cancel, and the second-derivative term gives the remaining term we need in Eq. (17.443), after canceling factors of dt. 33cf. Andrei N. Borodin and Paavo Salminen, op. cit., p. 103.
17.11 Feynman–Kac Formula 17.11.2 Proof: Diffusion with Damping¶
generate is
2\partial 2 x f −V (x, t)f, (17.447) and we want to show that this PDE is solved by
** f0[x + W(t)] exp − Z t dt′ V [x + W(t′), t −t′]
. (17.448) The procedure here is somewhat more complicated than for the simple diffusion in Eqs. (17.443) and (17.444). To keep the calculation organized, consider the quantity34
− Z t dt′ V [x + W(t′), t′′ −t′] , (17.449) where we assume f(x, t) to satisfy the PDE (17.447), with x and t′′ effectively fixed parameters. This is something like the solution (17.448) with the initial condition replaced by the time-dependent solution. We will now compute the differential of M:
−\partial tf[x + W(t), t′′ −t] dt + \partial xf[x + W(t), t′′ −t] dW + 1 2\partial 2 x f[x + W(t), t′′ −t] dt e− R t 0 dt′ V [x+W (t′),t′′−t] + f[x + W(t), t′′ −t] V [x + W(t), t′′ −t] e− R t 0 dt′ V [x+W (t′),t′′−t] dt. (17.450) The first, third, and fourth terms here vanish together since we assumed f(x, t) to satisfy (17.447). Thus,
− Z t dt′ V [x + W(t′), t′′ −t′] dW. (17.451)
average behavior of our solution. We will make use of this as follows. Note that evaluating \langle \langle M(t)\rangle \rangle at t = t′′, we have DD M(t′′) EE = ** f[x + W(t′′), 0] exp
− Z t′′ dt′ V [x + W(t′), t′′ −t′]
, (17.452)
gives the general solution DD M(0) EE
(17.453)
DD M(0) EE = DD M(t′′) EE , (17.454) which directly leads to Eq. (17.448). 34The proof here follows the basic idea of Richard Durrett, Stochastic Calculus: A Practical Introduction (CRC Press, 1996), pp. 137-41.
Chapter 17. Stochastic Processes 17.11.3 Proof: Diffusion with Source
generate is
2\partial 2 x f + g(x, t), (17.455) and we want to show that this PDE is solved by
**
Z t dt′ g[x + W(t′), t −t′] ++ . (17.456) Again, consider the quantity35
Z t dt′ g[x + W(t′), t′′ −t′], (17.457) with f(x, t) satisfying Eq. (17.455). The differential is
2\partial 2 x f[x + W(t), t′′ −t] dt
(17.458) Again, we have that M(t) is a martingale, d DD M(t) EE = 0, (17.459) which upon integration gives DD M(0) EE = DD M(t′′) EE , (17.460) where DD M(0) EE
(17.461) and DD M(t′′) EE = **
Z t′′ dt′ g[x + W(t′), t′′ −t′] ++ , (17.462) thus establishing the desired solution (17.456). 17.11.4 Proof: Diffusion with Damping and Source Now we return to the general case, Eqs. (17.438) and (17.440). In analogy with the simpler cases, we consider the candidate quantity
− Z t dt′ V [x + W(t′), t′′ −t′] + Z t dt′ g[x + W(t′), t′′ −t′] exp
− Z t′ dt′′′ V [x + W(t′′′), t′′ −t′′′] ! . (17.463) 35Richard Durrett, op. cit., pp. 130-6.
17.11 Feynman–Kac Formula The differential is¶
−\partial tf[x + W(t), t′′ −t] dt + \partial xf[x + W(t), t′′ −t] dW + 1 2\partial 2 x f[x + W(t), t′′ −t] dt exp − Z t dt′ V [x + W(t′), t′′ −t′] −f[x + W(t), t′′ −t] V [x + W(t), t′′ −t] exp − Z t dt′ V [x + W(t′), t′′ −t′] dt + g[x + W(t), t′′ −t] exp − Z t dt′′′ V [x + W(t′′′), t′′ −t′′′] dt
(17.464) after using the PDE (17.438) as usual to eliminate terms. Once again, M(t) is a martingale, d DD M(t) EE = 0, (17.465) so that DD M(0) EE = DD M(t′′) EE , (17.466) where DD M(0) EE
(17.467) and DD M(t′′) EE = ** f[x + W(t), 0] exp
− Z t′′ dt′ V [x + W(t′), t′′ −t′] ! + Z t′′ dt′ g[x + W(t′), t′′ −t′] exp
− Z t′ dt′′′ V [x + W(t′′′), t′′ −t′′′]
, (17.468) which establishes the solution (17.440). 17.11.5 Other Forms Other forms of the Feynman–Kac theorem are common, for example, that specify final conditions for the diffusion equation, or that employ more complicated diffusions. As a simple example of the latter, consider
x f, (17.469) (generalized diffusion equation) which is solved by
DD f0[y(t)] EE
(Feynman–Kac path average) (17.470)
(trajectories for Feynman–Kac equation) (17.471) with no explicit time dependence in the SDE. Defining
(17.472)
Chapter 17. Stochastic Processes the differential is
2\partial 2 x f[y(t), t′ −t] (dy)2
- 1
x f[y(t), t′ −t] dt
(17.473)
DD M(0) EE
(17.474) and DD M(t′) EE = DD f[y(t′), 0] EE , (17.475) establishing (17.470). 17.11.5.1 General form for State-Dependent Diffusion For the same state-dependent drift and diffusion represented by the SDE (17.470), but adding decay and source terms to the SDE as in
x f −V (x, t)f + g(x, t), (PDE for generalized Feynman–Kac formula) (17.476) the corresponding Feynman–Kac formula becomes
** f0[y(t)] exp − Z t dt′ V [y(t′), t −t′] + Z t dt′ g[y(t′), t −t′] exp
− Z t′ dt′′ V [y(t′′), t −t′′]
. (generalized Feynman–Kac formula) (17.477) The proof for this goes in the same way by considering the martingale function
− Z t dt′ V [y(t′), t′′ −t′] + Z t dt′ g[y(t′), t′′ −t′] exp
− Z t′ dt′′′ V [y(t′′′), t′′ −t′′′] ! , (17.478)
17.11 Feynman–Kac Formula where the differential is¶
−\partial tf[y(t), t′′ −t] dt + \partial xf[y(t), t′′ −t] dy + 1 2\partial 2 x f[y(t), t′′ −t] dy2 \times exp − Z t dt′ V [y(t′), t′′ −t′] −f[y(t), t′′ −t] V [y(t), t′′ −t] exp − Z t dt′ V [y(t′), t′′ −t′] dt + g[y(t), t′′ −t] exp − Z t dt′′′ V [y(t′′′), t′′ −t′′′] dt = −\partial tf[y(t), t′′ −t] dt + h
i \partial xf[y(t), t′′ −t] + 1
x f[y(t), t′′ −t] dt exp − Z t dt′ V [y(t′), t′′ −t′] −f[y(t), t′′ −t] V [y(t), t′′ −t] exp − Z t dt′ V [y(t′), t′′ −t′] dt + g[y(t), t′′ −t] exp − Z t dt′′′ V [y(t′′′), t′′ −t′′′] dt
(17.479) We have already implemented the diffusion equation (17.476) as usual to show that this is a martingale,
explicitly on time, because then the time dependences do not match in the correct way to permit the diffusion equation to cancel all the deterministic terms of dM(t). Then setting DD M(0) EE = DD M(t′′) EE , (17.480) where DD M(0) EE
(17.481) and DD M(t′′) EE = ** f[y(t′′), 0] exp
− Z t′′ dt′ V [y(t′), t′′ −t′] ! + Z t′′ dt′ g[y(t′), t′′ −t′] exp
− Z t′ dt′′′ V [y(t′′′), t′′ −t′′′]
, (17.482) and dropping the primes on t′′, we arrive at the result (17.478). 17.11.5.2 Time-Dependent Drift and Diffusion: First Form Again, for
(trajectories for Feynman–Kac equation) (17.483) the approach from the previous section does not carry through because the explicit time dependence of \alpha and \beta runs forward in time, but the time dependence of f runs backwards in time. However, both should be the same, as they appear in the PDE. Thus consider the alternate martingale function
− Z t dt′ V [y(t′′ −t′), t′′ −t′] + Z t dt′ g[y(t′′ −t′), t′′ −t′] exp
− Z t′ dt′′′ V [y(t′′ −t′′′), t′′ −t′′′] ! , (17.484)
Chapter 17. Stochastic Processes with differential
−\partial tf[y(t′′ −t), t′′ −t] dt + \partial xf[y(t′′ −t), t′′ −t] dy + 1 2\partial 2 x f[y(t′′ −t), t′′ −t] dy2 \times exp − Z t dt′ V [y(t′′ −t′), t′′ −t′] −f[y(t′′ −t), t′′ −t] V [y(t′′ −t), t′′ −t] exp − Z t dt′ V [y(t′′ −t′), t′′ −t′] dt + g[y(t′′ −t), t′′ −t] exp − Z t dt′′′ V [y(t′′ −t′′′), t′′ −t′′′] dt = −\partial tf[y(t′′ −t), t′′ −t] dt − h \alpha(y, t′′ −t) dt + \beta[y(t), t′′ −t] dW(t′′ −t) i \partial xf[y(t), t′′ −t] + 1
x f[y(t), t′′ −t] dt exp − Z t dt′ V [y(t′′ −t′), t′′ −t′] −f[y(t′′ −t), t′′ −t] V [y(t′′ −t), t′′ −t] exp − Z t dt′ V [y(t′′ −t′), t′′ −t′] dt + g[y(t′′ −t), t′′ −t] exp − Z t dt′′′ V [y(t′′ −t′′′), t′′ −t′′′] dt
(17.485) Note that all the evolution now is explicitly backwards, including the time dependence of the trajectories y(t′′ −t). In the last step, we used the PDE
x f −V (x, t)f + g(x, t), (PDE for generalized Feynman–Kac formula) (17.486) which gives the relevant diffusion equation. As usual, setting DD M(0) EE = DD M(t′′) EE , (17.487) where DD M(0) EE
(17.488) and DD M(t′′) EE = ** f[y(0), 0] exp
− Z t′′ dt′ V [y(t′′ −t′), t′′ −t′] ! + Z t′′ dt′ g[y(t′′ −t′), t′′ −t′] exp
− Z t′ dt′′′ V [y(t′′ −t′′′), t′′ −t′′′]
, (17.489) gives the generalized Feynman–Kac formula
** f0[y(0)] exp − Z t dt′ V [y(t −t′), t −t′] + Z t dt′ g[y(t −t′), t −t′] exp
− Z t′ dt′′ V [y(t −t′′), t −t′′]
. (generalized Feynman–Kac formula) (17.490)
17.11 Feynman–Kac Formula The time dependence can be simplified somewhat by writing
** f0[y(0)] exp − Z t dt′ V [y(t′), t′] + Z t dt′ g[y(t′), t′] exp − Z t t−t′dt′′ V [y(t′′), t′′]
. (generalized Feynman–Kac formula) (17.491) However, the form (17.490) is somewhat better for interpreting the formula, because it explicitly traces each
weights the average in the first term by the initial distribution f0 at the corresponding state y(0). The second term corresponds to the creating of a trajectory by the source g(x, t) at time t′ (summed over all possible
integrated over the temporal extent of each path. The backwards-propagating nature of the solutions is reflected in the above PDE, where the drift coefficient has a minus sign compared to the earlier, ‘‘forward- propagating’’ PDE (17.476). (The sign change does not apply to the diffusion term, since the backward propagation is conditioned on a final condition, not an initial condition.) Note also that while more general than Eq. (17.477) in accounting for explicit time dependence in the SDE coefficients, Eq. (17.477) is much more convenient, from the point of view of simulation—the trajectories here have to be propagated backwards from a particular final point (in which case the SDE acts as an anticipating SDE), and weighted according to the corresponding zero-time point. 17.11.5.3 Time-Dependent Drift and Diffusion: Second Form An alternative, ‘‘forward-time’’ version of the generalized Feynman–Kac formula arises by considering the forward-propagating martingale function
− Z t dt′ V [y(t′), t′] + Z t dt′ g[y(t′), t′] exp
− Z t′ dt′′ V [y(t′′), t′′] ! . (17.492) The differential is
\partial tf[y(t), t] dt + \partial xf[y(t), t] dy + 1 2\partial 2 x f[y(t), t] dy2 exp − Z t dt′ V [y(t′), t′] −f[y(t), t] V [y(t), t] exp − Z t dt′ V [y(t′), t′] dt + g[y(t), t] exp − Z t dt′′ V [y(t′′), t′′] dt = \partial tf[y(t), t] dt + h
i \partial xf[y(t), t] + 1
x f[y(t), t] dt exp − Z t dt′ V [y(t′), t′] −f[y(t), t] V [y(t), t] exp − Z t dt′ V [y(t′), t′] dt + g[y(t), t] exp − Z t dt′′ V [y(t′′), t′′] dt
(17.493) where in the last step we used the PDE
x f −V (x, t)f + g(x, t). (PDE for generalized Feynman–Kac formula) (17.494)
Chapter 17. Stochastic Processes Then as we have done so many times, we set DD M(0) EE = DD M(t) EE , (17.495) where DD M(0) EE
(17.496) and DD M(t) EE = ** f[y(t), t] exp − Z t dt′ V [y(t′), t′] + Z t dt′ g[y(t′), t′] exp
− Z t′ dt′′ V [y(t′′), t′′]
, (17.497) leading to the alternate generalized formula
** f[y(t), t] exp − Z t dt′ V [y(t′), t′] + Z t dt′ g[y(t′), t′] exp
− Z t′ dt′′ V [y(t′′), t′′]
. (generalized Feynman–Kac formula) (17.498) Like Eq. (17.491), this form is in a sense less convenient for simulation than Eq. (17.477), because it obtains the initial distribution from the final distribution. The reverse-time nature of this formula is reflected in the negative time derivative in the PDE (17.494). Note also that the diffusion equation (17.494) has the form of the Kolmogorov backward equation (17.144). Thus, the derivation of the above is effectively an alternate derivation of the Kolmogorov backward equation, for the evolution of P(x, t|x0, t0) with t0, given a final distribution at t. This is implicit in the form of the solution (17.498), which gives the initial distribution at t0, given the final distribution at t, such that the time derivative in the PDE (17.494) can be regarded here as a derivative with respect to t0, and the distribution f(x, t) should be regarded as equivalent to P(x, t|x0, t0). 17.12 Sojourn Times 17.12.1 Wiener Process As an example application of the Feynman–Kac formula, we wish to consider the sojourn time of a Wiener path W(t) past a boundary at d, which is defined as the functional
Z t dt′ \Theta[W(t′) −d]
(17.499) (sojourn time) or the total time that W(t) spends across the boundary. In the plot below, this counts the portion of time that the path is highlighted in green.
17.12 Sojourn Times t Wo(t) d We will only consider the case d \ge 0, since we can get the d < 0 probabilities as complements to the probabilities we calculate.
(17.500) Then the Laplace transform of fTs(x) is Z t
exp {−sTs[W(t); d]}
=
exp −s Z t dt′ \Theta[W(t′) −d] . (17.501) Our goal will be to compute the expectation value on the right via the Feynman–Kac formula, and then to obtain fTs(x) by inverting the Laplace transform.36 Consider now the driven diffusion equation
2\partial 2 x f −V (x)f −\lambdaf + g(x), (17.502) where V (x) is the occupation function we seek here (i.e., the step function)—we will take it to be the constant s past the barrier, and 0 before it:
(17.503) The Feynman–Kac formula (17.440) gives the solution of this equation as
** f0[x + W(t)] exp −\lambdat − Z t dt′ V [x + W(t′)] + Z t dt′ g[x + W(t′)] exp
−\lambdat′ − Z t′ dt′′ V [x + W(t′′)]
. (17.504)
f0(x), so we will take this initial condition to be zero. Thus, we have the steady state
dt exp −\lambdat − Z t dt′ V [x + W(t′)]
= Z \infty dt e−\lambdat ** exp −s Z t dt′ \Theta[x + W(t′) −d]
, (17.505) 36here we are adapting the method of Gerard Hooghiemstra, ‘‘On explicit occupation time distributions for Brownian pro- cesses,’’ Statistics & Probability Letters 56, 405 (2002) (doi: 10.1016/S0167-7152(02)00037-8).
Chapter 17. Stochastic Processes
This contains our desired expectation value from
expectation value. This is then the solution of the steady-state version of Eq. (17.502):
2\partial 2 x f(x) −s\Theta(x −d)f(x) + 1. (17.506) For x > d, the ODE is
(17.507)
that for x > d,
p
(17.508) or
Ae− p
(x > d) Be \sqrt
\lambda (x < d) (17.509) for some undetermined constants A and B, where we have chosen the bounded solution on either side of d.
Ae− p
\sqrt
\lambda, (17.510)
− p
p
\sqrt 2\lambdaBe \sqrt 2\lambda d. (17.511) The solution of these two equations fixes the coefficients as A = se p
p
, B = − se− \sqrt 2\lambda d
p
. (17.512)
Z \infty dt e−\lambdat ** exp −s Z t dt′ \Theta[W(t′) −d]
\lambda = 1 \lambda − se− \sqrt 2\lambda d
p
. (17.513) Then using Eq. (17.501) on the left-hand side, Z \infty dt e−\lambdatDD exp {−sTs[W(t)]} EE = 1 \lambda − se− \sqrt 2\lambda d
p
. (17.514) Now we can use the formula for the inverse Laplace transform [Problem 15.4]
2\pii
0+−i\infty ds es\tau L y, (17.515) where L y is the Laplace transform of y(t):
Z \infty dt e−st y(t). (17.516) Then we have DD exp {−sTs[W(t)]} EE = 2\pii
0+−i\infty
" \lambda − se− \sqrt 2\lambda d
p
. (17.517)¶
17.12 Sojourn Times Evaluating the first term is simple by rotating the integration direction and then completing the contour integral around the great half-plane: 2\pii
0+−i\infty
2\pii Z \infty−i0+
\lambda = 1, (t > 0). (17.518) The second term is more involved, and to evaluate it we will stick to ‘‘recipes’’ for known Laplace transforms. We will need37 L e−st/2
\sqrt \lambda
p
, (17.519) where I\nu(x) is a modified Bessel function, and38 L 1 \sqrt
\sqrt \lambda e−k \sqrt \lambda, (17.520) along with the convolution theorem for Laplace transforms,39 L Z t
= L [f] L [g] , (17.521) which combine to give the Laplace transform appropriate to the second term in Eq. (17.514), L "Z t
p
\sqrts
¶
e−k \sqrt \lambda \lambda
p
, (17.522) provided we take k = − \sqrt 2 d and insert an overall factor of −s. Thus, Eq. (17.514) becomes Z t
r s \pi Z t
(17.523) after using Eq. (17.501) to replace the ensemble average on the left-hand side. Now we will use an integral representation of the Bessel function,40
r x 2\pi Z 1 −1 du e−ux = r 2x \pi ex Z 1 dv e−2vx = r 2x \pi ex t Z t dv e−2vx/t, (17.524)
becomes Z t
\pit Z t d\tau Z t
r \tau t −\tau
\pi Z t d\tau e−s\tau −1
p
\pi Z t
p
= Z t
\pi Z t dx e−sx e−d2/2(t−x) p x(t −x) , (17.525) 37Milton Abramowitz and Irene A. Stegun, Handbook of Mathematical Functions (Dover, 1965), p. 1024, Eq. (29.3.54). 38Milton Abramowitz and Irene A. Stegun, op. cit., p. 1026, Eq. (29.3.84). 39Milton Abramowitz and Irene A. Stegun, op. cit., p. 1020, Eq. (29.2.8). 40Milton Abramowitz and Irene A. Stegun, op. cit., p. 376, Eq. (9.6.18).
Chapter 17. Stochastic Processes where in the last step we set \tau = x, and we have also named the integral
\pi Z t dx e−d2/2(t−x) p x(t −x) , (17.526) which we will evaluate shortly. Eq. (17.525) now equates (finite-time) Laplace transforms. Differentiating with respect to time undoes them, which then yields
\pi e−d2/2(t−x) p x(t −x) . (17.527) Clearly, this result is normalized, since integration over x produces canceling terms of I(d, \tau) from the \delta function and the integral. The \delta function indicates that there is a (possibly) finite probability for having zero sojourn time, whereas positive sojourn times have infinitesimal probability, as we expect for a distribution. In particular, the probability for having zero sojourn time is just
ϵ\rightarrow 0 Z ϵ dx fTs(x)
(17.528) However, we have already calculated this: the probability to not sojourn across the boundary at d is equivalent to the non-crossing probability for a boundary at d. This is the complement of the crossing probability from Eq. (17.372), and thus
d \sqrt 2t = erf d \sqrt 2t . (17.529) Thus, in computing the boundary-crossing probability before, we have essentially used a path-integral method to evaluate the integral (17.526), with the result
d \sqrt 2t . (17.530) Putting this result into Eq. (17.525), we finally have the probability density for the sojourn time of W(t) beyond d:41
d \sqrt 2t
\pi p x(t −x)
(probability density for sojourn time of W(t) past d) (17.531) The probability density for the time the particle stays under the barrier at d is given by the replacement x −\rightarrow t −x. The cumulative probability is given by simply integrating this from 0 to x:
d \sqrt 2t + 1 \pi Z x dx′ e−d2/2(t−x′) p x′(t −x′)
(cumulative probability for sojourn time of W(t) past d) (17.532)
nonintersecting regions in which to sojourn. 41cf. Andrei N. Borodin and Paavo Salminen, op. cit., p. 156, formula 1.4.4, though note the absence of the first (\delta-function) term; Paul Lévy, ‘‘Sur certains processus stochastiques homogènes,’’ Compositio Mathematica 7, 283 (1940), p. 327, Eq. (58)
Applied Probability 32, 405 (1995), Eq. (3.a) (doi: 10.2307/3215296); Lajos Takács, ‘‘Sojourn Times,’’ Journal of Applied Mathematics and Stochastic Analysis 9, 415 (1996), Eq. (3) (doi: 10.1155/S1048953396000366).
17.12 Sojourn Times The moments are given simply by integrating the density (17.531) as DD T n s EE = Z t dx xnfTs(x) = Z t dx xn e−d2/2(t−x) \pi p x(t −x)
(17.533)
becomes DD T n s EE = 2tn \pi Z 1 d\sigma
(17.534) Thus, for example, the mean sojourn time is DD Ts EE = 1 2(d2 + t) erfc d \sqrt 2t − r td2
(17.535) Again, for d < 0, this expression for d \ge 0 can be adapted by subtracting it from t. 17.12.1.1 Arcsine Laws
\pi Z x dx′ p x′(t −x′) . (17.536) Evaluating this integral gives
\pi sin−1 rx t
(Levy’s arcsine law) (17.537) which is known as Lévy’s arcsine law.42 This gives the probability distribution of the time a Wiener path spends on a particular side of its starting point. The density following from Eq. (17.536) is
\pi p x(t −x) . (17.538) (density for Levy arcsine law) Since this density diverges at x = 0 and x = t, it says that the path is most likely to spend essentially all its time on one side or the other of the origin. 17.12.1.2 Alternate Arcsine Law: Last Crossing time The same arcsine law also applies to other statistics related to Wiener processes. To explore the first statistic, suppose we fix a time t, and ask what was the most recent time before t when W(t) crossed the origin?43 That is, consider the statistic
(17.539)
possible values x of W(\tau), using the Gaussian probability density for W(\tau): P
= Z \infty −\infty
(17.540) 42Paul Lévy, ‘‘Sur un problème de M. Marcinkiewicz,’’ Comptes rendus hebdomadaires des séances de l’Académie des sci- ences 208, 318 (1939) (http://gallica.bnf.fr/ark:/12148/bpt6k3160g/f300.image); Paul Lévy, ‘‘Sur certains processus stochastiques homogènes,’’ op. cit. 43The proof here is from Jim Pitman’s course notes, http://www.stat.berkeley.edu/~pitman/s205s03/lecture18.pdf, Theorem 18.4.
Chapter 17. Stochastic Processes
Since the Gaussian probability density is even, we can change to evaluating only half of the integral x \ge 0: P
= 2 Z \infty dx fW (\tau)(x) P W(t′) < x ∀t′\in [0,t−\tau) . (17.541) Note we have also changed the second probability statement: it was that a Wiener process, starting at x, should not touch the origin for a time interval of length t −\tau. This is the same as the probability that a Wiener process, starting at 0, should not touch x over an interval of the same duration. The latter probability is just the complement of a crossing probability, which we can obtain from the complement of Eq. (17.372). Putting in this probability as an error function, and also putting in the Gaussian density, we have P
= 2 Z \infty dx \sqrt
x p 2(t −\tau) ! . (17.542) Evaluating this integral, we have P
= 2 \pi sin−1 r\tau t . (second arcsine law for last crossing time) (17.543) Note that this implies the density (17.538) for \tauL0(t), which is someone odd in that it is a symmetric function over [0, t]. This means that the most recent crossing most likely occurred close to t′ = t or close to t′ = 0, and the symmetry is peculiar, given that W(t′) itself does not share the symmetry. Note that while the last crossing is well-defined, asking something like when was the next-to-last crossing occur is much more complicated: the crossings of W(t′) form a set of measure zero, but are nevertheless uncountable, much like the Cantor set. Thus, the last crossing is not isolated, but rather has other, arbitrarily close crossings. 17.12.1.3 Alternate Arcsine Law: Maximum-Value Time The remaining statistic governed by the arcsine law is the time \taumax at which a Wiener path attains its maximum value over the time interval [0, t]. To see this, it is first necessary to make some comments about the maxima of a Wiener process.44 First, note that over any time interval [a, b], the maximum of W(t) over the interval does not occur at either endpoint (with unit probability). This follows at the lower endpoint a
and the crossing probability (17.372): over an interval of any finite duration, the probability of crossing a barrier just above the starting point is arbitrarily close to 1. The maximum does not occur at the upper endpoint b for the same reason (i.e., the same argument under a time-reversal). By subdividing the [a, b] into smaller and smaller intervals, each of which has its own local maximum, we can see that the local maxima are dense (but countable). Also by subdividing [a, b] into two intervals (which do not overlap except at the boundary point), each subinterval will have a local maximum. But because the Wiener process progresses independently in each interval, the local maximum can take on any value, and in particular the probability for attaining the same local maximum in each interval is zero. The point is that the maximum over [a, b] is unique, so that \taumax is well-defined. Now consider a Wiener process W(t), with the maximum-value process
(17.544) representing the running maximum value attained by W(t) up to time t. Defining the difference process
(17.545) we will aim to show that R(t) is statistically equivalent45 to the absolute-value process or reflection Brownian motion |W(t)| (about which we will have more to explore in Section 17.13.4.1). To do this, 44Peter Mörters and Yuval Peres, Brownian Motion (Cambridge, 2010) (ISBN: 0521760186). 45This argument follows that of Peter Mörters and Yuval Peres, op. cit., at least until near the end.
17.12 Sojourn Times consider a fixed time t0 > 0, and an evolution over a time interval ∆t. The Wiener process changes in this interval by
(17.546) and the maximum changes by
max
(17.547) The idea will be to show that R(t0 + ∆t) is statistically equivalent to |R(t0) + ∆W0(∆t)|. First, note that the maximum at t0 + ∆t is
n M(t0), W(t0) + ∆M0(∆t) o . (17.548) That is, the final maximum is either the initial maximum, or the initial value of W(t) plus the maximum added over the interval ∆t, whichever is larger. Now writing
= max n M(t0), W(t0) + ∆M0(∆t) o − h W(t0) + ∆W0(∆t) i = max n R(t0), ∆M0(∆t) o −∆W0(∆t). (17.549) Now to analyze the two cases here. If R(t0) is the larger value in the maximum, then this reduces to
maximum value that we are subtracting from R(t0) is less than R(t0), by the definition of the maximum and ∆M0(∆t). In the other case, ∆M0(∆t) > R(t0), which means that R(t0) −∆W0(∆t) crosses through zero, and so this case requires more careful treatment. In fact it is useful to identify the time \deltatmax at which the
a local maximum, according to its definition. Then continuing the evolution, R(t0 + ∆t) is the magnitude of the difference between W(t0 + \deltatmax) and W(t0 + ∆t). So R(t0 + ∆t) is still positive as required, and if we apply the reflection principle to ∆W0(∆t) at this time t0 + \deltatmax, the reflected version of R(t0 + ∆t) is exactly equivalent to a Wiener process that happened to cross through zero at time t0 + \deltatmax. Returning to the main point, since M(t) −W(t) behaves like |W(t)|, the last zero of which correspond to the time \taumax when W(t) achieves its maximum value. Since |W(t)| has the same set of zeros as W(t), and the last zero has the arcsine distribution, the time \taumax of the maximum likewise has the arcsine distribution (17.543). Somewhat counterintuitively, then, the Wiener process tends to achieve its maximum value near the beginning or near the end of the time interval in question. 17.12.2 Standard Brownian Bridge Now we will repeat the above analysis to compute the probability distributions for the sojourn time of a Brownian bridge above a boundary at d.46 As before, in the plot here, this counts the portion of time that
46this is the first calculation in Gerard Hooghiemstra, op. cit.
Chapter 17. Stochastic Processes t Bo(t) d
functional
Z 1 dt \Theta[B(t) −d], (sojourn time for standard Brownian bridge) (17.550) or the total time that B(t) spends across the boundary. The calculation here is similar to the calculation for the Wiener path, and we will refer to that calculation frequently to save effort. The main difference is that
with respect to k to introduce a \delta function that will pin the endpoint of the Wiener path to make a bridge. This will lead to somewhat more complicated algebra, but nothing conceptually different. The procedure is the same as for the Wiener path up to Eq. (17.505), which we will leave as
Z \infty dt e−\lambdat ** g[x + W(t)] exp −s Z t dt′ \Theta[x + W(t′) −d]
. (17.551)
2\partial 2 x f(x) −s\Theta(x −d)f(x) + eikx, (17.552) which for x > d is
(17.553)
2k2eikx
2k2eikx
2k2eikx
(17.554) so that for x > d,
p
(17.555)
17.12 Sojourn Times or¶
Ae− p
2eikx
(x > d) Be \sqrt
2eikx
(x < d), (17.556) where A and B must be determined by requiring continuity of f and f ′ at x = d. Solving the two resulting equations gives A = \sqrt 2 s \sqrt 2\lambda −ik e p
\sqrt
\sqrt
, B = − \sqrt 2 s p
e− \sqrt
\sqrt
\sqrt
. (17.557)
Z \infty dt e−\lambdat ** eikW (t) exp −s Z t dt′ \Theta[W(t′) −d]
(17.558) Integrating with respect to k and then dividing through by 2\pi, Z \infty dt e−\lambdat ** \delta[W(t)] exp −s Z t dt′ \Theta[W(t′) −d]
= 1 2\pi Z \infty −\infty dk B +
. (17.559)
bridges). We will treat this carefully to obtain the correct normalization for the ensemble average. In the following, we will use t′ as a dummy variable and t as the temporal endpoint, so that 0 \le t′ \le t. We have an ensemble average of a functional F[W(t′)] (note that \delta[W(t)] only operates on the final value W(t), and so is not itself a functional),
\delta[W(t)] F[W(t′)]
= X W (t′) P[W(t′)] \delta[W(t)] F[W(t′)] = X
dW F[W(t′)], (17.560) where we have expanded into a sum over all possible paths, weighted by the probability P[W(t′)] of W(t′) to occur, and then implemented the \delta function in selecting the probabilities for Wiener paths satisfying
The factor of 1/dW comes from thinking of
\delta[W(t)] F[W(t′)]
= X {W (t′)}
dW F[W(t′)] = X
\sqrt 2\pit dW F[W(t′)] = X
\sqrt 2\pit F[W(t′)] = X {Bt(t′)} P[Bt(t′)] \sqrt 2\pit F[Bt(t′)] = \sqrt 2\pit
F[Bt(t′)]
, (17.561)
Chapter 17. Stochastic Processes where we have switched to Brownian bridges, which are equivalent to pinned Wiener paths, and we are using
\sqrt
Changing the ensemble average in Eq. (17.559) to Brownian bridges, we have Z \infty dt \sqrt t e−\lambdat ** exp −s Z t dt′ \Theta[Bt(t′) −d]
= \sqrt 2\pi Z \infty −\infty dk B +
, (17.562)
\sqrt T B(t/T) for a bridge that closes at time T. Then using Eq. (17.501) on the left-hand side, adapted for the bridge, and (17.557) on the right-hand side, Z \infty dt \sqrt t e−\lambdat
exp {−sTs[Bt]}
= \sqrt 2\pi Z \infty −\infty dk
\sqrt 2 s p
e− \sqrt
\sqrt
\sqrt
= r\pi
\sqrt 2\lambda d \sqrt \lambda \sqrt
\sqrt \lambda \sqrt
\sqrt \lambda ! . (17.563) To invert the Laplace transform here, we will use the formula47 L 1 t e−st/2I1(st/2)
\sqrt
\sqrt \lambda \sqrt
\sqrt \lambda (17.564)
\sqrt 2 d) and the Laplace convolution rule (17.521) to write L "Z t d\tau p
= e−2 \sqrt 2\lambda d \sqrt \lambda \sqrt¶
\sqrt \lambda \sqrt
\sqrt \lambda ! , (17.565) which takes care of the second term. The first term on the right-hand side of (17.563) is also covered by the
exp {−sTs[Bt]}
= 1 − Z t d\tau 1 \tau r t
(17.566)
and using Eq. (17.501) on the left-hand side, Z 1
Z 1 d\tau
(17.567) We will proceed as before by using the integral representation of the Bessel function,48
\pi Z 1 −1 du p 1 −u2 e−ux = 4x \pi ex Z 1 dv p v(1 −v) e−2vx, (17.568)
Z 1
\pi Z 1 d\tau Z 1 dv r v(1 −v) 1 −\tau
(17.569) and letting v = x/\tau, Z 1
\pi Z 1 d\tau Z \tau dx 1 \tau 2 r x(\tau −x) 1 −\tau
= Z 1 dx e−sx \delta(x −0+) −2s \pi Z 1 dx e−sx Z 1 x d\tau 1 \tau 2 r x(\tau −x) 1 −\tau
(17.570) 47Milton Abramowitz and Irene A. Stegun, op. cit., p. 1024, Eq. (29.3.52). 48Milton Abramowitz and Irene A. Stegun, op. cit., p. 376, Eq. (9.6.18).
17.12 Sojourn Times Note how we changed the integration limits when interchanging the order of integration in the last step, in order to integrate over the same triangular area in the (x, \tau)-plane. We can regard fTs(x) here as vanishing for x > 1, since it is nonsensical to consider the standard bridge past t = 1. If we regard the last term in the same way, we can extend the upper integration limits to \inftyand write out the Laplace transforms as
−sL " \pi Z 1 x d\tau 1 \tau 2 r x(\tau −x) 1 −\tau
. (17.571) Recall [Eq. (5.144)] that the Laplace transform of a derivative satisfies¶
(17.572) so that
−L " \partial x 2\sqrtx \pi Z 1 x d\tau r \tau −x 1 −\tau
\tau 2 # −I(d) = L
−L " \partial x 2\sqrtx \pi Z 1 x d\tau r \tau −x 1 −\tau
\tau 2 # , (17.573) where
x\rightarrow 0 2\sqrtx \pi Z 1 x d\tau r \tau −x 1 −\tau
\tau 2 . (17.574) Thus, we can write the sojourn-time density as
" 2\sqrtx \pi Z 1 x d\tau r \tau −x 1 −\tau
\tau 2 # . (17.575) Again, the [1 −I(d)] coefficient of the \delta function is associated with the probability of having a zero sojourn time. This is the probability of a bridge to not touch the boundary. We have already calculated the touching probability in Eq. (17.380), so
(17.576) Thus, our first complete expression for the sojourn-time density is
h 1 −e−2d2i
" 2\sqrtx \pi Z 1 x d\tau r \tau −x 1 −\tau
\tau 2 #
(probability density for bridge sojourn time) (17.577) with cumulative probability49
\pi Z 1 x d\tau r \tau −x 1 −\tau
\tau 2
(cumulative probability distribution for bridge sojourn time) (17.578) though we will continue a bit more in massaging the expression here into nicer forms. Next, we change variables via
(17.579) with the result
h 1 −e−2d2i
4(1 −x) \pi e−2d2/(1−x) Z \infty d\sigma \sigma2
. (probability density for bridge sojourn time) (17.580) 49Gerard Hooghiemstra, op. cit., Eq. (6).
Chapter 17. Stochastic Processes Note that here we can identify
x\rightarrow 0 4(1 −x) \pi e−2d2/(1−x) Z \infty d\sigma \sigma2
= 4 \pi e−2d2 Z \infty d\sigma \sigma2
= 4 \pi e−2d2 Z \infty d\sigma \sigma2
= e−2d2, (17.581) which serves as a separate verification of the boundary-crossing result (17.380). Additionally, we can now integrate Eq. (17.580) from 0 to x to obtain the cumulative distribution for the sojourn time. In doing so, the integral of the last term, evaluated at the lower integration limit, will cancel the e−2d2 from the first term, and the result is50
\pi e−2d2/(1−x) Z \infty d\sigma \sigma2
(cumulative probability distribution for bridge sojourn time) (17.582) Finally, we can evaluate the \sigma integral, with the result51
r x(1 −x) 2\pi e−2d2/(1−x) −(1 −x + 4d2x) e−2d2 erfc r 2d2x 1 −x !
(cumulative probability distribution for bridge sojourn time) (17.583) The same integration applied to (17.580) gives the explicit probability density52
h 1 −e−2d2i
r 8d2(1 −x) \pix e−2d2/(1−x) + (1 −4d2) e−2d2 erfc r 2d2x 1 −x ! . (probability density for bridge sojourn time) (17.584)
\sqrt
replacements d −\rightarrow d/ \sqrt T and x −\rightarrow x/T. Also, recall that we explicitly assumed d \ge 0; for d < 0, this formula for fTs(x) can be regarded as the density for 1 −Ts. We can then compute the moments for the sojourn time as DD (Ts)nEE = Z 1 dx xnfTs(x) = 1 −n Z 1
= n Z 1
= 2n \pi Z 1 dx xn−1/2 Z 1 x d\tau r \tau −x 1 −\tau
\tau 2 = 2n \pi Z 1
Z \tau
\tau −x, (17.585) where we integrated by parts to integrate with respect to the cumulative distribution, for which we used the 50Gerard Hooghiemstra, op. cit., before Eq. (10). 51Lajos Takács, ‘‘The Distribution of the Sojourn Time for the Brownian Excursion,’’ Methodology and Computing in Applied Probability 1, 7 (1999), Eq. (6) (doi: 10.1023/A:1010060107265). 52cf. Andrei N. Borodin and Paavo Salminen, op. cit., p. 158, formula 1.4.8.
17.12 Sojourn Times form (17.578). After performing the x integration, the result is53 DD (Ts)nEE
Z 1
. (17.586) Notice that this moment formula also arises from the moment-generating function in the form (17.566), where the power-series expansion of e−s\tau/2I1(s\tau/2) in s yields the individual moments here. A somewhat easier form to handle arises by changing variables via \tau = 1 −\sigma2, with the result DD (Ts)nEE
Z 1
(17.587) (sojourn-time moments) The integral here does not have a simple general form, but it can be readily evaluated for particular n. For example, we have DD Ts EE = e−2d2 − r\pi 2 d erfc h\sqrt 2 d i DD (Ts)2EE
− \sqrt 2\pid 1 + 4d2 erfc h\sqrt 2 d i DD (Ts)4EE
− \sqrt 2\pid
erfc h\sqrt 2 d i (17.588) for the first, second, and fourth moments of the sojourn time, remembering that d \ge 0 in these expressions. 17.12.3 Brownian Bridge
so that in this section we will consider the sojourn-time functional
Z 1 dt \Theta[B0\rightarrow c(t) −d]. (sojourn time for Brownian bridge) (17.589)
is then essentially the same, with the result with exp(ikx) −\rightarrow exp[ik(x −c)] in (17.556), and exp(ikd) −\rightarrow
Z \infty dt e−\lambdat ** eik[W (t)−c] exp −s Z t dt′ \Theta[W(t′) −d]
(17.590) in place of (17.590), where we should keep in mind that the form of B is modified here. Then integrating with respect to k and dividing through by 2\pi gives Z \infty dt e−\lambdat ** \delta[W(t) −c] exp −s Z t dt′ \Theta[W(t′) −d]
= 1 2\pi Z \infty −\infty dk B + 2e−ikc
. (17.591) in place of Eq. (17.559). Thus we have introduced the correct modification to force the \delta function to pin the
53Gerard Hooghiemstra, op. cit., Eq. (10).
Chapter 17. Stochastic Processes Next, we should generalize the result (17.561) to the case of pinning W(t) to c:
\delta[W(t) −c] F[W(t′)]
= X W (t′) P[W(t′)] \delta[W(t) −c] F[W(t′)] = X
dW F[W(t′)] = X
dW F[W(t′)] = X
\sqrt 2\pit dW F[W(t′)] = X
\sqrt 2\pit F[W(t′)] = X {Bt(t′)} P[Bt(0\rightarrow c)(t′)] e−c2/2t \sqrt 2\pit F[Bt(0\rightarrow c)(t′)]
\sqrt 2\pit
F[Bt(0\rightarrow c)(t′)]
. (17.592) The change here is straightforward, and involves evaluating the Gaussian probability density for W(t) at c instead of 0. Then changing the average in Eq. (17.591) to encompass the appropriate Brownian bridges, we have Z \infty dt e−c2/2t \sqrt t e−\lambdat ** exp −s Z t dt′ \Theta[Bt(0\rightarrow c)(t′) −d]
= \sqrt 2\pi Z \infty −\infty dk B + 2e−ikc
, (17.593) in place of Eq. (17.562). Carrying out the following integration gives Z \infty dt e−c2/2t \sqrt t e−\lambdat
exp −sTs[Bt(0\rightarrow c)]
= \sqrt 2\pi Z \infty −\infty dk 2e−ikc
\sqrt 2 s p
e− \sqrt
\sqrt
\sqrt
= "r\pi \lambda e− \sqrt
\sqrt 2\lambda (2d−c) \sqrt \lambda \sqrt
\sqrt \lambda \sqrt
\sqrt \lambda !#
(17.594)
Z \infty dt e−c2/2t \sqrt t e−\lambdat
exp −sTs[Bt(0\rightarrow c)]
= \sqrt 4\pi e− \sqrt 2\lambda d− p
\sqrt
\sqrt \lambda
(17.595) The first expression (17.594) is the same result as in Eq. (17.563), except for a factor exp(−c2/2t) on the left-hand side, the factor of exp(− \sqrt 2\lambda c) in the first term on the right-hand side, and the replacement 2d −\rightarrow 2d −c in the second term on the right-hand side. To invert the Laplace transform, the second term on the right-hand side inverts in the same way with the replacement 2d −\rightarrow 2d−c, while Eq. (17.520) applies to the first term with k = \sqrt 2 |c|, such that Eq. (17.566) becomes DD exp −sTs[Bt(0\rightarrow c)] EE
Z t d\tau 1 \tau r t
(17.596)
17.12 Sojourn Times or that is, the last term on the right-hand side is multiplied by exp(c2/2t), and subject to the replacement
rest of the treatment in the previous section. In particular, Eq. (17.574) becomes
x\rightarrow 0 2\sqrtx \pi ec2/2 Z 1 x d\tau r \tau −x 1 −\tau
\tau 2
(17.597) which is the correct boundary-crossing probability (17.387) for the same bridge pinned to c. Then we can adapt the probability-density expressions with these modifications, with the result54
h 1 −e−2d(d−c)i
" 2\sqrtx \pi ec2/2 Z 1 x d\tau r \tau −x 1 −\tau
\tau 2 # = h 1 −e−2d(d−c)i
−\partial x 4(1 −x) \pi
Z \infty d\sigma \sigma2
= h 1 −e−2d(d−c)i
r 2(1 −x) \pix
- 1 −(2d −c)2 e−2d(d−c) erfc s (2d −c)2x 2(1 −x) !
(probability density for bridge sojourn time) (17.598) Similarly, the cumulative-probability expressions become
\pi ec2/2 Z 1 x d\tau r \tau −x 1 −\tau
\tau 2
\pi
Z \infty d\sigma \sigma2
r 2x(1 −x) \pi
− 1 −x + (2d −c)2x e−2d(d−c) erfc s (2d −c)2x 2(1 −x) ! , (cumulative probability distribution for bridge sojourn time) (17.599) and the moment formula (17.587) becomes DD (Ts)nEE
Z 1
(sojourn-time moments) (17.600) under the same replacements. For example, the explicit expression for the mean is DD Ts EE
− r\pi 8 (2d −c) ec2/2 erfc 2d −c \sqrt
(17.601) which reduces to the standard-bridge mean in Eqs. (17.588). 54cf. Andrei N. Borodin and Paavo Salminen, op. cit., p. 158, formula 1.4.8.
Chapter 17. Stochastic Processes
L e−st/2
1 −e−st 2s \sqrt \pit3
\sqrt
\sqrt \lambda , (17.602) where I\nu(x) is a modified Bessel function, along with56 L k \sqrt
\sqrt \lambda, (17.603) these transforms can be combined via the convolution formula (17.521) to give L "Z t d\tau k 4\pis p
¶
e−k \sqrt \lambda \sqrt
\sqrt \lambda . (17.604) Then using the shifted-argument version of Eq. (17.603), L k′ \sqrt \pit3 e−st e−k′2/4t
(17.605) we can again employ convolution to combine this with Eq. (17.604) to give L "Z t d\sigma Z t−\sigma d\tau kk′
p
e−s\sigma #
\sqrt
\sqrt
\sqrt \lambda . (17.606) Thus, setting k = \sqrt 2 d and k′ = \sqrt 2 (c −d), we may invert the Laplace transform in Eq. (17.595) to give
exp −sTs[Bt(0\rightarrow c)]
\sqrt t 2\pis Z t d\sigma Z t−\sigma
p
1 −e−s\tau e−s\sigma
(17.607) We can simplify this to t = 1 for a Brownian bridge that runs over a unit time interval:
exp −sTs[B(0\rightarrow c)]
2\pis Z 1 d\sigma Z 1−\sigma
p
1 −e−s\tau e−s\sigma
(17.608) Now using
s =
\sigma dx e−sx (17.609) to remove the difference in the final factor, while introducing a new integral, we can shift the order of integration according to Z 1 d\sigma Z 1−\sigma d\tau
\sigma dx = Z 1 d\sigma Z 1 \sigma dx Z 1−\sigma x−\sigma
Z 1 dx Z x d\sigma Z 1−\sigma x−\sigma d\tau, (17.610) and shifting \tau −\rightarrow \tau −\sigma, Eq. (17.608) becomes
exp −sTs[B(0\rightarrow c)]
2\pi Z 1 dx e−sx Z x d\sigma Z 1 x
p
(17.611) 55Milton Abramowitz and Irene A. Stegun, op. cit., p. 1024, Eq. (29.3.53). 56Milton Abramowitz and Irene A. Stegun, op. cit., p. 1026, Eq. (29.3.82).
17.12 Sojourn Times Now the Laplace transform in s may be inverted to give the probability density:
2\pi Z x d\sigma Z 1 x
p
(17.612) Changing variables via \tau −\rightarrow 1 −\tau gives
2\pi Z x d\sigma Z 1−x d\tau
p
(17.613)
d(c −d) 2\pi p x(1 −x) Z \infty du Z \infty dv uv e−d2v/2(1−x)−(c−d)2u/2x [uv −(1 −x)u −xv]3/2
(17.614) from where it is difficult to proceed with the integration. There is a second approach to inverting the Laplace transforms here that will lead to expression with only a single integral, in both the density and the Laplace transform of the density.57 As an alternate form of the right-hand side in Eq. (17.606), we can consider e−k \sqrt
\sqrt
\sqrt \lambda = 1 s \sqrt
\sqrt \lambda e−k \sqrt
= 1 s e−k \sqrt
−1
\lambda e−k \sqrt \lambda . (17.615) Using again the Laplace transform (17.603) and its shifted version (17.605), along with58 L k2 −2t \sqrt
\sqrt \lambda e−k \sqrt \lambda, (17.616) and the shifted version L k′2 −2t \sqrt \pit5 e−st e−k′2/4t
\sqrt
(17.617) the convolution theorem allows us to combine these transforms to form the right-hand side of Eq. (17.615), with the result L "Z t
8\pis p
¶
\sqrt
\sqrt
\sqrt \lambda . (17.618) Then again setting k = \sqrt 2 d and k′ = \sqrt 2 (c −d), we may invert the Laplace transform in Eq. (17.595) to give
exp −sTs[Bt(0\rightarrow c)]
\sqrt t \sqrt 2\pis Z t
p
(17.619) or for t = 1,
exp −sTs[Bt(0\rightarrow c)]
\sqrt 2\pis Z 1
p
(17.620) 57Vadim Linetsky, ‘‘Step Options,’’ Mathematical Finance 9, 55 (2001) (doi: 10.1111/1467-9965.00063). See in particular Eq. (C.9) and the discussion just before Eq. (A.9).
Chapter 17. Stochastic Processes Then writing the s-dependence as e−s\tau s = Z \infty \tau dx e−sx = Z \infty dx e−sx − Z \tau dx e−sx, (17.621) and then massaging the two resulting integrals (with the rest of the integrand suppressed for brevity) via Z 1 d\tau Z \infty dx e−sx − Z 1 d\tau Z \tau dx e−sx = Z \infty dx e−sx Z 1 d\tau − Z 1 dx e−sx Z 1 x
Z 1 dx e−sx Z x d\tau, (17.622) where we have made use of the fact that the sojourn-time density has its support on [0, 1]. Thus we can now invert the Laplace transform (17.620) of the sojourn-time density to obtain
\sqrt 2\pi Z x
p
(17.623) This integral can be performed analytically, with result59
r \pi (c −d)x + d(1 −x) p x(1 −x) exp c2 2 −(c −d)2 2x − d2 2(1 −x)
- 1 −(2d −c)2 e−2d(d−c) erfc
(c −d)(1 −x) + d x p 2x(1 −x) ! , (probability density for bridge sojourn time) (17.624) as can be verified by differentiating this expression to obtain the integrand of Eq. (17.623). This then gives the density for the boundary-crossing case of the bridge sojourn time. The cases with d < 0 can be generated from the cases with d \ge 0 by replacing x −\rightarrow 1 −x. And again, all of these expressions can be generalized to a bridge Bt(0\rightarrow c)(t) that is pinned to c at time T via the replacements c −\rightarrow c/ \sqrt t, d −\rightarrow d/ \sqrt t, and x −\rightarrow x/t. Overall factors of t must be restored to make the density come out with ‘‘units’’ of 1/t [including, e.g., an overall factor of \sqrt t that came from the path- pinning factor ec2/2 in Eq. (17.593)]. Also, a bridge starting at a different location a than 0 can be obtained by shifting c and d by −a. Writing these out explicitly from Eqs. (17.598) and (17.624),
h 1 −e−2(d−a)(d−c)/ti
r 2(t −x) \pit3x e(c−a)2/2t−(2d−a−c)2/2(t−x) + 1 t 1 −(2d −a −c)2 t e−2(d−a)(d−c)/t erfc s (2d −a −c)2x 2t(t −x) !
r \pi (c −d)x + (d −a)(t −x) p t3x(t −x) e(c−a)2/2t−(c−d)2/2x−(d−a)2/2(t−x)
- 1 t 1 −(2d −a −c)2 t e−2(d−a)(d−c)/t erfc
(c −d)(t −x) + (d −a) x p 2tx(t −x) ! , (probability density for bridge sojourn time) (17.625) 59cf. Andrei N. Borodin and Paavo Salminen, op. cit., p. 158, formula 1.4.8.
17.12 Sojourn Times while to obtain the other cases, we can change the signs of a, c, and d, while replacing x with t −x,
h 1 −e−2(a−d)(c−d)/ti
s 2x
- 1 t 1 −(a + c −2d)2 t e−2(a−d)(c−d)/t erfc r (a + c −2d)2(t −x) 2tx !
r \pi (a −d)x + (d −c)(t −x) p t3x(t −x) e(a−c)2/2t−(a−d)2/2x−(d−c)2/2(t−x)
- 1 t 1 −(a + c −2d)2 t e−2(a−d)(c−d)/t erfc
(a −d)(t −x) + (d −c) x p 2tx(t −x) ! , (probability density for bridge sojourn time, reflected cases) (17.626) which takes advantage of the fact that the problem is reflection symmetric if the occupation time is taken to be the non-occupation time. It is also useful to write out the moment-generating functions corresponding to these probability
replace that treatment with one that avoids obtaining a Bessel function, as we will do here. First, to set up the Laplace transforms, we will need to arrive at something of the form e−k \sqrt \lambda \sqrt \lambda \sqrt
\sqrt \lambda \sqrt
\sqrt \lambda ! = e−k \sqrt \lambda s \sqrt \lambda
p
= 2 s \sqrt \lambda − \sqrt
e−k \sqrt
\sqrt \lambda \sqrt \lambda . (17.627) Then using the Laplace-transform formula60 L \sqrt \pit3 e−st −1
\sqrt \lambda − \sqrt
(17.628) we can combine this with the transform formula (17.603) using the convolution theorem to obtain the formula L "Z t d\tau k 4\pi p
e−s(t−\tau) −1 #
\sqrt \lambda − \sqrt
e−k \sqrt \lambda. (17.629) Then adding this to the transform formula (17.520) to hit the last term on the right-hand side of Eq. (17.627), L " \sqrt
Z t d\tau k 2\pis p
e−s(t−\tau) −1 #
\sqrt \lambda \sqrt \lambda \sqrt
\sqrt \lambda \sqrt
\sqrt \lambda ! . (17.630) Changing k to k′ and combining this again with Eq. (17.520) gives L " \sqrt \pit
− Z t d\tau k′ 2\pis p
e−s(t−\tau) −1 # (\lambda) = e−k \sqrt \lambda \sqrt \lambda −e−k′\sqrt \lambda \sqrt \lambda \sqrt
\sqrt \lambda \sqrt
\sqrt \lambda ! . (17.631)
Chapter 17. Stochastic Processes
\sqrt 2 |c| and k′ = \sqrt 2(2d −c), we find61
exp −sTs[Bt(0\rightarrow c)]
\sqrt t(2d −c) \sqrt 2\pi s Z t d\tau p
1 −e−s(t−\tau) , (17.632) The non-integral part of the expression here is related to the boundary-touching probability of the bridge. Note that the normalization can be verified by taking the limit s −\rightarrow 0 and carrying out the resulting integral, recovering the unit-normalization result. Now putting this together with Eq. (17.619) and putting in a starting point of a, DD e−sTs EE
\sqrt t(2d −a −c) \sqrt 2\pi s Z t d\tau p
1 −e−s(t−\tau) DD e−sTs EE
\sqrt t \sqrt 2\pis Z t
p
DD e−sTs EE = 1 −e−2(a−d)(c−d)/t e−st
\sqrt t(a + c −2d) \sqrt 2\pi s Z t d\tau p
DD e−sTs EE
\sqrt t \sqrt 2\pis Z t
p
(sojourn-time moment generating function) (17.633) The last two cases here follow from the first two via the symmetry of the problem. Specifically, ‘‘flipping’’ the geometry by reversing the signs of a, c, and d allows us to calculate the ‘‘complimentary’’ generating function \langle \langle e−s(t−Ts)\rangle \rangle . Thus in addition to reversing these signs, we also need to reverse the sign of s in the first two expressions and then multiply through by e−st to obtain the desired generating function. In the fourth expression, we also changed integration variables, letting \tau −\rightarrow t −\tau, to bring it into a form more similar to the second expression. In this case it turns out to have the same form, but with a and c interchanged. The singularity of the integrand in Eqs. (17.633) can be problematic, and merits some further discus- sion. For example, in the case of the second expression with d = a, the integral naïvely becomes DD e−sTs EE
\sqrt t \sqrt 2\pis Z t d\tau (c −a) p
(17.634)
\tau = 0. Thus, some cases above may require some caution in their evaluation. One solution to regularize this example62 is to multiply through by s in Eq. (17.619), and then take the limit s −\rightarrow 0 to obtain ec2/2t \sqrt t \sqrt 2\pi Z t
p
(17.635) 61cf. Andrei N. Borodin and Paavo Salminen, op. cit., p. 158, formula 1.4.7. 62Vadim Linetsky, op. cit., in the discussion after Eq. (C.9).
17.12 Sojourn Times [That the result vanishes may be somewhat more obvious from performing the same maneuver in Eq. (17.618).] After restoring the initial point a, this is a useful counterpart to the second and fourth expressions in Eqs. (17.633). For example, multiplying Eq. (17.635) by e−st/s and subtracting from the second equation in Eqs. (17.633) leads to a similar expression, but with the replacement e−s\tau −\rightarrow e−s\tau −e−st. This cures the problematic case (17.634), which becomes DD e−sTs EE = \sqrt t \sqrt 2\pis Z t d\tau (c −a) p
e−s\tau −e−sti
(17.636) after this regularization; this cures the divergence at t = \tau, which is now cut off by the difference in exponentials. Note that this is the same expression that follows from the third expression in Eqs. (17.633).
the second equation in Eqs. (17.633); this again leads to a similar expression, but with the replacement e−s\tau −\rightarrow e−s\tau −1, cutting off the divergence at \tau = 0. However, the integrals are still pathological when
DD e−sTs EE = 1 −e−st st
(17.637) To see how this limit follows directly from Eq. (17.636), note that as c −\rightarrow a, the factor of (c −a) makes the integrand small everywhere except for the range of small \tau—the divergent factor \tau −3/2 is important here, and although it is cut off by the Gaussian factor e−(c−a)2/2\tau, it is still important for \tau ∼\sqrtc −a. Thus, since the significant part of the integrand moves towards vanishingly small \tau, we can expand the integrand in \tau to obtain DD e−sTs EE = \sqrt 2\pist Z t d\tau (c −a) \tau 3/2 e−(c−a)2(\tau −1−t−1)/21 −e−st + O
This integral may then be carried out, with result DD e−sTs EE
1 −e−st st erfc c −a \sqrt 2t
(17.639) At this point, the limit c −\rightarrow a yields Eq. (17.637). In general, the moments must then be computed by integrating the probability density or differentiating the moment-generating functions listed above. Although we wrote out an integral for c \le d and d \ge 0 in
work out a relatively nice expression by integrating the local time. The result is given in Eq. (17.724) as DD Ts[Ba\rightarrow c(t); d] EE = t 2 + sgn(2d −a −c) t h e−2[(d−a)(d−c) \Theta(d−a) \Theta(d−c)+(a−d)(c−d) \Theta(a−d) \Theta(c−d)]/t −1 i − r \pit 8 (2d −a −c) e(c−a)2/2t erfc |d −a| + |d −c| \sqrt 2t , (mean sojourn time) (17.640) after restoring the initial point a of the Brownian bridge. While it may not be obvious from this expression, when d is in the interval (a, c), the second term does not contribute, and the last term leads to a decreasing, straight-line dependence on d. 17.12.4 Path-Pinning Normalization as a Constrained Integration In deriving sojourn-time statistics, we made use of the relation (17.592) to remove the delta function that
\delta[W(T) −c] F[W(t)]
\sqrt 2\piT
F BT (0\rightarrow c)(t) . (17.641)
Chapter 17. Stochastic Processes We also used this in restricted form in Eqs. (17.561). In both cases we used probabilistic arguments to justify this relation. However, it is also useful to see how this arises as a constrained integration problem, by considering the explicit probability measure of the path. To begin, suppose that we consider N discrete
integration over the multidimensional Gaussian probability density:
\delta[W(T) −c] F[W(t)]
= Z d∆W0 . . . d∆WN−1
− 2∆T N−1 X j=0 ∆Wj 2 \delta N−1 X j=0 ∆Wj −c F[W(t)]. (17.642) Note that the path functional F[W(t)] is still written in continuous notation, since its details are not important to this calculation. Now we will define shifted increments that move linearly towards the ‘‘target’’ c,
N , (17.643) which will become the increments of the Brownian bridge. Changing variables, the integral becomes
\delta[W(T) −c] F[W(t)]
= Z d∆B0 . . . d∆BN−1
− 2∆T N−1 X j=0 ∆Bj −c N 2 \delta N−1 X j=0 ∆Bj F BT (0\rightarrow c)(t) = Z d∆B0 . . . d∆BN−1 e−c2/2T
− 2∆T N−1 X j=0 ∆Bj 2 + c T N−1 X j=0 ∆Bj
N−1 X j=0 ∆Bj F BT (0\rightarrow c)(t) . (17.644) Now carrying out the integration over ∆BN−1 will remove the delta function and enforce the replacement ∆BN−1 = − N−2 X j=0 ∆Bj. (17.645) However, some care must be taken in the resulting integral. The delta-function integration formula (17.65) reads Z
I h−1(0) f(q) |\nabla h| dS, (17.646) where the reduced integral involves the Euclidean norm |\nabla h| of the gradient of the constraint function h(q). In Eq. (17.644), the constraint function is a simple sum over all the ∆Bj, and the derivatives are taken with respect to each ∆Bj (i.e., each of the N derivatives in the gradient has unit magnitude). The corresponding norm |\nabla h| is then simply \sqrt N. Thus, the result of carrying out the ∆BN−1 integral in Eq. (17.644) gives
17.13 Local Time the desired result:¶
\delta[W(T) −c] F[W(t)]
= Z d∆B0 . . . d∆BN−2 e−c2/2T \sqrt
− 2∆T N−2 X j=0 ∆Bj 2 − 2∆T N−2 X j=0 ∆Bj 2 \times F BT (0\rightarrow c)(t)
\sqrt 2\piT Z d∆B0 . . . d∆BN−2
− 2∆T N−2 X j=0 ∆Bj 2 − 2∆T BN−1 2 \times F BT (0\rightarrow c)(t) .
\sqrt 2\piT
F BT (0\rightarrow c)(t) . (17.647) In the second-to-last expression here, there are N −1 steps of variance ∆T, in addition to a constraint that the variance of BN−1 is ∆T (BN−1 is close to zero rather than c because we defined the Bj to be
from the path measure, and part came from the integration over a delta function. 17.13 Local Time We will define the local time of a stochastic process y(t) at displacement d as
Z t dt′ \delta[y(t′) −d]. (17.648) (local time) The integrand here only ‘‘activates’’ when y(t) passes through d, and the local time is a measure of how much time y(t) spends at the displacement d, but normalized so that the answer is not merely zero. Recalling that we defined the sojourn time for y(t) as
Z t dt′ \Theta[y(t′) −d], (17.649) we can immediately deduce that
(17.650) (local time) so that we may simply adapt our sojourn-time results to obtain local-time statistics. One useful aspect of the local time arises in calculating functionals of stochastic processes of the form Z t
Z t dt′ Z
Z da F(a) Z t dt′ \delta[y(t′) −a]. (17.651) The definition of the local time then implies Z t
Z \infty −\infty da F(a) ℓ[y(t); a], (17.652) (local-time density formula) so that the local time acts as an occupation density for y(t). The local time is also commonly thought of as a time-dependent process, here through the time dependence of the process itself. Intuitively, the local time ‘‘accumulates’’ as the stochastic process continues, so ℓ[y(t); a] is a nondecreasing function of time.
Chapter 17. Stochastic Processes As an alternate representation of the local time, recall the property of the \delta function
X x0\in f −1(0) \delta(x −x0) |f ′(x0)| , (17.653) where the sum is over all real roots x0 of f(x). Then the definition (17.648) becomes
Z t dt′ X
\delta(t′ −td) |y′(td)| , (17.654) or carrying out the integral,
X
|y′(td)|. (intersection representation of local time) (17.655) Thus, the local time is given by summing over the intersections of the process y(t) with the boundary at d, where at each intersection the contribution is the reciprocal of the ‘‘speed’’ |y′| during the intersection. Intuitively, this makes sense, since the greater the speed, the less the time spent at the level d during the
form:
∆t\rightarrow 0 X {j:(Wj+1−d)(Wj−d)<0} ∆t |∆Wj|. (17.656) Note that as ∆t −\rightarrow 0, the contribution from each intersection in the sum decreases as ∆t1/2. As we will show, the local time can converge to a nonzero value; evidently, this means that the smaller step size ‘‘reveals’’ extra intersections in the neighborhood of each intersection to compensate for this decrease. It is also interesting to consider possible generalizations of the local time. For example, consider the functional
Z t dt′ \delta′[y(t′) −d], (17.657) which we can see is related to the local time via a −\partial d derivative, as the local time is related to the sojourn time. Using the composition rule (Problem 17.13)
X x0\in f −1(0) \delta′(x −x0) f ′(x0)|f ′(x0)| + f ′′(x0) \delta(x −x0) |f ′(x0)|3 , (17.658) the local-time derivative becomes
Z t dt′ X
\delta′(t′ −td) y′(td)|y′(td)| + y′′(td) \delta(t′ −td) |y′(td)|3 . (17.659) The first term under the integral vanishes so long as the velocity in the denominator does not vanish (something that occurs with zero probability, and which we also ignored in the local-time analysis), and we have
X
y′′(td) |y′(td)|3 . (17.660) (local-time derivative)
y′′(t) as O(∆t−3/2), so the summand here is of order unity. However, recall from our discussion of local-time intersections in Eq. (17.655) that the number of intersection times td grows as ∆t−1/2. Thus, while we can assign an ensemble average to this statistic, DD ℓ′[y(t); d] EE
DD ℓ[y(t); d] EE = \partial 2 d DD Ts[y(t); d] EE , (17.661) evidently the variance of the derivative statistic is arbitrarily large (see Problem 17.18).
17.13 Local Time 17.13.1 Wiener Process To compute the probability distribution for the local time of the Wiener process, we will follow closely the procedure of Section 17.12.1 for the sojourn time of the Wiener process. Let fℓ(x) denote the probability
(17.662) Then the Laplace transform of fℓ(x) is Z \infty
exp {−sℓ[W(t); d]}
=
exp −s Z t dt′ \delta[W(t′) −d] . (17.663) Note that fℓ(x) is not limited in domain to x < t as was the case for the sojourn time, but the domain is limited to x > 0. Consider then the driven diffusion equation
2\partial 2 x f −V (x)f −\lambdaf + g(x), (17.664) where V (x) is the occupation function, which here is a delta function:
(17.665)
dt exp −\lambdat − Z t dt′ V [x + W(t′)]
= Z \infty dt e−\lambdat ** exp −s Z t
. (17.666) This is then the solution of the steady-state version of Eq. (17.664):
2\partial 2 x f(x) −s\delta(x −d)f(x) + 1. (17.667) For x̸ = d, the ODE is
(17.668) Setting h = f −1/\lambda, we have h′′ = 2\lambdah, so that for x̸ = d, h(x) ∝e\pm \sqrt 2\lambda x, (17.669) or choosing the bounded solutions,
Ae− \sqrt
\lambda (x > d) Be \sqrt
\lambda (x < d) (17.670)
B = Ae−2 \sqrt 2\lambda d. (17.671)
− \sqrt 2\lambdaAe− \sqrt 2\lambda d −2s Ae− \sqrt
\lambda = \sqrt 2\lambdaBe \sqrt
\sqrt 2\lambdaAe− \sqrt 2\lambda d, (17.672)
Chapter 17. Stochastic Processes
as A = − se \sqrt 2\lambda d \lambda \sqrt
, B = − se− \sqrt 2\lambda d \lambda \sqrt
. (17.673)
Z \infty dt e−\lambdat ** exp −s Z t dt′ \delta[W(t′) −d]
\lambda = 1 \lambda − se− \sqrt 2\lambda d \lambda \sqrt
, (17.674) where we have assumed d > 0. Then using Eq. (17.663) on the left-hand side, Z \infty dt e−\lambdat
exp {−sℓ[W(t); d]}
= 1 \lambda − se− \sqrt 2\lambda d \lambda \sqrt
. (17.675) Now using the Laplace-transform formulae
\lambda (17.676) and63 L −eakea2t erfc a \sqrt t + k \sqrt t + erfc k \sqrt t
ae−k \sqrt \lambda \lambda \sqrt
, (17.677) which becomes with k = d \sqrt
\sqrt 2, L −esdes2t/2 erfc s p
d \sqrt 2t + erfc d \sqrt 2t
se−d \sqrt 2\lambda \lambda \sqrt
. (17.678) Thus, Eq. (17.675) becomes Z \infty
d \sqrt 2t + esdes2t/2 erfc s p
d \sqrt 2t (17.679) after using Eq. (17.663) to replace the ensemble average on the left-hand side. Now using Z \infty dx e−sx \sqrt
2 esd es2t/2 erfc "r t s + d t # , (17.680) we can invert the Laplace transforms on both sides, with the result64
|d| \sqrt 2t
r
(local-time probability density) (17.681) which is normalized in view of r
d \sqrt 2t = 1 −erf d \sqrt 2t . (17.682) 63Milton Abramowitz and Irene A. Stegun, op. cit., p. 1027, Eq. (29.3.89). 64cf. Andrei N. Borodin and Paavo Salminen, op. cit., p. 155, formula 1.3.4; A. N. Borodin, ‘‘Brownian local time,’’ Russian Mathematical Surveys 44, 1 (1989), p. 5 (doi: 10.1070/RM1989v044n02ABEH002050); Jim Pitman, ‘‘The distribution of local times of a Brownian bridge,’’ in Séminaire de Probabilités XXXIII, Jacques Azéma, Michel Émery, Michel Ledoux, and Marc Yor, Eds. (Springer, 1999), p. 388, Eq. (1) (doi: 10.1007/BFb0096528).
17.13 Local Time Again, the \delta-function term here gives the boundary-noncrossing probability, with the crossing probability (17.372) appearing, which is the same as the probability to have zero local time at d:
d \sqrt 2t = erf d \sqrt 2t . (17.683) We have also inserted absolute-value symbols for d: while we assumed d > 0 for this derivation, the result should be exactly symmetric in d owing to the same symmetry in W(t). The cumulative probability is again given by integrating 0 to x:65
Z x
|d| \sqrt 2t + erf x + |d| \sqrt 2t −erf |d| \sqrt 2t , (17.684) or
x + |d| \sqrt 2t (17.685) (local-time cumulative density) after cancelling terms. 17.13.2 Standard Brownian Bridge The local-time distributions for the standard Brownian bridge B(t) are special cases of the results in the following section, so we just quote them here.
h 1 −e−2|d|2i
x + 2|d|
(local-time probability density for B(t)) (17.686) is the probability density for the standard bridge, while
(local-time cumulative density for B(t)) (17.687) is the cumulative density. Note that as in Eq. (17.585), it is not difficult to compute moments of the local time using DD ℓnEE = n Z \infty
= n Z \infty dx xn−1e−(x+2|d|)2/2. (17.688) Thus, for example, DD ℓ EE = r\pi 2 erfc \sqrt 2 |d| (17.689) for the mean local time of a standard Brownian bridge. 17.13.3 Brownian Bridge
as in Section 17.13.1 up to Eq. (17.666), where we will leave g(x) in the Feynman–Kac formula:
Z \infty dt e−\lambdat ** g[x + W(t)] exp −s Z t
. (17.690) 65cf. Lajos Takács, ‘‘On the Local Time of the Brownian Motion,’’ The Annals of Applied Probability 5, 741 (1995), Eq. (3) (doi: 10.1214/aoap/1177004703).
Chapter 17. Stochastic Processes
steady-state diffusion equation
2\partial 2 x f(x) −s\delta(x −d)f(x) + eik(x−c). (17.691) For x̸ = d, the ODE is
(17.692)
−2eik(x−c) + 2k2eik(x−c)
(17.693) so that for x̸ = d, h(x) ∝e\pm \sqrt 2\lambda x, (17.694) or picking the bounded solutions in each domain,
Ae− \sqrt
(x > d) Be \sqrt
(x < d), (17.695)
B = Ae−2 \sqrt 2\lambda d. (17.696)
− \sqrt 2\lambdaAe− \sqrt 2\lambda d −2s Ae− \sqrt
= \sqrt 2\lambdaBe \sqrt
\sqrt 2\lambdaAe− \sqrt 2\lambda d, (17.697)
A = − se \sqrt 2\lambda d eik(d−c)
\sqrt
, B = − se− \sqrt 2\lambda d eik(d−c)
\sqrt
. (17.698)
Z \infty dt e−\lambdat ** eik[W (t)−c] exp −s Z t dt′ \delta[W(t′) −d]
e−ikc
(17.699) where we have again assumed d > 0 [but note that if we assume d < 0, the following treatment is the same, but with d −\rightarrow −d, so we will simply replace d by |d| in the factor exp(− \sqrt 2\lambda d)]. Integrating with respect to k and dividing through by 2\pi introduces the \delta function that pins the Wiener path to c at time t: Z \infty dt e−\lambdat ** \delta[W(t) −c] exp −s Z t dt′ \delta[W(t′) −d]
= 1 2\pi Z \infty −\infty dk B + e−ikc
. (17.700)
17.13 Local Time Using Eq. (17.592) to change the Wiener-path average into a bridge average, Z \infty dt e−c2/2t \sqrt t e−\lambdat ** exp −s Z t dt′ \delta[Bt(0\rightarrow c)(t′) −d]
= \sqrt 2\pi Z \infty −\infty dk B + e−ikc
. (17.701) Then using Eq. (17.663) on the left-hand side and carrying out the k integration, Z \infty dt e−c2/2t \sqrt t e−\lambdat
exp −sℓ[Bt(0\rightarrow c); d]
= \sqrt 2\pi Z \infty −\infty dk e−ikc
s e− \sqrt 2\lambda|d| eik(d−c)
\sqrt
= r\pi \lambda e− \sqrt 2\lambda|c| −
\sqrt
\sqrt \lambda \sqrt
. (17.702) Now using the Laplace-transform formulae66 L 1 \sqrt
\sqrt \lambda e−k \sqrt \lambda, (17.703) and67 L eakea2t erfc a \sqrt t + k \sqrt t
e−k \sqrt \lambda \sqrt \lambda \sqrt
, (17.704)
\sqrt
\sqrt 2, L es(d+|c−d|)es2t/2 erfc s p
\sqrt 2t
e−(d+|c−d|) \sqrt 2\lambda \sqrt \lambda \sqrt
\sqrt = \sqrt 2 e−(d+|c−d|) \sqrt 2\lambda \sqrt \lambda \sqrt
. (17.705) Thus, Eq. (17.701) becomes Z \infty
r \pit 2 s ec2/2tes(|d|+|c−d|)es2t/2 erfc
s r t
\sqrt 2t ! (17.706) after using Eq. (17.663) to replace the ensemble average on the left-hand side. Now again using Z \infty dx e−sx \sqrt
2 esd es2t/2 erfc "r t s + d t # , (17.707) with the derivative rule
(17.708) so that − Z \infty dx e−sx (x + d) \sqrt
\sqrt
2 esd es2t/2 erfc "r t s + d t # , (17.709) we can invert the Laplace transforms on both sides, with the result68
h 1 −e[c2−(|d|+|c−d|)2]/2ti
t
(local-time probability density for Bt(0\rightarrow c)(t′)) (17.710) 66Milton Abramowitz and Irene A. Stegun, op. cit., p. 1026, Eq. (29.3.84). 67Milton Abramowitz and Irene A. Stegun, op. cit., p. 1027, Eq. (29.3.90). 68cf. Andrei N. Borodin and Paavo Salminen, op. cit., p. 155, formula 1.3.8; A. N. Borodin, ‘‘Brownian local time,’’ Russian Mathematical Surveys 44, 1 (1989), p. 6 (doi: 10.1070/RM1989v044n02ABEH002050).
Chapter 17. Stochastic Processes The \delta-function term here gives the boundary-noncrossing probability in the case of c < d, where the expo- nential part reduces to exp[−2d(d −c)], in agreement with Eq. (17.387). The corresponding cumulative probability is
Z x
h
0 , (17.711) or
(local-time cumulative density for Bt(0\rightarrow c)) (17.712) after cancelling terms. 17.13.3.1 Moments Note that as in Eqs. (17.585) and (17.688), we can compute moments of the local time via the cumulative as probability function DD ℓnEE = n Z \infty
= n Z \infty dx xn−1 e[c2−(x+|d|+|c−d|)2]/2t (17.713) for n > 0. For example, the mean is
ℓ
= r \pit 2 ec2/2t erfc |d| + |c −d| \sqrt 2t . (17.714)
the argument |d| + |c −d| reduces to |c|, in which case the mean local time is independent of the boundary location, which seems peculiar.69 Since the probability (17.712) and density (17.710) depend on c and d in exactly the same way, these probabilities are also invariant to shifts of the interface, provided the shift keeps the interface d between 0 and c.70 17.13.3.2 Moment-Generating Function Since the moments are fairly straightforward to calculate, it shouldn’t come as much surprise that the moment-generating function is easy to obtain. Writing this out, DD e−sℓEE = Z \infty dx e−sx fℓ(x) = s Z \infty
= 1 −s Z \infty dx e−sx
= 1 −s Z \infty dx e−sxe[c2−(x+|d|+|c−d|)2]/2t, (17.715) where we used Eq. (17.712) for the cumulative density. Note that in the first step, we discarded the boundary
69Zhiyi Chi, Vladimir Pozdnyakov, and Jun Yan, ‘‘On expected occupation time of Brownian bridge,’’ Statistics and Probability Letters 97, 83 (2015) (doi: 10.1016/j.spl.2014.11.009). 70Peter Howard and Kevin Zumbrun, ‘‘Shift invariance of the occupation time of the Brownian bridge process,’’ Statistics and Probability Letters 45, 379 (1999) (doi: 10.1016/S0167-7152(99)00080-2).
17.13 Local Time¶
component separately). Carrying out the final integral, we find71 DD e−sℓEE = 1 −s r \pit
\sqrt 2t (local-time moment-generating function for Bt(0\rightarrow c)) (17.716) for the moment-generating function. 17.13.3.3 Application to the Sojourn Time The local-time mean can also be used to obtain an expression for the mean of the sojourn time for a general Brownian bridge. By taking the expectation value of Eq. (17.650), we can simply write the mean sojourn time as the integral
Ts[y(t); d]
= Z \infty d dx
ℓ[y(t); x]
. (17.717) To handle this integral, we can use the integral formula Z
a −
x −b a erfc(ax −b), (17.718)
Ts[y(t); d]
= r \pit 2 ec2/2t Z \infty d dx erfc 2x −c \sqrt 2t = t 2e−2d(d−c)/t − r \pit 8 ec2/2t(2d −c) erfc 2d −c \sqrt 2t . (17.719) In the case d \ge 0 and c \ge d, we can divide the integrand into the parts with x \le c and x \ge c, the latter of which we just did:
Ts[y(t); d]
= r \pit 2 ec2/2t Z c d dx erfc |x| \sqrt 2t +
Ts[y(t); c]
= r \pit 2 ec2/2t(c −d) erfc c \sqrt 2t + t 2 − r \pit 8 ec2/2tc erfc c \sqrt 2t = t 2 − r \pit 8 ec2/2t(2d −c) erfc c \sqrt 2t . (17.720)
Ts[y(t); d]
= t 2 e−2[d(d−c)\Theta(d)\Theta(d−c)]/t − r \pit 8 (2d −c) ec2/2t erfc |d| + |c −d| \sqrt 2t . (17.721) For the case d \le 0, we can change the signs of both c and d to obtain the mirror image, and replace the sojourn time by the sojourn time subtracted from t. The net result for d \le 0 is
Ts[y(t); −d]
= t −t 2 e−2[d(d−c)\Theta(−d)\Theta(c−d)]/t − r \pit 8 (2d −c) ec2/2t erfc |d| + |c −d| \sqrt 2t . (17.722) Equations (17.721) and (17.722) can then be combined into a single case as
Ts[y(t); d]
= t 2 + sgn(d) t h e−2[d(d−c)\Theta(d)\Theta(d−c)+d(d−c)\Theta(−d)\Theta(c−d)]/t −1 i − r \pit 8 (2d −c) ec2/2t erfc |d| + |c −d| \sqrt 2t . (17.723) 71cf. Andrei N. Borodin and Paavo Salminen, op. cit., p. 155, formula 1.3.7.
Chapter 17. Stochastic Processes
Ts[y(t); d]
= t 2 + sgn(2d −c) t h e−2[d(d−c)\Theta(d)\Theta(d−c)+d(d−c)\Theta(−d)\Theta(c−d)]/t −1 i − r \pit 8 (2d −c) ec2/2t erfc |d| + |c −d| \sqrt 2t , (17.724) which makes the expression symmetric about c/2. 17.13.4 Local Time and Discontinuities in Stochastic Processes Because of the relation of the local time to the delta function, the local time tends to show up in the theory of stochastic processes especially when discontinuities are involved. We will consider two examples of this here: the absolute-value process, and a diffusion process with a discontinuity in the diffusion rate. 17.13.4.1 Reflected Brownian Motion
This is an example of reflected Brownian motion, in the sense that this is Brownian motion on (0, \infty), but when the process encounters the origin and attempts to cross it, it is ‘‘reflected’’ back into the positive real axis. The It¯o differential for this process is, expanding to second order using the It¯o chain rule (17.193),
dW dW + 1 d2|W| dW 2 dW 2
(17.725) where we have used d|x|
d|x|2
(17.726) thinking of the signum function as twice the Heaviside function plus a constant offset. Then integrating this transformed SDE from 0 to t,
Z t sgn[W(t′)] dW(t′) + ℓ[W(t); 0], (17.727) (Tanaka formula) where we used the local-time definition (17.648). This is the Tanaka formula.72 This can also be written as
Z t sgn[W(t′) −d] dW(t′) + ℓ[W(t); d], (17.728) (Tanaka formula) if W(t) is shifted by d at the start of the derivation. This says that in It¯o calculus, the reflecting effect of the absolute value is reflected in a deterministic ‘‘pushing’’ term that activates whenever the process hits the origin. 17.13.4.2 Discontinuous Diffusion Rate Consider the It¯o-form SDE
(17.729) representing state-dependent diffusion without an explicit drift term. We want to consider the case where \beta(y) has a discontinuity (or a number of isolated discontinuities). As a particular model, consider
\beta>, y > 0 \beta<, y < 0, (17.730) 72Andrei N. Borodin and Paavo Salminen, op. cit., p. 43.
17.13 Local Time with the point \beta(0) being, say, the average of the left- and right-hand values. As we argued before in Section 17.3.5.1, the probability density of particles governed by this SDE evolves by the Fokker–Planck equation
2\partial 2 y \beta2(y)P(y, t), (17.731)
Making the Lamperti transform (Section 17.7.4.2), we construct the function
Z y dy′
y > 0
y < 0 (17.732) (which is invertible provided \beta≷have the same sign), we can transform our SDE to the form of Eq. (17.359), dz = −1
(17.733) Evaluating the derivative of \beta as the derivative of a step function, dz = −1
(17.734) Now we have
\delta(y)
2\delta(y)
, (17.735) where we think of the delta function as a limit of centered distributions, so that the derivative of S(y) is the average of the derivatives on either side of the discontinuity. Thus, we have the fully transformed SDE dz =
(17.736) In this form, where we have scaled away the variance, the discontinuity arises as a deterministic ‘‘kick’’ term that acts only at the discontinuity’s location. Using the local-time definition (17.648), we can integrate this to find
ℓ[W(t); 0], (17.737) (skew Brownian motion) so that z(t) has the form of a regular Wiener path, with a shift that depends on the time the particle spent at the discontinuity (and thus being kicked by the local gradient). The motion of z(t) is called skew Brownian motion. 17.13.4.3 Skew Brownian Motion Skew Brownian motion73 is equivalent to ordinary Brownian motion, except that the probability is skewed towards one direction at the origin. Suppose we take a random walk of steps taken every ∆t, and probability density
x − \sqrt ∆t
x + \sqrt ∆t , (17.738) so that the probability of a step of size x = + \sqrt ∆t is p, and the probability is 1−p of stepping in the opposite direction. The second moment of each step is Z \infty −\infty
(17.739) 73Skew Brownian motion was introduced by Kiyosi Itô and Henry P. McKean, Jr., Diffusion Processes and their Sample Paths (Springer-Verlag, 1974). See also J. M. Harrison and L. A. Shepp, ‘‘On Skew Brownian Motion,’’ The Annals of Probability 9, 309 (1981) (doi: 10.1214/aop/1176994472); and Antoine Lejay, ‘‘On the constructions of the skew Brownian motion,’’ Probability Surveys 3 413 (2006) (doi: 10.1214/154957807000000013).
Chapter 17. Stochastic Processes
gives a bias of Z \infty −\infty
\sqrt ∆t −(1 −p) \sqrt
\sqrt ∆t (17.740)
an inhomogeneous walk where we take a biased stepping probability p in the region x \in [− \sqrt ∆t/2, \sqrt ∆t/2), and an unbiased probability p = 1/2 outside this region, then we have the discrete-step equation
x \sqrt ∆t (2p −1) \sqrt ∆t, (17.741)
stochastic component of the random walk by an equivalent Gaussian process in view of the limit ∆t −\rightarrow 0. Note that this process has the correct mean
\sqrt ∆t/2, \sqrt ∆t/2)(x), (17.742) where 1A(x) is the indicator function for x \in A, and the process has the correct second moment,
∆x2 =
∆W 2 = ∆t, (17.743) provided we ignore terms of order ∆W ∆t and ∆t2. Then as ∆t −\rightarrow 0, we can also note that \sqrt ∆t rect x \sqrt ∆t −\rightarrow \delta(x), (17.744) so that in this limit we obtain the SDE
(17.745) (skew Brownian motion) with solution
(17.746) (skew Brownian motion) These equations are equivalent to Eqs. (17.736) and (17.737) if we identify p = \beta<
, 1 −p = \beta>
. (17.747) As an example at the extreme end of this range of probabilities, note that the reflected Brownian motion in
the sgn into dW, which is possible due to the symmetry of dW). In this limit, the process is always reflected upward at the origin. In terms of the inhomogeneous diffusion, this corresponds to the limit where \beta< ≫\beta>, though the diffusion picture does not necessarily show the same repulsion from the y < 0 region; this is more an artifact of the transformation (17.732), which ‘‘smooshes’’ the y < 0 region up against the origin if \beta < is large. 17.13.4.4 Reflections, Images, and Transition Densities for Skew Brownian Motion The original construction of skew Brownian motion74 began with reflected Brownian motion, and then introduced probabilities p and 1−p of reflecting |W(t)| upwards or downwards, respectively, on each encounter with the origin. In this sense,
(17.748) (‘‘reflection coefficient’’ for diffusion) 74Kiyosi Itô and Henry P. McKean, Jr., op. cit.
17.13 Local Time acts as a ‘‘reflection coefficient’’ from above the interface for the Wiener path. The reflection occurs always
a refractive interface, the reflection coefficient r has the form of the Fresnel coefficient for the electric field
Fresnel transmission coefficient from above]. Indeed, the reflection principle can give the transition density for the trajectories,75 making the con- nection to a reflection coefficient more explicit. To compute the transition densities for skew Brownian motion, we will use a procedure similar to the calculation for boundary crossings of Brownian bridges in Section 17.8.2. To begin, let’s compute the probability density for a skew-Brownian-motion process z(t) starting at x0 > 0 and ending at x > 0. Let \tau0 denote the first-crossing time through the boundary x = 0. Then we start with
(17.749) That is, the path either hits the origin or it doesn’t before time t. We can compute the first term using the Reflection Principle, similar to the calculation in Section 17.8:
(17.750) where \phi(x) is the transition density for W(t) from 0,
\sqrt
(17.751) In changing from skew Brownian to regular Brownian motion, we used the ratio of probabilities for taking
= p
(17.752)
We can evaluate second term in Eq. (17.749) using the method of images, which works in the same way as the method of images for potentials due to charges in the presence of boundary conditions. Here,
We can achieve this by considering the density 75J. B. Walsh, ‘‘A diffusion with a discontinuous local time,’’ Société Mathématique de France Astérisque 52-53, 37 (1978); J. F. Le Gall, ‘‘One-dimensional stochastic differential equations involving the local times of the unknown process,’’ in Stochastic Analysis and Applications: Proceedings of the International Conference held in Swansea, April 11-15, 1983, A. Truman and D. Williams, Eds. (Springer-Verlag, 1984), p. 51; G. J. M. Uffink, ‘‘A random walk method for the simulation of macrodispersion in a stratified aquifer,’’ in Relation of Groundwater Quantity and Quality (Proceedings of the Hamburg Symposium, August 1983) (IAHS Publ. no. 146, 1985), p. 103; Antoine Lejay, op. cit., around Eq. (42); Thilanka Appuhamillage, Vrushali Bokil, Enrique Thomann, Edward Waymire, and Brian Wood, ‘‘Occupation and local times for skew Brownian motion with applications to dispersion across an interface,’’ Annals of Applied Probability 21, 183 (2011), Eq. (2.2) (doi: 10.1214/10-AAP691); Eric M. LaBolle and Graham E. Fogg, ‘‘Reply [to ‘‘Comment on ‘Diffusion theory for transport in porous media: Transition-probability densities of diffusion processes corresponding to advection-dispersion equations’ by Eric M. LaBolle et al.’’],’’ Water Resources Research 36, 823 (2012) (doi: 10.1029/1999WR900325).
Chapter 17. Stochastic Processes diffusing from a distribution initially localized at x0, balanced by a distribution of negative density at −x0, such that perfect cancellation occurs at x = 0:
(17.753) Putting this all together in Eq. (17.749), we then have
(17.754) again for x > 0 and x0 > 0. For x0 > 0 and x < 0, the second term in Eq. (17.749) corresponds to an impossible event, and by a similar Reflection-Principle argument, the second term is 2(1 −p) \phi(x + x0):
(17.755)
z > 0, z0 > 0
z < 0, z0 > 0, (transition density for skew Brownian motion) (17.756) with \phi(x) defined in Eq. (17.751), and the reflection coefficient r given by Eq. (17.748). The same formula applies for z0 < 0 with minor changes,
z < 0, z0 < 0
z > 0, z0 < 0, (transition density for skew Brownian motion) (17.757) where both r and z change sign. If we translate these results back into the language of the discontinuous-diffusion SDE (17.729), we can invert the Lamperti transform (17.732) to obtain
z\beta>, z > 0 −z\beta<, z < 0. (17.758) Then transforming the probability density, Eq. (17.756) becomes
\phi y −y0 \beta> dy \beta>
\beta> , y > 0, y0 > 0 (1 −r) \phi
\beta< dy \beta< , y < 0, y0 > 0, (transition density for discontinuous diffusion) (17.759) where note that in the y < 0 case, the y0 was determined by x0 in the Reflection-Principle argument for x > 0, and should correspond to the same distance from the origin as y0/\beta> to properly match the boundary
at the same time only with this rescaling). Again, the y0 < 0 case can be obtained by exchanging \beta≷ everywhere, including in the reflection coefficient r. 17.13.4.5 Stratonovich Discontinuous Diffusion We began the discussion of a discontinuous diffusion rate with the It¯o-form diffusion (17.729), and the choice of It¯o calculus was important in generating the ‘‘pumping’’ effect at the discontinuity. If we instead consider a Stratonovich-form diffusion SDE
(17.760) corresponding to the It¯o SDE dy = 1
(17.761)
17.14 Bessel Processes then things are different; in fact, the Lamperti-transformed version of this SDE is simply dz = dW, (17.762) as we discussed in Section 17.7.4.2. The equivalent Fokker–Planck equation from Eq. (17.130) is
2\partial 2 y \beta2(y)P(y, t), (17.763) or commuting a derivative,
(equivalent Fokker–Planck equation) (17.764) Note that by changing variables to z according to dz dy = 1 \beta , (17.765) which is the same as the transformation (17.732), we can rewrite the Fokker–Planck equation as a standard diffusion equation,
2\partial 2 z \rho, (17.766) (equivalent Fokker–Planck equation) where
(17.767) This simplification is the same as the result dz = dW from the Lamperti transformation of the Stratonovich- diffusion SDE. 17.14 Bessel Processes A Bessel process in d dimensions is simply the ‘‘radius’’ or Euclidian norm of a vector Wiener process in d dimensions: Rd :=
W(t)
d = v u u t d X \alpha=1 [W\alpha(t)]2 . (17.768) (Bessel process) As we will see, this process is useful for characterizing the motion of a vector Wiener process, particularly in the sense of whether or not it returns to its starting point. Also, note that this is the multidimensional generalization of the reflecting Wiener process of Section 17.13.4.1, to which this reduces for d = 1. To begin, we can try to work out an SDE for Rd. First, writing out the differential for Rd, dRd = 2Rd d X \alpha=1
\alpha = d 2Rd dt + 1 Rd d X \alpha=1
(17.769) Unfortunately, it is difficult to proceed directly from this point. However, we can also work more simply with the square of the Bessel process, and write out the differential: dR 2 d = d X \alpha=1
\alpha
d X \alpha=1
(17.770) On the other hand, we can write out the differential as dR 2
(17.771)
Chapter 17. Stochastic Processes From Eq. (17.769), we can find
R 2 d d X \alpha=1 W 2
(17.772)
in Eq. (17.771), we find dRd = d −1 2Rd dt + 1 Rd d X \alpha=1
(17.773) What remains is to analyze the last term, which turns out to be equivalent to a Wiener process. Defining
Rd d X \alpha=1
(17.774)
d ˜W 2 = R 2 d d X \alpha=1 W 2
(17.775) The higher-order moments of d ˜W vanish also, so we can see that ˜W is another representation of a Wiener process. Then dropping the tilde, we have the SDE dRd = d −1 2Rd dt + dW (17.776) (Bessel-process SDE) to directly generate the dynamics of the Bessel process. Note that because of the divergence at Rd = 0, we should be careful to specify that the transformation to this equation is only valid provided Rd̸ = 0. That is, if the Bessel process hits the origin, then this equation becomes invalid afterwards. [For example, in d = 1, we are missing the delta-function term that arises in Eq. (17.725).] We will soon return to the question of whether or not this actually happens. 17.14.1 Generator The term ‘‘Bessel process’’ comes specifically for the generator. Recall that an SDE of the form dx = \alpha dt + \beta dW implies a diffusion-type equation, which can be summarized by the generator [see Eq. (17.136)]
\partial 2 x . (17.777) For the Bessel-process SDE (17.776), the generator is G = d −1
2\partial 2 x , (17.778)
d −1 2x u′ + 1
(17.779) [Note that this can be viewed as a Laplace transform of the backward Kolmogorov equation (17.144).] The idea is that this differential equation is equivalent to a Bessel equation, as we will now see. First, changing variables according to v = xmu, (17.780)
17.14 Bessel Processes¶
2x−mv′′ + d −2m −1 x−m−1v′ −m(d −m −2)
(17.781) Then multiplying through by 2xm+2, x2v′′ + (d −2m −1)xv′ −
v = 0. (17.782)
˜x2v′′ + ˜xv′ − ˜x2 + m2 v = 0. (17.783) This is the modified Bessel equation, and the solutions are the modified Bessel functions Im(˜x) and Km(˜x).
ez/\sqrtz, while Km(z) −\rightarrow 0. So the Km(z) represent physical solutions in this problem. 17.14.2 Brownian Recurrence One useful question we can address with Bessel processes is, given a Wiener path in d dimensions starting at the origin and running for arbitrarily large times, how far away does it go, and does it return to the origin? 17.14.2.1 One Dimension In d = 1, we have already answered this question by analyzing boundary-crossing probabilities. That is, from Eq. (17.372) we know that as t −\rightarrow \infty, the probability for crossing a boundary at distance d away converges to unity for any value of d. This means that the path wanders arbitrarily far away from the origin, and once it has wandered any given distance away, it is guaranteed to return to the origin. Thus Brownian motion in d = 1 dimension is said to be recurrent. 17.14.2.2 Two Dimensions The d = 2 case turns out to be marginally recurrent, as we will see. We already know that R2 will wander arbitrarily far from the origin, because the projection onto one dimension already does this, as we know from the d = 1 case. However, the question of returning to the origin is somewhat different. To analyze this case, it is convenient to define L := log R2. (17.784) Then the SDE (17.776) for this case, dR2 = 2R2 dt + dW, (17.785) will simplify as follows. Expanding the logarithm to second order gives dL = dR2 R2 −(dR2)2 2R 2 , (17.786) which with Eq. (17.785) becomes dL = 2R 2 dt + 1 R2 dW −dW 2 2R 2 = dW R2 . (17.787) Then since R2 = eL, dL = e−L dW. (17.788)
Chapter 17. Stochastic Processes This transformation simplifies things by removing the drift, at the expense of multiplicative noise. However, we will use the martingale property of this noise, which means that the ensemble average vanishes:
(17.789) This means that \langle \langle L\rangle \rangle = \langle \langle log R2\rangle \rangle is a constant of the motion. Since we know that L will take on arbitrarily large values, this essentially means that L will also have to take on arbitrarily negative values to keep the mean constant. This means that R2 will come arbitrarily close to 0. To make this argument more formal,76 consider two positive radii R< and R> that bound the initial value R0 of R2: R< \le R0 \le R>. Since \langle \langle L\rangle \rangle is constant, log R0 =
log R2(t)
. (17.790) The reasoning that led to Eq. (17.433) applies in a similar way here. Let \tau be the first time at which R2(t) achieves either R< or R>. Then we can consider a paths up to its stopping time \tau, at which point L takes the value log R2(\tau); when averaging over all continuations of the path up to time t > \tau, the value log R(\tau) is unchanged. Now extending the average over all paths up to the stopping time, we see that it is equivalent to write Eq. (17.790) in terms of the hitting time as log R0 =
log R2(\tau)
(17.791) after taking the limit t −\rightarrow \infty. Now writing the expectation value in terms of probabilities for first hitting R< or R>,
(17.792)
log R> . (17.793) Now fixing R0 and R< \le R0, we can let R> −\rightarrow \inftyto obtain
(17.794) That is, the probability for R2(t) to ‘‘hit infinity before hitting R<’’ (which means never hitting R<) is zero. Thus, for any inner bound R> > 0, R2(t) is guaranteed to hit it in finite time. However, the same argument with fixed R> and R< −\rightarrow 0 gives that R2(t) will never hit the origin in finite time. Thus, so to speak, the origin is nonrecurrent, but any disc surrounding the origin is recurrent. 17.14.2.3 Three and More Dimensions In three and more dimensions, it turns out that Brownian motion is nonrecurrent, and in fact the path wanders ‘‘far away’’ in the sense that Rd(t) −\rightarrow \infty. It is sufficient to show this in three dimensions, again because any projection of a d-dimensional motion onto three dimensions will satisfy these conditions. However, it isn’t much more difficult to treat any d \ge 3 directly. In this case, the simplifying transformation for the SDE (17.776) has the form P := aR b d , (17.795) for constants a and b to be determined. Computing the differential of P (and expanding to second order), dP = abR b−1 d dRd + ab(b −1) R b−2 d (dRd)2, (17.796) 76See Ioannis Karatzas and Steven E. Shreve, Brownian Motion and Stochastic Calculus (Springer, 1988), pp. 161-3 (doi:
17.14 Bessel Processes and then using Eq. (17.776) for dRd,¶
R b−2 d dt + abR b−1 d dW. (17.797) We can then force the drift term to vanish by setting b = 2 −d, (17.798) in which case
d dW. (17.799) Now since P = aR2−d d
(17.800) where c := 1 −d 2 −d = d −1 d −2. (17.801) Finally, if we choose
(17.802) we end up with the SDE dP = −P c dW, (17.803) which again has the martingale property \langle \langle dP\rangle \rangle = 0. In what follows, it is useful to remember that for d > 2, P goes as an inverse power of Rd.
outer boundaries satisfying R< \le R0 \le R>, with \tau the first passage time through either boundary. Then since P0 =
P(\tau)
, (17.804) this becomes R2−d = DD R2−d d (\tau) EE (17.805) directly in terms of the radial coordinate. We can rewrite the expectation value as before in terms of the probabilities for first crossing the inner and outer boundaries as R2−d = R2−d <
(17.806) Now if R> −\rightarrow \infty, the last term vanishes, and
R< R0 d−2 . (probability for achieving inner bound R<) (17.807) This gives the probability for the path, starting from R0 > 0, to hit R< \le R0. Note that this probability converges to zero as R< −\rightarrow 0. Thus, the probability for a path starting at R0 > 0 to hit the origin is zero (and thus the origin is nonrecurrent, since once a path moves away from the origin, it will not return). Each path starting at R0 has a (pathwise) minimum radius, and the probability density for this minimum is given by differentiating (17.807), with the result
R d−2 . (probability density for minimum radius) (17.808) Now to show that Rd −\rightarrow \inftyfor d \ge 3, first fix some upper limit R>. We want to show that after some time
that the path does not return below R>.
Chapter 17. Stochastic Processes 17.15 Differentiation of Stochastic Quantities We have so far studied stochastic processes and ensemble averages over stochastic problems. In cases where we compute quantities like the sojourn time or the local time, analytic expressions are available for the ensemble average, generating function, and so on. But for more complicated statistics, analytic expressions may not be available, and numerical simulations are needed to compute such averages. Another important class of ensemble-average problems that we have not yet discussed is the differen- tiation of ensemble-averaged statistics. Since Wiener processes are nondifferentiable, the differentation of numeric quantities with underlying Wiener processes can be tricky. Recall, for example, that the sojourn time can be differentiated to yield the local time [see Eq. (17.650)]; this differentiation then carries over to the mean sojourn and local times, but differentiation of the underlying Wiener process causes larger fluctu- ations in the ensemble-average computation on a pathwise basis [see Eq. (17.656), due to the variation in the number of intersections and the possibility of small values of ∆Wj]. In a further differentiation of the local-time process, the derivative of the ensemble average is well-defined, but the pathwise derivative has arbitrarily large fluctuations in the continuum limit [see Eq. (17.660) and the following discussion]. Thus, it is useful to consider some general approaches to handling derivatives of stochastic processes, which can be handled efficiently in certain cases with a variety of tricks. As a byproduct, we will also briefly study the derivative of a stochastic process with respect to a stochastic quantity, which gives the generalizes variational calculus to stochastic processes. Such derivatives are useful in, for example, computing Casimir forces from energy path integrals (Chapter 21). These techniques are also widely employed in financial mathematics, where ensemble averages over stochastic trajectories are used to price financial derivatives (e.g., financial options). In such pricing, derivatives of the price with respect to parameters such as the volatility or starting price yields the ‘‘sensitivity’’ of the price with respect to parameter variations. These techniques also apply in optimization problems, where the ‘‘payoff’’ to optimize is estimated via Monte-Carlo simulations. 17.15.1 Parameter Differentiation and Likelihood Estimation Consider an ensemble average of the form
DD \Phi(Z) EE Z, (17.809) (model ensemble average)
quantities that we wish to average, and the \lambda on the left-hand side indicates a parameter dependence of the ensemble average. In some important cases we can arrange things such that the only dependence on the right-hand side on the parameter \lambda is in the probability distribution for Z. That is, we can write the ensemble average as
Z dz \Phi(z) f\lambda(z), (17.810) where f\lambda(z) is the probability density for Z, and the \lambda subscript emphasizes the dependence on \lambda. This approach can work even in some cases where it seems that the payoff function should represent the parameter dependence. For example, suppose we differentiate the sojourn time of Wiener paths with respect to the boundary distance d, which is a part of the functional rather than the paths. However, this quantity can be computed equivalently by differentiating with respect to the starting point of the paths. Now the derivative can be written simply as
Z
(17.811) and by multiplying and dividing the integrand by f\lambda(z), we can rewrite this as an ensemble average:
f\lambda(Z)
Z =
Z . (likelihood-ratio derivative estimator) (17.812)
17.15 Differentiation of Stochastic Quantities Of course, this generalizes readily to higher derivatives, \partial m
\Phi(Z) \partial m
f\lambda(Z)
Z , (likelihood ratio mth derivative estimator) (17.813) although of course without the logarithmic form. Thus, the derivative here appears simply as a reweighted version of the original ensemble average (17.809). Of course, ‘‘good’’ behavior of this ensemble average is not guaranteed, but the hope is that the logarithmic weight will not cause much in the way of fluctuations. In particular, this approach should give an advantage when the probability density is smooth, whereas the payoff function is not (e.g., it may have a discontinuity or singularity as in the sojourn or local time, such that a small change in the boundary location d can produce a large change in the payoff value). One particularly nice property of this expression is that the weight is universal for any payoff function \Phi(Z), because the parameter dependence lies entirely with the probability distribution. 17.15.1.1 Likelihood Now for a brief digression to explain some terminology. The likelihood of a parameter \lambda, given some observed outcome Z, is defined as the probability of the outcome given the particular parameter value:
(17.814) In our ensemble average, the likelihood is just the probability density f\lambda(Z). Then the ratio of two like- lihoods f\lambda′(Z)/f\lambda(Z) is a common quantity in statistics to evaluate the relative plausibility of two models given observed data. Then the likelihood-ratio derivative \partial \lambda′[f\lambda′(Z)/f\lambda(Z)] appears in the ensemble average (17.812), so a Monte-Carlo evaluation of Eq. (17.812) is known as a likelihood-ratio estimator for the derivative.77 17.15.1.2 Application: Differentiation of Brownian-Bridge Path Averages When computing statistics related to sojourn or local times, or also boundary-crossing and escape probabil- ities, typically path-averages functional of the form
DD \Phi[V (x)] EE x(t), (17.815) (model ensemble average)
average. As noted before, derivatives with respect to a boundary position d can be regarded as derivatives with respect to the initial point x0 of the path. Thus, we will specifically consider Brownian-bridge paths
(17.816)
if we consider the time-sliced path in N steps of duration ∆t = t/N, the probability density of the paths in the average (17.815) have the x0-dependent factor f[x(t)] ∝e−(x1−x0)2/2∆t e−(x0−xN−1)2/2∆t, (17.817) 77Martin I. Reiman and Alan Weiss, ‘‘Sensitivity analysis via likelihood ratios,’’ Proceedings of the 1986 Winter Simulation Conference, J. Wilson, J. Henriksen, and S. Roberts, Eds., p. 285 (1986) (doi: 10.1145/318242.318450). P. W. Glynn, ‘‘Stochastic approximation for Monte Carlo optimization,’’ Proceedings of the 1986 Winter Simulation Conference, J. Wilson, J. Henriksen, and S. Roberts, Eds., p. 356 (1986) (doi: 10.1145/318242.318459). Reuven Y. Rubinstein, ‘‘Sensitivity Analysis and Performance Extrapolation for Computer Simulation Models,’’ Operations Research, 37, 72 (1989) (doi: 10.1287/opre.37.1.72). P. W. Glynn, ‘‘Likelihood Ratio Derivative Estimators for Stochastic Systems,’’ Proceedings of the 1989 Winter Simulation Conference, E. A. MacNair, K. J. Musselman, and P. Heidelberger, Eds., p. 374 (1989) (doi: 10.1109/WSC.1989.718702).
Chapter 17. Stochastic Processes
\partial x0f[x(t)] f[x(t)]
∆t −(x0 −xN−1) ∆t
∆t , (17.818) where we are defining the shorthand
. (17.819) Thus, the first derivative of the functional (17.815) becomes
** \Phi[V (x)] 2(¯x1 −x0) ∆t ++ x(t) . (differentiated ensemble average) (17.820)
It is also straightforward to derive an expression for high-order derivatives. First, rewriting Eq. (17.817) as f[x(t)] ∝e−(x0−¯x1)2/∆t e−(xN−1−x1)2/4∆t (17.821)
\sqrt ∆t, and use the Hermite-polynomial definition
y e−y2 (17.822) we can write out the derivative \partial m x0 f[x(t)] f[x(t)]
y f(y) f(y)
x0 −¯x1 \sqrt ∆t
¯x1 −x0 \sqrt ∆t . (17.823) Then the derivative of the functional (17.815) becomes \partial m
** \Phi[V (x)] Hm ¯x1 −x0 \sqrt ∆t
x(t) . (differentiated ensemble average) (17.824)
derivatives may be computed merely by reweighting via x1 and xN−1. Of course, each derivative supplies an additional factor ∆t−1/2 ∝N 1/2, making the fluctuations larger for each successively higher derivative. Note that the expression (17.824) also applies to ordinary Wiener processes if ∆t is replaced by 2∆t, and ¯x1 is replaced by x1. Although in principle the expression (17.824) applies to arbitrarily high derivatives, it may only prac- tical for relatively small m, depending on the nature of the solution to the path average. The factor of (∆t)−m/2 is something we can put aside for the moment, since this problem can at least be tamed. Recall that the harmonic-oscillator eigenfunctions have the (normalized) form
p
Hm(x) e−x2/2, (17.825) which grow in width as \sqrtm. Roughly, this means that the weighting factor in Eq. (17.824) favors increasingly large mean steps |¯x1| with increasing m. This width can match the original probability density quite poorly for large m, in which case it could be advantageous to absorb the Hermite polynomial into the probability measure for the paths—this avoids rare events in the Gaussian tails from making large contributions. How- ever, doing so introduces a normalization factor that grows exponentially with m. It is only in cases where the ensemble-average derivative grows similarly that the relative accuracy in a Monte-Carlo calculation does not suffer (an average that depends on x−s for s > 0 is one case where this is okay).
17.15 Differentiation of Stochastic Quantities 17.15.1.3 Variance Reduction Now to deal with the factor of (∆t)−m/2 that causes Eq. (17.824) to become numerically inefficient for larger derivatives. Again, this is because, generally speaking, simulations of path integrals and path averages are accurate in the limit N −\rightarrow \infty, but the fluctuations of the path average grow as N m/2. This obviously complicates accurate simulation of the path averages difficult. The key idea is that under certain conditions, the ‘‘payoff’’ functional \Phi[V (x)] is approximately inde- pendent of the path samples near the beginning (x1, x2, . . .) and end (xN−1, xN−2, . . .) of the path. This
the boundary. These parts of the path thus contributes very little to the mean sojourn time. Under these conditions, we can successively integrate over the early coordinates (x1, x2, . . .) and the late coordinates (xN−1, xN−2, . . .) to perform a ‘‘partial average’’ over paths. More specifically, we can hold x0 and x2 fixed, while averaging over all possible values of x1. Then while holding x3 fixed, average over all possible values of x2, and so on until the error associated with the averaging is no longer negligible. To arrive at the appropriate, partially averaged version of Eq. (17.824), consider the integral I := Z dxj Hm
x0 −¯xjk p
! e−(x0−xN−k)2/2k∆t \sqrt 2\pik∆t e−(xj−x0)2/2j∆t
\sqrt 2\pi∆t , (17.826) which consists of the Gaussian and Hermite-polynomial weights for a step of size j∆t at the beginning of
for the Gaussian step of size ∆t from xj to xj+1; and the partial-averaging integral over xj. The mean displacement of the first and last steps here is ¯xjk := kxj + jxN−k j + k . (17.827) Then completing the squares on the first product of Gaussian factors and then the second product, we can rewrite this as I = Z dxj Hm
x0 −¯xjk p
! e−(x0−¯xjk)2(j+k)/2jk∆t p
e−(xN−k−xj)2/2(j+k)∆t p
\sqrt 2\pi∆t = Z dxj Hm
x0 −¯xjk p
! e−(x0−¯xjk)2(j+k)/2jk∆t p
p
p
p
Z d¯xjk Hm
x0 −¯xjk p
! e−(x0−¯xjk)2(j+k)/2jk∆t
p
Z d¯xjk Hm
x0 −¯xjk p
! e−(x0−¯xjk)2(j+k)/2jk∆t
(17.828) where we have defined the offset
, (17.829)
Chapter 17. Stochastic Processes and we used
k
k ¯xjk − (j + 1)
kxj+1
k ¯xjk −¯x(j+1)k , (17.830) where now
. (17.831) Then using the convolution formula78
Hn x \beta
=
x p
!
(17.832) with
s 2k2∆t
s 2jk∆t (j + k) (17.833) such that
\beta p
s
(17.834) we find
p
m/2 \times Hm
x0 −¯x(j+1)k p
!
(17.835) Rearranging to make the recursion more clear, we have j + k jk m/2 I =
(j + 1)k m/2 Hm
x0 −¯x(j+1)k p
!
p
p
=
(j + 1)k m/2 Hm
x0 −¯x(j+1)k p
! \times e−(x0−xN−k)2/2k∆t \sqrt 2\pik∆t
p
(17.836) where we are basically ‘‘unravelling’’ the first square that we completed in Eq. (17.828). The same argument may be made for integrating over xN−k by the time-symmetry of the Brownian bridge. Thus the expression
\partial m
j + k 2jk∆t
\Phi[V (x)] Hm
¯xjk −x0 p
x(t) = t1 + t −t2 2t1(t −t2)
\Phi[V (x)] Hm
¯xjk −x0 p
x(t) , (differentiated ensemble average, with partial average) (17.837) 78Daniel A. Steck, Classical and Modern Optics (2006), available online at http://steck.us/teaching.
17.15 Differentiation of Stochastic Quantities as we have just proven by induction, where again ¯xjk := kxj + jxN−k j + k
t1 + t −t2 (17.838) in both discrete and continuous notation, where t1 = j∆t and t−t2 = k∆t, with a total running time t of the
for all 0 < i < j and N −k < j < N. In the simpler, symmetric case where j = k, we have \partial m
** \Phi[V (x)] Hm ¯xjj −x0 \sqrtj∆t
x(t)
** \Phi[V (x)] Hm
¯xjj −x0 p t1/t
x(t) , (differentiated ensemble average, with symmetric partial average) (17.839) where ¯xjj = xj + xN−j
. (17.840) The advantage here is clear, especially in the last expression. In Eq. (17.824), the statistical fluctuations in the path average grew as N m/2. But since typically j and k are chosen for some fixed times j and k (e.g., as in the sojourn-time example), the fluctuations here are independent of N. Of course, since t1 < t/2, the fluctuations still increase by a factor of (t/t1)1/2 for each derivative, but this is a huge improvement when working with large N. One remaining detail is to examine the statistics of ¯xjk in the Hermite polynomial, to make sure it does not cause any problem. For a discrete Brownian bridge Bn of unit running time, the covariance is given by [Problem 17.10]
Bn Bm
N 2 , (17.841) which gives
N (17.842) and
N (17.843) for the symmetric case. Then setting x0 = 0, we thus have Var[¯xjj] = 1 h Var(xj) + Var(xN−j) + 2 Cov(xj, xN−j) i = ∆t 2N h j(N −j) + j2i = j∆t 2 . (17.844) Thus the typical ¯xjj is of the order \sqrtj∆t, which is precisely what we see in the denominator of the Hermite polynomial. Again, to use this partially averaged path integral as a numerical technique, it can be helpful to incorporate the Hermite polynomial into the sampling distribution for the path itself. To simplify the math somewhat, we will stick to the symmetric expression (17.839) for the path integral. We know from our discussion of the Brownian bridge above that we may choose the path coordinates at times j∆t = t1 and
Chapter 17. Stochastic Processes they are correlated because they are part of a bridge, with covariance j2∆t2 = t 2 1 . We can decouple this correlation by working with ¯xjj and \deltaxjj := xj −xN−j
, (17.845) such that xj = ¯xjj + \deltaxjj and xN−j = ¯xjj −\deltaxjj. We already computed Var[¯xjj], and we also need
h Var(xj) + Var(xN−j) −2 Cov(xj, xN−j) i = ∆t 2N h j(N −j) −j2i
2N . (17.846) Of course, the simplification comes because of the vanishing covariance:
h Var(xj) −Var(xN−j) i = 0. (17.847) Thus the ensemble average in the path integral (17.839) implies a probability measure
" e−(¯xjj−x0)2/j∆t
p
¶
\sqrt 2\pi∆t
\sqrt 2\pi∆t , (17.848) where the bracketed factor gives the partially averaged first and last steps, and the rest of the factors give the small, in-between steps. Note that there is a factor of 1/2 that is now omitted from the bracketed factor, owing to the choice of ¯xjj and \deltaxjj as path variables (before transforming to these variables the corresponding normalization factors had a factor of 2 associated with each \pi); this factor of 1/2 is absorbed into the integration measure via the Jacobian derivative in the variable transformation. Then the Hermite polynomial may be grouped with the first Gaussian factor
( \etam
Hm ¯xjj −x0 \sqrtj∆t e−(¯xjj−x0)2/j∆t
p
)
\sqrt 2\pi∆t
\sqrt 2\pi∆t , (17.849) where the factor \eta−1 m :=
Z \infty −\infty dx Hm x \sqrtj∆t
\sqrt\pi Z \infty dx |Hm(x)| e−x2 (17.850) normalizes the Hermite–Gaussian probability density for ¯xjj. Then to summarize, we should take
p j∆t ¯z,
r j(1 −2j/N)∆t z, (17.851) where z is standard normal, and ¯z is chosen from the density
(17.852) Then with the path measure (17.849), the path average (17.839) becomes \partial m
m (j∆t)−m/2 ** \Phi[V (x)] sgn Hm(¯z) ++ Hm , (differentiated ensemble average, with Hermite–Gauss path) (17.853)
17.15 Differentiation of Stochastic Quantities where note that we kept the sign of the Hermite polynomial in the path-average argument rather than the measure, and we have included the factor of 1/2 that we discussed after Eq. (17.848). Then with the new path measure, xj and xN−j can be calculated via xj = ¯xjj + \deltaxjj and xN−j = ¯xjj −\deltaxjj, while the rest of the path can be constructed as a Brownian bridge connecting xj to xN−j in time t −2t1. 17.15.2 Stochastic Differentiation For a stochastic process y(t) satisfying the SDE
(17.854) the notion of a stochastic derivative arises by considering the effect on y(t) of a perturbation to dW(t0)
for our purposes can be written as the partial derivative79
\partial y(t) \partial [dW(t0)]. (17.855) (stochastic derivative) Another common, related notation is Dy, which is essentially the collection of all stochastic derivatives
This derivative generalizes the notion of a functional derivative to stochastic processes. Thus, Dy is essentially the functional derivative \deltay/\delta(dW), regarding \delta(dW)(t) as a perturbation to dW(t). And while Dy acts as a gradient, Dt0y acts as a directional derivative in the ‘‘direction’’ \delta(t −t0), in the same sense as the functional derivative via a perturbation only at time t0: \deltay/\delta(dW)(t0). 17.15.2.1 Example: Stochastic Derivative of the Wiener Process As a first example of stochastic differentiation, we will begin by computing Dt0W(t). Since
Z t dW(t′), (17.856) the function to differentiate is a simple sum of increments dW(t′) and we are differentiating with respect to one of them. Then we can use \partial [dW(t)]
(17.857)
and zero otherwise. Thus,
Z t
Z t dt′ \delta(t0 −t′), (17.858) or simply
(17.859) provided, of course, 0 \le t0 \le t. Otherwise this derivative vanishes. 17.15.2.2 Example: Stochastic Derivative of Additive Diffusion As a slightly more sophisticated example, let’s now compute Dy, where
(17.860) 79See Arturo Kohatsu–Higa and Miquel Montero, ‘‘Malliavin Calculus in Finance,’’ in Handbook of Computational and Numerical Methods in Finance Svetlozar T. Rachev, Ed. (Springer, 2004), pp. 111-74, especially p. 130. This derivative is more generally known as the Malliavin derivative. For a more rigorous and sophisticated development, see David Nualart, The Malliavin Calculus and Related Topics, 2nd ed. (Springer, 2006), Section 1.2, p. 24.
Chapter 17. Stochastic Processes We can do this by generalizing the solution from the previous section, which leads straightforwardly to
(17.861) Thus, a simple relabeling gives the result80 Dy = \beta. 17.15.2.3 Variation Process Going back to the stochastic derivative Dt0y(t), where y(t) satisfies the prototype SDE (17.854), we can think of the derivative itself as a stochastic process (considering t to be the time variable, with fixed t0). We can refer to this as the derivative process or the variation process. This describes the modified evolution of y(t) as a linearized displacement, given a perturbation to dW(t0). That is, the perturbed evolution is
perturbation to dW(t0). Taking the stochastic derivative of Eq. (17.854) by applying the Dt0 operator, we find
(17.862) Using Eq. (17.857) again in the last term, and dividing through by Dt0y(t), we arrive at the SDE d[Dt0y(t)] Dt0y(t)
Dt0y(t) \delta(t −t0) dt, (SDE for variation process) (17.863) which yields the evolution of the variation processs.
differential equation d[Dt0y(t)] Dt0y(t)
(17.864)
Rewriting the left-hand side of the homogeneous SDE, d log[Dt0y(t)] + d[Dt0y(t)] 2[Dt0y(t)]2
(17.865)
dt + \beta′[y(t)] dW(t). (17.866) Integrating from t0 to t,
Z t t0
dt′ + Z t t0 \beta′[y(t′)] dW(t′). (17.867) Exponentiating this equation leads to the expression
Z t t0
dt′ + Z t t0 \beta′[y(t′)] dW(t′) , (solution for variation process) (17.868) which solves the SDE (17.863).
where it is not obvious that W(h) can be interpreted as the integral of h dW, but see Example 1.1 on p. 3 in the corresponding lecture notes at http://www.math.wisc.edu/~kurtz/NualartLectureNotes.pdf; for a more direct comparison, see Giulia Di Nunno, Bernt Øksendal, and Frank Proske, Malliavin Calculus for Lévy Processes with Applications to Finance (Springer, 2009), p. 29, Eq. (3.8).
17.15 Differentiation of Stochastic Quantities 17.15.3 Integration by Parts The notion of a stochastic derivative also gives a kind of integration by parts, which can be helpful in transforming certain path averages. As an example, consider the model path average
DD \delta[W(t) −c] EE W (t). (17.869) (model ensemble average) That is, we are averaging over Wiener paths W(t), which only ‘‘score’’ when the endpoint W(t) matches the desired ending point c. Since we know the explicit probability density for W(t), we can compute this explicitly as
Z \infty −\infty dx e−x2/2t \sqrt
\sqrt 2\pit . (17.870) However, as we will see, the general idea here will apply to more general stochastic processes, where we may not have an analytic expression for the density of paths. From the standpoint of numerical simulation, Eq. (17.869) is very difficult, because only a subset of paths of zero measure contribute at all to the average, each with an infinite ‘‘score.’’ In a numerical simulation, the delts function could, for example, be changed into a function of finite height and width, but the variance of a sample path average would still be very large for any function narrow enough to give an
Eq. (17.592), but this presumes that we know the distribution of W(t). In this case we do, but we would like to explore more general methods for cases where the distribution is not known. The approach we will try here is to rewrite the path functional as
2\partial c DD sgn[W(t) −c] EE W (t), (17.871) which removes the delta function at the expense of introducing a derivative. Now switching to shifted Brownian bridges
(17.872) with source point x0, we can write the path functional as
2\partial c DD sgn[x(t) −c] EE x(t)
x0=0 = 1 2\partial x0 DD sgn[x(t) −c] EE x(t)
x0=0 , (17.873) Then using the likelihood-ratio estimator (17.812),
DD sgn[x(t) −c] \partial x0 log f[x(t)] EE x(t)
x0=0 , (17.874) where f[x(t)] is the probability density of the paths x(t), which we regard as containing all the dependence on x0. For Wiener paths,
2t −log \sqrt 2\pit, (17.875) and thus
t
t . (17.876) Then Eq. (17.874) becomes
2t DD sgn[W(t) −c] W(t) EE W (t). (17.877) (alternative path average) We can verify by direct integration with respect to the probability desnity that this gives the correct answer. However, unlike the original path functional (17.869), every path W(t) contributes value to the average. The ‘‘cost’’ is a weight of W(t)/2t in addition to the integral of the delta function, but everything is well-behaved
Chapter 17. Stochastic Processes here from a numerical standpoint. Note that other antiderivatives of the delta function are possible: the Heaviside function is an obvious choice, but then half of the paths contribution no ‘‘information’’ to the path average. However, in deriving Eq. (17.877), we again used the explicit probability density for W(t), as an illustration of the method. To perform a more general calculation, we cannot do this. Consider now the model path integral with more general stochastic processes:
DD \delta[y(t) −c] EE y(t)
(17.878) (model ensemble average) To replace the delta function here, first we will compute the chain-rule derivative of the step function:
(17.879) Solving for the delta function and integrating over time gives
2t Z t dt0 Dt0sgn[y(t) −c] Dt0y(t) . (17.880) Now if we include an explicit probability measure for dW(t0), we can integrate by parts with respect to dW(t0) to write Z dW(t0) e−dW 2(t0)/2dt \sqrt 2\pidt
Z dW(t0) e−dW 2(t0)/2dt \sqrt 2\pidt dW(t0) dt sgn[y(t) −c]. (17.881) Thus, under an ensemble average integration by parts amounts to writing
dt sgn[y(t) −c]. (17.882) If we then introduce an expectation value in Eq. (17.880) and perform this integration by parts, we find DD \delta[y(t) −c] EE = 1 2t Z t dt0 ** Dt0sgn[y(t) −c] Dt0y(t) ++ = 1 2t Z t dt0 ** dW(t0) dt sgn[y(t) −c] Dt0y(t) ++ , (17.883) noting that Dt0y(t) is independent of dW(t0). Simplifying, we have81
2t ** sgn[y(t) −c] Z t dW(t0) Dt0y(t) ++ , (17.884) (alternate path functional)
that the derivative (17.868) in the denominator of the integral reduces to the simpler form
Z t t0 dt′ \alpha′[y(t′)] , (17.885) so that Eq. (17.884) becomes
2t ** sgn[y(t) −c] Z t dW(t0) exp − Z t t0 dt′ \alpha′[y(t′)]
, (17.886) such that a drifting path still introduces a nontrivial weighting function that measures the drift-induced focusing or divergence of paths. 81A more general version of this argument is given in Peter K. Friz, ‘‘Malliavin Calculus in Finance,’’ Section 7.1, available at http://www.math.nyu.edu/phd_students/frizpete/finance/my_case_lecture/malliavin_lecture.pdf. See also Eric Fournié, Jean-Michel Lasry, Jérôme Lebuchoux, and Pierre-Louis Lions, ‘‘Applications of Malliavin calculus to Monte-Carlo methods in finance. II,’’ Finance and Stochastics 5, 201 (2001), Section 4.1 (doi: 10.1007/PL00013529).
17.15 Differentiation of Stochastic Quantities 17.15.3.1 Digression: Faddeev–Popov Approach to Conditional Averages An alternative approach, based on the method of Faddeev and Popov, can work to give regularized expressions for the conditional path averages in the previous section. To start with the simpler case of Eq. (17.869),
DD \delta[W(t) −c] EE W (t), (17.887) (model ensemble average) as a first step, let’s write out the probability measure explicitly in discrete form with N path steps:
N 2\pit N/2 Z
−N 2t N−1 X j=0 ∆Wj 2 \delta W(t) −c . (17.888) Note that we can regard W(t) here as an explicit sum:
N−1 X j=0 ∆Wj. (17.889) Now the key is to note that because we are integrating over all possible values of ∆Wj, the integral A(c) is invariant if we shift the variable ∆Wj. In particular, suppose that we introduce a scalar quantity \omega, and shift each coordinate according to ∆Wj −\rightarrow ∆Wj −\omega N . (17.890) Then A(c) is independent of \omega, and we can regard this shift as a kind of gauge transformation with gauge parameter \omega. Implementing this gauge transformation, Eq. (17.888) becomes
N 2\pit N/2 Z
−N 2t N−1 X j=0 ∆Wj −\omega N 2 \delta W(t) −c −\omega . (17.891) Now suppose that we introduce some function g(\omega), with the only requirement that it be normalized such that Z
(17.892) Then we don’t change the value of Eq. (17.891) by multiplying by g(\omega) and integrating, thus taking a linear combination of the same value. The delta function disappears in the result:
N 2\pit N/2 Z
−N 2t N−1 X j=0 ∆Wj −W(t) −c N 2 g W(t) −c . (17.893) In the exponentiated summand, we can multiply out the quadratic factor to obtain ∆Wj −W(t) −c N 2 = ∆Wj 2 −2∆Wj W(t) −c N + W(t) −c 2 N 2 , (17.894) and thus when summed over j, we find N−1 X j=0 ∆Wj −W(t) −c N 2 = N−1 X j=0 ∆Wj 2 + c2 N − W(t) 2 N . (17.895) In this case the path integral (17.893) becomes
** exp " W(t) 2 −c2 2t # g W(t) −c ++ W (t) (17.896)
Chapter 17. Stochastic Processes after again hiding the path measure in the ensemble average. At this point we are free to choose g(\omega) in order to simplify or improve the statistical behavior of this ensemble average. For example, if we choose g(\omega) to be a Gaussian of variance t,
\sqrt
(17.897) then the path average (17.896) becomes
\sqrt 2\pit DD ecW (t)/tEE W (t), (17.898) (alternate path integral) which attains the correct value of e−c2/2t/ \sqrt 2\pit. On the other hand, if we take g(\omega) to be centered at −c instead of at 0,
\sqrt
(17.899) the path integral (17.896) takes on the average value directly, with zero variance. This illustrates the advantage of adapting g(\omega) to the problem as well as can possibly be done. Of course, this approach becomes more complicated for the more general diffusion problem (17.878):
DD \delta[y(t) −c] EE y(t)
(17.900) (model ensemble average) The general approach, however, will be the same. Writing out the details of the path integral, we have
N 2\pit N/2 Z
−N 2t N−1 X j=0 ∆Wj 2 \delta y(t) −c (17.901) when integrating in terms of the unshifted path variables, where
Z t dt′ \alpha y(t′) + Z t dW(t′) \beta y(t′) . (17.902) Implementing the same gauge transformation (17.890), the result is
N 2\pit N/2 Z
−N 2t N−1 X j=0 ∆Wj −\omega N 2 \delta y(t) −c −\omega N N−1 X j=0 \partial y(t) \partial ∆Wj , (17.903) where in the delta function we have expanded to O(N −1). Then we can define the shorthand for the last term in the argument of the delta function, A := 1 N N−1 X j=0 \partial y(t) \partial ∆Wj = 1 t Z t dt0 Dt0y(t), (17.904) in both discrete and continuous notation. This quantity is written more explicitly using the expression
via Eq. (17.885) to give A = 1 t Z t dt0 exp Z t t0 dt′ \alpha′[y(t′)] . (17.905) Then we can proceed by introducing the function g(\omega) and integrating over all \omega, with the result
N 2\pit N/2 Z
−N 2t N−1 X j=0 ∆Wj −y(t) −c N|A| 2 1 |A| g y(t) −c |A| , (17.906)
17.15 Differentiation of Stochastic Quantities¶
again gives
** |A| exp " y(t) −c 2t|A|
2W(t) − y(t) −c |A| !# g y(t) −c |A|
y(t) . (17.907) Again, we can choose the form of g(\omega) to simplify the problem, but for a general drift function \alpha[y], we cannot choose it to arrive directly at the answer. Choosing the gaussian form (17.897) again gives the result
\sqrt 2\pit ** |A| exp " y(t) −c t|A|
W(t) − y(t) −c |A|
y(t) , (17.908) which does not make the path integral particularly simple, but at least it casts the path contribution in
make the probably better choice (17.899), which leads to
\sqrt 2\pit ** |A| exp " y(t) −c t|A|
W(t) −c − y(t) −c |A|
y(t) , (alternate path integral) (17.909)
Note in the path integrals how y(t) is accompanied by a factor of A; intuitively, from the discussion in Section 17.15.2.3, this factor gives a time-averaged measure of the divergence of trajectories in the vicinity of y(t), and thus accounts for any focusing or defocusing effects due to gradients in \alpha(y).
Chapter 17. Stochastic Processes 17.16 Exercises Problem 17.1 Geometric Brownian motion describes the time-dependence price S of a stock according to the SDE dS = S
(17.910) where µ represents the steady growth of the stock value, and \sigma represents the stock volatility (which is assumed to be constant within this model). Show that
(17.911) satisfies the above SDE. Problem 17.2 Show that the vector SDE (17.145)
(17.912) is equivalent to the Fokker–Planck equation (17.146),
2\partial i\partial jDij(x, t)f(x, t), (17.913) where
(17.914) Problem 17.3 By formally resumming the Taylor series expansion, compute exp(dN). Problem 17.4 Recall the Poisson distribution in terms of the single parameter \lambda has the form
n! , (17.915) where n is a nonnegative integer. Also, suppose that N is a Poisson-distributed random variable. (a) Show that P(n) is normalized.
Problem 17.5 Consider a monochromatic field with added, white (Gaussian) frequency noise,
(17.916) so that the instantaneous frequency is d\phitotal dt
dt , (17.917)
17.16 Exercises¶
process by
(17.918) Use the Wiener–Khinchin theorem to show that the spectrum is Lorentzian with full width at half maximum \gamma. Problem 17.6
I := Z dy1 . . . dyN−1 exp −N N X j=1 (yj −yj−1)2 . (17.919) (a) Evaluate this integral by using the recurrence (17.305) yn = y′ n \sqrtcn
cn . (17.920) to decouple the integral into a product of Gaussian integrals.
y(t) and W(1), to evaluate this integral, to show the consistency of these approaches. Problem 17.7 Show that the integral expression (17.335)
Z t dW(t′) 1 −t′ (17.921) solves the SDE (17.334) dB = − B 1 −t dt + dW. (17.922) Problem 17.8 Derive the correlation function
B(t) B(t′)
(17.923) (a) using the definition (17.287)
(17.924) (b) using the definition (17.335)
Z t dW(t′) 1 −t′ . (17.925) Problem 17.9 Show that if B(t) is a standard Brownian bridge, the covariance of two Brownian-bridge increments
Chapter 17. Stochastic Processes Problem 17.10 (a) Using the recurrence for the finite Brownian bridge [Eq. (17.310)], B0 = 0 Bn = zn s N −n
N −n N −n + 1 Bn−1, n = 1, . . . , N −1 BN = 0, (17.926) derive an explicit formula for the Brownian-bridge samples Bn in terms of (only) the standard-normal deviates zn. (b) Show that
Bn Bm
N 2 , (17.927) and show that this is consistent with what you know for the continuous-time limit. (c) Show that
∆Bn ∆Bm
(N −n −1) N(N −n) + min(n, m) N 2[N −min(n, m)] + (\deltanm −1) \sqrt N(N −n)(N −m) . (17.928) and show that this is consistent with what you know for the continuous-time limit, if we take
(17.929) Problem 17.11 For the standard Brownian bridge B(t): (a) Use the statistics derived in Section 17.7 derived for the Brownian bridge to write down the
Problem 17.12 Given a state-dependent diffusion function of the form
(17.930) where the diffusion rate is \sigma0 for y < d and \sigma1 for y > d, show that an explicit solution (17.349) for the SDE (17.338)
(17.931)
1 −\sigma1 \sigma0 \Theta[B(t) −d], (17.932) where \Theta(x) is the Heaviside function and B(t) is a standard Brownian bridge. Problem 17.13 (a) Show that
X x0\in f −1(0) \delta(x −x0) |f ′(x0)| , (17.933)
17.16 Exercises where the sum is over the (isolated) zeros of f(x). Good ways to proceed are: to write the delta function as a derivative of the step function \Theta[f(x)], or to treat the delta function as a limit of a more ‘‘reasonable’’ function, such as a box function of width a and height 1/a. (b) Show that
X x0\in f −1(0) \delta′(x −x0) f ′(x0)|f ′(x0)| + f ′′(x0) \delta(x −x0) |f ′(x0)|3 . (17.934) The same suggested approaches apply here. Problem 17.14
(17.935) Using only this expression (along with basic probability theory and your own cunning) to obtain the analogous expression [Eq. (17.372)]
d \sqrt 2t (17.936)
Hint: rescale the bridge probability to account for the alternate time interval, and interpret this as a
Problem 17.15
some time 0 \le t′ \le T, and show that it approaches unity as T −\rightarrow \infty. Problem 17.16 Rederive the probability density 17.710
h 1 −e[c2−(|d|+|c−d|)2]/2ti
t
(17.937) for the local time of a Brownian bridge Bt(0\rightarrow c)(t) pinned to c at time t. Use an alternate derivation:
with respect to c to obtain the pinning delta function.
single extra integral over k. Also, note that by avoiding the c-derivative and completing the calculation, it is possible to obtain the joint probability of W(t) occupying [0, \infty) and the local time ℓ[W(t); d] taking on a particular value. Problem 17.17 (a) Using the calculation of the probability density for the local time ℓ[W(t); d] as a template, derive
Chapter 17. Stochastic Processes the following formula for the expectation value82 \lambda Z \infty dt e−\lambdat
exp −sℓ[W(t); d] −s′ℓ[W(t); d′] W (t) = 1 − s h\sqrt
1 −e− \sqrt 2\lambda|d−d′|i e− \sqrt
1 −e− \sqrt 2\lambda|d−d′|i e− \sqrt 2\lambda|d′| \sqrt
\sqrt
−ss′e−2 \sqrt 2\lambda|d−d′| , (17.938) giving the dual moment-generating function for the two local times ℓ[W(t); d] and ℓ[W(t); d′], given that the process W(t) is stopped exponentially at rate \lambda. (b) Repeat for a Brownian bridge Bt(0\rightarrow c)(t), pinned at the exponentially distributed stopping time t
\lambda Z \infty dt e−\lambdat e−c2/2t \sqrt 2\pit
exp −sℓ[Bt(0\rightarrow c)(t′); d] −s′ℓ[Bt(0\rightarrow c)(t′); d′] Bt(0\rightarrow c)(t′) = r \lambda 2 e− \sqrt 2\lambda|c| 1 − s h\sqrt
e− \sqrt 2\lambda|d| −s′e− \sqrt
e− \sqrt 2\lambda(|c−d|−|c|) \sqrt
\sqrt
−ss′e−2 \sqrt 2\lambda|d−d′| − s′ h\sqrt
e− \sqrt 2\lambda|d′| −se− \sqrt
e− \sqrt 2\lambda(|c−d′|−|c|) \sqrt
\sqrt
−ss′e−2 \sqrt 2\lambda|d−d′| , (17.939) where the Gaussian factor on the left-hand side indicates this expectation value was taken with respect to a joint distribution for the local times and for the stopping point of B(t). Problem 17.18 (a) Using the representation
a\rightarrow 0+
a (17.940) and the results of Problem 17.17, derive the following formula for the expectation value \lambda Z \infty dt e−\lambdat
exp −sℓ′[W(t); d] W (t) = 1 −e− \sqrt 2\lambda|d|, (17.941) giving the moment-generating function for the local-time derivative ℓ′[W(t); d] and ℓ[W(t); d′], given that the process W(t) is stopped exponentially at rate \lambda. (b) Invert the \lambda and s Laplace transforms in this result to obtain an expression for the probability density fℓ′(x) of ℓ′[W(t); d]. You should obtain something a bit funny; what gives? Problem 17.19
(a) Using the probability method of Eqs. (17.592), show that
\delta[BT (\tau) −c] F[BT (t)]
F BT (t)
, (17.942) 82Andrei N. Borodin and Paavo Salminen, Handbook of Brownian Motion—Facts and Formulae, 2nd ed. (Birkhäuser, 2002), p. 177, Eq. (1.18.1). 83This is similar to, but not quite the same as, Borodin and Salminen, op. cit., p. 177, Eq. (1.18.5).
17.16 Exercises where the conditional ensemble average on the right-hand side refers to Brownian bridges pinned to c at time \tau (0 < \tau < T), and
\sqrt
(17.943) is the centered Gaussian distribution with variance \sigma2. (b) Prove the same result by adapting the explicit-integration method of Section 17.12.4. Problem 17.20
DD \delta[y(t) −c] EE y(t)
(17.944)
the path integral converges to the same value.