Derivation of Wiener process first passage times using probability generating function?

Question

I would like to find the distribution of first passage times in a simple Wiener process using the idea of probability generating functions. Thus there will be, at certain point, a limiting step to go from the discrete method of pgf's to a distribution function. I have seen the derivation of this distribution using the reflection principle, but I am keen to do it by the pgf method, although I am not sure if it can work. I am doing the following:

Let $\tau(m)$ be the first time a level $m>0$ is reached after n steps and $S_0=0$: $$ \tau(m)=min\{n\geq0:S_n=m\}\ $$ The pgf of $\tau(m =1)$ is $$ f_{\tau}(s)=\sum_{n=1}^{\infty}s^nP\{\tau=n\} $$ Similarly for m>1. Then noticing that the pgf of a sum of random variables is the product of the pgf's I get the following $$ f_{\tau}(s)=\frac{(1-\sqrt{1-s^2})}{s} $$ Now I suspect I need to do some sort of limiting step to get a distribution for $\tau$, but I do not know how proceed from here. My aim is to recover the well known distribution for first passage times $f(\tau)=\frac{a}{\sqrt{2\pi}}t^{-3/2}exp(-a^2/2t)$
Any advice will be highly appreciated

Answer 1

You can indeed get the distribution of first passage times of the Wiener process in this way, but it's a bit more effort than the standard way. All references in this answer are to Analytic Combinatorics, which you can download for free.

As Robert showed, the probability generating function for the first passage time $\tau(n)$ to $n$ is

$$g_n(s) = \left(\frac{1-\sqrt{1-s^2}}{s}\right)^n\;.$$

I'm using $n$ instead of $m$ to match the variables in Theorem VIII.8 (p. 587). First we need to bring this into a form in which we can apply the theorem, that is, we need to satisfy conditions $\mathbf L_1$ and $\mathbf L_2$ in Section VIII.8.1 (p. 586):

$$g_n(s)=s^n\left(\frac{1-\sqrt{1-s^2}}{s^2}\right)^n=:s^nB(s^2)^n$$

with

$$B(z)=\frac{1-\sqrt{1-z}}z\;.$$

To find the distribution of the first passage time $t$ at $a$, we need to take the step sizes $\Delta x$ and $\Delta t$ to zero while keeping $\Delta x^2/\Delta t=:\beta$ constant. Then $n=a/\Delta x$ and $N=t/(2\Delta t)=\beta t/(2\Delta x^2)$, where the factor $2$ arises from $z=s^2$. We need to find the asymptotic behaviour of the coefficients of $B(z)$ as $\Delta x$ goes to zero and thus both $n$ and $N$ go to infinity.

First note that the spread $T$ defined in Section VIII.8.1 is infinite in this case, so we don't have to worry about restrictions on the ratio $\lambda=N/n=\beta t/(2a\Delta x)$, which goes to infinity as $\Delta x\to0$. We have to find the root $\zeta$ of

