# The distribution of the first hitting time for the Constant Elasticity of Variance process.

The Constant Elasticity of Variance (CEV) process is a one dimensional diffusion process given by the following stochastic differential equation.

$$$$d X_t = \mu X_t \cdot dt + \sigma X_t^\beta \cdot d B_t \tag{1}$$$$

where $$\mu,\sigma,\beta$$ are positive real parameters and $$B_t$$ is the Brownian motion. In what follows we assume that the starting value of the process reads $$X_0 = x > 0$$ and that $$\beta > 1$$. The infinitesimal generator of this process reads $${\mathfrak G}_z := \mu z d/d z + \sigma^2/2 z^{2 \beta} d^2/d z^2$$ . The eigen-functions of this operator $${\mathfrak G}_z \phi^{\pm} (z) = \lambda \phi^{\pm}(z)$$ to the eigenvalue $$\lambda > 0$$ are given below:

$$\begin{eqnarray} \phi^{+}(z) &=&U\left(\frac{\lambda}{2(-1+\beta) \mu}, 1+ \frac{1}{2(-1+\beta)}, \frac{\mu z^{2-2\beta} }{(-1+\beta) \sigma^2} \right) \tag{2a} \\ \phi^{-}(z) &=& L_{-\frac{\lambda}{2(-1+\beta) \mu}}^{(\frac{1}{2(-1+\beta)})} \left( \frac{\mu z^{2-2\beta} }{(-1+\beta) \sigma^2}\right) \tag{2b} \end{eqnarray}$$

where $$U$$ is the confluent hypergeometric function and $$L$$ is the generalized Laguerre polynomial. We have checked that $$U(z)$$ is a strictly increasing function of $$z$$.

Now, by using the theory of diffusion processes, see section 4.6 pages 128-134 in

Ito, K.; McKean, H. P. jun., Diffusion processes and their sample paths, Berlin-Heidelberg-New York: Springer-Verlag. XVII, 321 p. (1965). ZBL0127.09503.

, we have found the Laplace transform of the first hitting time $$\tau_y := inf(s>0:X_s = y)$$ of a horizontal barrier $$b$$ by this process. The quantity in question reads:

$$\begin{eqnarray} E_x \left[ e^{-\lambda \tau_y} \right] = \frac{\phi^{+}(x)}{\phi^{+}(y)} \quad \mbox{for x \le y} \tag{3} \end{eqnarray}$$

Now, by inverting the Laplace transform in $$(3)$$, by using the Bromwich integral and then the Cauchy theorem, we have expressed the probability density function of the first hitting time $$n_x(t;y) := P_x\left( \tau_y \in dt\right)/dt$$ as follows:

$$\begin{eqnarray} n_x(t;y) = \sum\limits_{p=1}^\infty \underbrace{ \frac{U\left( -\zeta_p^{(y;\mu,\sigma,\beta)}, 1+ \frac{1}{2(-1+\beta)}, \frac{\mu x^{2-2\beta}}{(-1+\beta) \sigma^2} \right)}{ \zeta_p^{(y;\mu,\sigma,\beta)} U^{(1,0,0)}\left( -\zeta_p^{(y;\mu,\sigma,\beta)}, 1+ \frac{1}{2(-1+\beta)}, \frac{\mu y^{2-2\beta}}{(-1+\beta) \sigma^2} \right) } }_{{\mathfrak w}_p^{(y;\mu,\sigma,\beta)}} \cdot \underline{(2 (-1+\beta) \mu \zeta_p^{(y;\mu,\sigma,\beta)} ) \cdot e^{-2 (-1+\beta) \mu \cdot \zeta_p^{(y;\mu,\sigma,\beta)} \cdot t}} \tag{5} \end{eqnarray}$$

As we can see the quantity in $$(5)$$ is an infinite linear combination of exponential distributions with weights $$\left( {\mathfrak w}_p^{(y;\mu,\sigma,\beta)} \right)_{p=1}^\infty$$ that sum up to unity. Here $$\left( \zeta_p^{(y;\mu,\sigma,\beta)} \right)_{p=1}^\infty$$ are zeros of the function $${\mathbb R}_+ \ni \lambda \rightarrow U(-\lambda,1+ \frac{1}{2(-1+\beta)}, \frac{\mu y^{2-2\beta}}{(-1+\beta) \sigma^2}) \in {\mathbb R}$$.