$$\zeta\frac{B'(\zeta)}{B(\zeta)}=\lambda\;,$$

which turns out to have the nice form

$$\zeta=1-\frac1{(2\lambda+1)^2}$$

with

$$B(\zeta)=1-\frac1{2(\lambda+1)}\;.$$

We also need

$$\xi=\frac{\mathrm d^2}{\mathrm d\zeta^2}\left(\log B(\zeta)-\lambda\log\zeta\right)=\frac{B''(\zeta)}{B(\zeta)}+O\left(\lambda^2\right)=(2\lambda)^3+O\left(\lambda^2\right)\;.$$

Now we can apply Theorem VIII.8 to obtain the distribution function as the limit

$$ \begin{align} \lim_{\Delta x\to0}\frac{[s^{2N}]g_n(s)}{\Delta t} &= \lim_{\Delta x\to0}\frac{[z^{N-n/2}]B(z)^n}{\Delta t} \\ &= \lim_{\Delta x\to0}\frac{[z^N]B(z)^n}{\Delta t} \\ &= \lim_{\Delta x\to0}\frac{B(\zeta)^n}{\zeta^{N+1}\sqrt{2\pi n\xi}\Delta t} \\ &= \lim_{\Delta x\to0}\frac{\left(1-\dfrac1{2(\lambda+1)}\right)^n}{\left(1-\dfrac1{(2\lambda+1)^2}\right)^{N+1}\sqrt{2\pi (2\lambda)^3\dfrac a{\Delta x}}\,\Delta t} \\ &= \lim_{\Delta x\to0}\frac{\left(1-\dfrac{a\Delta x}{\beta t}\right)^{a/\Delta x}}{\left(1-\left(\dfrac{a\Delta x}{\beta t}\right)^2\right)^{\beta t/2(\Delta x^2)+1}\sqrt{2\pi\left(\dfrac{\beta t}{a\Delta x}\right)^3\dfrac a {\Delta x}}\,\Delta x^2/\beta} \\ &= \frac{\exp\left(-\dfrac{a^2}{\beta t}\right)}{\exp\left(-\dfrac{a^2}{2\beta t}\right)\sqrt{\dfrac{2\pi\beta t^3}{a^2}}} \\ &= \frac{a\exp\left(-\dfrac{a^2}{2\beta t}\right)}{\sqrt{2\pi\beta t^3}}\;. \end{align} $$

Then setting $\beta=1$ gives the desired result; it appears you forgot a square root around $2\pi$.

Answer 2

Let $g_m(s) = E[s^{\tau(m)}]$ be the probability generating function for the first passage time $\tau(m)$ to $m$. Since $\tau(0) = 0$ with probability $1$, $g_0(s) = 1$. Otherwise, there must be a first step which with equal probabilities takes you to $\pm 1$, so $g_m(s) = s (g_{m+1}(s) + g_{m-1}(s))/2$. The general solution of this homogeneous linear difference equation is $g_m(s) = A {r_+}^m + B {r_-}^m$ where $r_\pm = \frac{1 \pm \sqrt{1-s^2}}{s}$ and $A$, $B$ are constants. However, since we need $\lim_{s \to 0+} g_m(s)$ to exist, $r_+$ is excluded. Since $g_0(s) = 1$, we conclude that $$g_m(s) = {r_-}^m = \left(\frac{1-\sqrt{1-s^2}}{s}\right)^m$$

Answer 3

This is to expand on some comments on @joriki's answer, which might complete it. The story began when I (somewhat foolishly) suggested:

In the case at hand, a shortcut might be to show that $g_{a\sqrt{n}}(\mathrm e^{-\lambda/n})\to\mathrm e^{-a\sqrt{2\lambda}}=\mathrm E(\mathrm e^{-\lambda\tau_a})$, where $\tau_a$ denotes the first hitting time of the level $a$ by a standard Brownian motion, and to identify the RHS as the Laplace transform of the correct density.

Prompted (quite rightfully) for some less cursory explanations, I then added:

Using $s=\mathrm e^{-\lambda/n}$ and $k=a\sqrt{n}$ in the expression of $g_k(s)$ stated in your post yields $1-s^2\sim2\lambda/n$, $(1-\sqrt{1-s^2})^{a\sqrt{n}}\to\mathrm e^{-a\sqrt{2\lambda}}$ and $s^{a\sqrt{n}}\to1$, hence $g_k(s)\to\mathrm e^{-a\sqrt{2\lambda}}$. Now, using $x=a\sqrt{2\lambda}$ and introducing $H(x)=\int\limits_0^{+\infty}\mathrm e^{-x^2s/2-1/(2s)}\mathrm ds/\sqrt{2\pi s^3}$, one sees that $H(0)=1$ (use $s=1/t^2$) and $H'(x)=-H(x)$ (use $x^2s=1/u$) hence $H(x)=\mathrm e^{-x}$, QED.

To which @joriki (pitilessly but legitimately) replied:

Where does the exponent in the integral come from?

The following is my best try to answer this question.

To sum up, at this point one knows that $\mathrm e^{-a\sqrt{2\lambda}}=\mathrm E(\mathrm e^{-\lambda\tau_a})=\int\limits_0^{+\infty}\mathrm e^{-\lambda t}f_a(t)\mathrm dt$ for some density $f_a$ and the goal is to determine $f_a$. The changes of variables $x=a\sqrt{2\lambda}$ and $t=sa^2$ yield that the function $f$ defined by $f(s)=a^2f_a(a^2s)$ is independent of $a$ and solves $$ \mathrm e^{-x}=\int\limits_0^{+\infty}\mathrm e^{-x^2s/2}f(s)\mathrm ds. $$ My problem here is that I know what is the solution $f$ of this identity, as a distribution canonically associated to the linear Brownian motion, but that what @joriki asked for is a route to find the density $f$. Such a route might be as follows.

First, multiply both sides by $\mathrm e^x$, then complete the square in the exponential. This yields $2x-x^2s=-xz^2+1/s$ with $z=\sqrt{xs}-1/\sqrt{xs}$. The transformation $s\mapsto\zeta(s)=\sqrt{xs}-1/\sqrt{xs}$ is a diffeomorphism from $(0,+\infty)$ to $\mathbb R$ with inverse $z\mapsto\sigma(z)$ and one asks that $$ 1=\int\mathrm e^{-xz^2/2}\mathrm e^{1/(2s)}f(s)\mathrm ds=\int\limits_{-\infty}^{+\infty}\mathrm e^{-xz^2/2}\mathrm e^{1/(2\sigma(z))}f(\sigma(z))\sigma'(z)\mathrm dz. $$ This is a Gaussian integral with respect to $z$ hence, for every nonnegative function $\eta$ such that $\eta(z)+\eta(-z)=1$, the identity $$ \sqrt{x}\eta(z)=2\sqrt{2\pi}\mathrm e^{1/(2\sigma(z))}f(\sigma(z))\sigma'(z), $$ yields a solution $f$. Thus, $$ f(s)=\frac{\sqrt{x}}{\sqrt{2\pi}}\mathrm e^{-1/(2s)}\zeta'(s)2\theta(xs)=\frac1{\sqrt{2\pi}}\mathrm e^{-1/(2s)}\frac1{s\sqrt{s}}(1+xs)\theta(xs), $$ for some function $\theta$ such that $\theta(1/u)+\theta(u)=1$. Choosing $\theta(u)=1/(1+u)$ yields, as required, a solution $f$ independent of $x$, namely, $$ f(s)=\frac1{\sqrt{2\pi}}\mathrm e^{-1/(2s)}\frac1{s\sqrt{s}}. $$