Now we took the following process parameters $$\mu,\sigma,\beta = 3/2,1/2,5/2$$ and the first value and the barrier $$x,y = 3/2, 5/2$$ and we plotted the quantity $$(5)$$ below. We also verified the normalization numerically. Here we go:

{\[Mu], \[Sigma], \[Beta]} = {3/2, 1/2, 5/2};
(*Here x\[LessEqual]y*)
{x, y} = {3/2, 5/2};

SetOptions[FindRoot, WorkingPrecision -> mprec, PrecisionGoal -> prec];
mzeros = \[Lambda] /.
Table[FindRoot[
HypergeometricU[-\[Lambda],
1 + 1/(2 (-1 + \[Beta])), (\[Mu] y^(
2 - 2 \[Beta]))/((-1 + \[Beta]) \[Sigma]^2)] == 0, {\[Lambda],
n}], {n, 0, 50}];
mzeros = Sort[#[[1]] & /@ Tally[mzeros]];

ts = Array[# &, {300}, {1/100, 3}];
vals = {#,
Total[Table[
HypergeometricU[-mzeros[[p]],
1 + 1/(2 (-1 + \[Beta])), (\[Mu] x^(
2 - 2 \[Beta]))/((-1 + \[Beta]) \[Sigma]^2)]/
\!$$\*SuperscriptBox[\(HypergeometricU$$,
TagBox[
RowBox[{"(",
RowBox[{"1", ",", "0", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[-mzeros[[p]],
1 + 1/(2 (-1 + \[Beta])), (\[Mu] y^(
2 - 2 \[Beta]))/((-1 + \[Beta]) \[Sigma]^2)] (2 (-1 + \
\[Beta]) \[Mu]) Exp[-2 (-1 + \[Beta]) \[Mu] mzeros[[p]] #], {p, 1,
Length[mzeros]}]]} & /@ ts;
ListPlot[vals, PlotRange :> All,
AxesLabel -> {"t", "\!$$\*SubscriptBox[\(n$$, $$x$$]\)(t,y)"}]
f = Interpolation[vals];
NIntegrate[f[xi], {xi, 0.01, 3}]


As you can see the distribution in question has a correct shape and a correct normalization. The negative values close to the origin are due to a numerical error.

Having said all this my question would be how do we evaluate the limit of $$\beta \rightarrow 1_+$$. In this case the process tends towards the geometric Brownian motion and as such we should have:

$$$$\lim_{\beta \rightarrow 1_+} n_x(t;y) \stackrel{(??)}{=} \frac{\left| \log(\frac{y}{x} )\right|}{\sqrt{2 \pi t^3} \sigma} e^{-\frac{1}{2 \sigma^2 t} \left[ \log(\frac{y}{x} - (\mu - \frac{\sigma^2}{2} ) t\right]^2}$$$$

as shown in a previous question on a similar topic.

How do we work out this limit analytically in our framework?

##### Update:

We have verified numerically that the Laplace transform $$(3)$$ approaches the correct limits when $$\beta \rightarrow 1_+$$. See code below:

    {lmb, mu, sig} = RandomReal[{0, 1}, 3, WorkingPrecision -> 50];
x = RandomReal[{0, 1}, WorkingPrecision -> 50];
y = RandomReal[{x, 2}, WorkingPrecision -> 50];
NN = 100;
HypergeometricU[lmb/mu NN, NN, NN (2 mu/sig^2) x^(-1/NN)]/
HypergeometricU[lmb/mu NN, NN, NN (2 mu/sig^2) y^(-1/NN)]
(x/y)^((-2 mu + sig^2 + 2 Sqrt[2 lmb sig^2 + (mu - sig^2/2)^2])/(
2 sig^2))

0.275584165705622319278498109409236772998453 + 0.*10^-50 I

0.2729011283963510407220223125184479965801097988605


We are going to prove that the Laplace transform $$(3)$$ approaches the correct limit when $$\beta \rightarrow 1_+$$. We denote $${\mathfrak N}:= 1/(2(-1+\beta)$$ and $$\theta := 2 \mu/\sigma^2$$ and we have:

$$\begin{eqnarray} &&\left( \frac{x}{y} \right)^{-\frac{\lambda}{\mu}}E_x \left[ e^{-\tau_y} \right] = \\ &&\left( \frac{x}{y} \right)^{-\frac{\lambda}{\mu}} \frac{U(\frac{\lambda}{\mu} {\mathfrak N}, {\mathfrak N}, {\mathfrak N} \theta x^{-1/{\mathfrak N}})}{ U(\frac{\lambda}{\mu} {\mathfrak N}, {\mathfrak N}, {\mathfrak N} \theta y^{-1/{\mathfrak N}})} = \\ % && \frac{ \int\limits_0^\infty e^{-t} t^{\frac{\lambda}{\mu} {\mathfrak N}-1} \left( 1 + t \frac{x^{1/{\mathfrak N}}}{{\mathfrak N} \cdot \theta} \right)^{{\mathfrak N}(1- \lambda/\mu-1/{\mathfrak N})} dt } { \int\limits_0^\infty e^{-t} t^{\frac{\lambda}{\mu} {\mathfrak N}-1} \left( 1 + t \frac{y^{1/{\mathfrak N}}}{{\mathfrak N} \cdot \theta} \right)^{{\mathfrak N}(1- \lambda/\mu-1/{\mathfrak N})} dt } =\\ && \left( \frac{x}{y} \right)^{-\frac{\lambda}{\mu}} \cdot \frac{ % \int\limits_0^1 \left( e^{- \frac{\theta (1-v) x^{-1/{\mathfrak N}}}{v}} \cdot \frac{(1-v)^{\frac{\lambda}{\mu}}}{v} \right)^{\mathfrak N} \cdot (1-v)^{-1} dv % } { % \int\limits_0^1 \left( e^{- \frac{\theta (1-v) y^{-1/{\mathfrak N}}}{v}} \cdot \frac{(1-v)^{\frac{\lambda}{\mu}}}{v} \right)^{\mathfrak N} \cdot (1-v)^{-1} dv % } = \\ && \left( \frac{x}{y} \right)^{-\frac{\lambda}{\mu}} \cdot % \frac{ % \exp(-{\mathfrak N} f(x,v^{(*)})) \sqrt{\frac{2 \pi}{{\mathfrak N} f(x,v^{(*)})} } % } { % \exp(-{\mathfrak N} f(y,v^{(*)})) \sqrt{\frac{2 \pi}{{\mathfrak N} f(y,v^{(*)})} } % } \stackrel{=}{{\mathfrak N} \rightarrow \infty} \\ % && \left( \frac{x}{y} \right)^{-\frac{\lambda}{\mu}} \cdot % \frac{ % \exp(-{\mathfrak N} f(x,v^{(*)})) % } { % \exp(-{\mathfrak N} f(y,v^{(*)})) % } = \\ && \left( \frac{x}{y} \right)^{-\frac{\lambda}{\mu}} \cdot % \exp \left(\frac{1}{2} \left(-\frac{\sqrt{4 \theta \lambda +(\theta -1)^2 \mu }}{\sqrt{\mu }}+\theta -1\right) {\mathfrak N} \cdot \left(x^{-1/{\mathfrak N}}-y^{-1/{\mathfrak N}}\right)\right) \stackrel{=}{{\mathfrak N} \rightarrow \infty} \\ % && \left(\frac{x}{y}\right)^{\sqrt{\frac{\theta \lambda }{\mu }+\left(\frac{1}{2}-\frac{\theta }{2}\right)^2}-\frac{\theta }{2}-\frac{\lambda }{\mu }+\frac{1}{2}} \end{eqnarray}$$

Now, in the third line from the top we used the integral representation of the confluent hypergeometric function . In the fourth line from the top we substituted for $$v := (1+ t x^{1/{\mathfrak N}}/({\mathfrak N} \theta))^{-1}$$. In the fifth line we applied the stationary phase approximation with $$f(x,v):= \theta (1-v)/v x^{-1/{\mathfrak N}} - \log((1-v)^{\lambda/\mu}/v)$$. Here $$v^{(*)} := \frac{\sqrt{\mu } \sqrt{\theta ^2 \mu +4 \theta \lambda -2 \theta \mu +\mu }-\theta \mu -\mu }{2 (\lambda -\mu )}$$ and it is the stationary point of both integrals, i.e. that in the numerator and in the denominator, meaning $$d f(x,v)/dv = d f(y, v)/d v = 0$$. In the sixth line we dropped the higher order terms in $${\mathfrak N}$$. In the next line we evaluated the expression from the line directly above. In the final line we took the limit $${\mathfrak N} \rightarrow \infty$$.

As always, we include a Mathematica code snippet that verifies all the steps taken:

prec = 50;
NN = 100 + 1/10;
SetOptions[NIntegrate, WorkingPrecision -> 15, PrecisionGoal -> 10];
{lmb, mu} = RandomReal[{0, 1}, 2, WorkingPrecision -> prec];
th = RandomReal[{1, 2}, WorkingPrecision -> prec];
x = RandomReal[{0, 1}, WorkingPrecision -> prec];
y = RandomReal[{x, 2}, WorkingPrecision -> prec];
(x/y)^(-lmb/mu) HypergeometricU[lmb/mu NN, NN, NN (th) x^(-1/NN)]/
HypergeometricU[lmb/mu NN, NN, NN (th) y^(-1/NN)]

(*Use the integral representation from \
https://en.wikipedia.org/wiki/Confluent_hypergeometric_function#\
Integral_representations and substitute in it for z t.*)
NIntegrate[
Exp[-t] t^(lmb/mu NN - 1) (1 + t/(NN th ) x^(1/NN))^(
NN (1 - lmb/mu - 1/NN)), {t, 0, Infinity}]/
NIntegrate[
Exp[-t] t^(lmb/mu NN - 1) (1 + t/(NN th ) y^(1/NN))^(
NN (1 - lmb/mu - 1/NN)), {t, 0, Infinity}]

(*Substitute in the above for v= (1+ t/(NN th )x^(1/NN))^(-1)*)
(x/y)^(-lmb/
mu) NIntegrate[    (E^(-((th (1 - v) x^(-1/NN))/v)) (1 - v)^(lmb/
mu)/v)^NN (1 - v)^-1, {v, 0, 1}]/
NIntegrate[  (E^(-((th (1 - v) y^(-1/NN))/v)) (1 - v)^(lmb/mu)/v)^
NN (1 - v)^-1 , {v, 0, 1}]

(*Apply the stationary point approximation*)
v =.; f[x_, v_] := (th (1 - v) x^(-1/NN))/v - Log[(1 - v)^(lmb/mu)/v];
vstar1 = (v /. FindRoot[D[f[x, v], v] == 0, {v, 1/2}]);
vstar2 = (v /. FindRoot[D[f[y, v], v] == 0, {v, 1/2}]);
vstar = (-mu - mu th +
Sqrt[mu] Sqrt[mu + 4 lmb th - 2 mu th + mu th^2])/(2 (lmb - mu));

(*
num=Exp[-NN f[x,vstar]]Sqrt[2 Pi/(NN f[x,vstar])];
denom=Exp[-NN f[y,vstar]]Sqrt[2 Pi/(NN f[y,vstar])];
*)
num = Exp[-NN f[x, vstar]];
denom = Exp[-NN f[y, vstar]];
(x/y)^(-lmb/mu) num/denom

(x/y)^(-lmb/mu) Exp[
1/2 (-1 + th - Sqrt[mu (-1 + th)^2 + 4 lmb th]/Sqrt[
mu]) NN (x^(-1/NN) - y^(-1/NN))]

(x/y)^((-1/2 th + 1/2 + Sqrt[ lmb th/(mu) + (-1/2 th + 1/2)^2]) -
lmb/mu)