The probability density function of a certain diffusion process in a particular limit.

Question

Let $\mu_0 \in {\mathbb R}$ and let $\theta < 1$. Now, consider a following diffusion process $d X_t = \mu_0 X_t^{2\theta- 1} \cdot dt + X_t^\theta dB_t$ where $B_t$ is a Brownian motion. Here $X_0 = x_0$ is the starting value of the process. Below is a simulation of those processes with $\mu_0,x_0 = 1/3, 1$ and $\theta = -2,-5/4,-1/2,1/4,1$ from violet to red respectively. Here we go:

(*Monte Carlo simulation of the generalized diffusion model.*)
SetOptions[ListPlot, LabelStyle :> {15, FontFamily -> "Arial"}, 
  BaseStyle :> {15, FontFamily -> "Bold"}, PlotMarkers -> Automatic];
NN = 1000;
x = 1; dt = 1/(24 256); mu0 = 1/3;
Bs = Last@RandomVariate[NormalDistribution[], {1, NN}];

(*Generalized diffusion model.*)
ths = Array[# &, 5, {-2, 1}];
X2 = Table[
   Drop[FoldList[#1 + mu0 #1^(2 ths[[ii]] - 1) dt + 
       Sqrt[dt] #1^ths[[ii]] #2 &, x, Bs], -1], {ii, 1, Length[ths]}];
cols = Table[
   ColorData["Rainbow", i/(Length[ths] - 1)], {i, 0, Length[ths] - 1}];
ListPlot[X2, PlotStyle -> cols, Joined -> True, 
 PlotMarkers -> Automatic, AxesLabel -> {"t", "X[t]"}, 
 PlotLabel -> 
  "Monte Carlo simulation. mu0=" <> ToString[N@mu0] <> " x0= " <> 
   ToString[x], PlotLegends -> Automatic]

Now, the probability density of the process $\rho_{X_t} (x,t) := P\left( x \le X_t \le x+dx \right)/dx $ satisfies the Fokker-Planck equation below:

\begin{equation} \frac{\partial }{\partial t} \rho_{X_t}(x,t) = -\frac{\partial }{\partial x} \left[ \mu_0 x^{2\theta-1} \rho_{X_t}(x,t) \right] + \frac{1}{2} \frac{\partial^2}{\partial x^2} \left[ x^{2 \theta} \rho_{X_t}(x,t) \right] \tag{1} \end{equation}

Its solution reads as follows:

\begin{equation} \rho_{X_t}(x,t) = \frac{x_0^{\frac{1}{2}-\mu _0} x^{\mu _0+\frac{1}{2} (1-4 \theta )} e^{-\frac{x^{2-2 \theta }+x_0^{2-2 \theta }}{2 t (1-\theta )^2}} I_{\frac{2 \mu _0-1}{2 \theta -2}}\left(\frac{\left(x x_0\right){}^{1-\theta }}{t (1-\theta )^2}\right)}{t (1-\theta )} \tag{2} \end{equation}

where $I_\nu()$ is the modified Bessel function. Again, the code below verifies that $(2)$ satisfies equation $(1)$ both when $\theta < 1$ and when $\theta = 1$. We have:

In[473]:= 
th =.; mu0 =.; x0 =.; x =.; t =.; Clear[NN];
n = (2 mu0 - 1)/(2 th - 2);
Clear[p]; 
p[x_, t_] :=   (x)^((-4 th + 1)/2 + mu0) (x0)^(1/2 - mu0) 2/(1 - th)
   1/(2 t) Exp[-(2/(1 - th)^2) (x^(2 - 2 th) + x0^(2 - 2 th))/(
    4 t)] BesselI[n, 2/(1 - th)^2 (x x0)^(1 - th)/(2 t)];
Clear[plim]; 
plim[x_, t_] := 
 1/x 1/Sqrt[2 Pi t] Exp[-(Log[x/x0] - (mu0 - 1/2) t)^2/(2 t)]
(*Check if forward Kolmogorov equation satisfied.*)
rhs = (-D[ mu0 x^(2 th - 1) p[x, t], x] + 
    1/2 D[x^(2 th) p[x, t], {x, 2}] - D[p[x, t], t]);
rhs // PowerExpand // FullSimplify

rhs = (-D[ mu0 x plim[x, t], x] + 1/2 D[x^(2 ) plim[x, t], {x, 2}] - 
    D[plim[x, t], t]);
rhs // PowerExpand // FullSimplify


Out[478]= 0

Out[480]= 0

Those probability densities in equation $(2)$ are being plotted as a function of time $t= 1/50 + j/50$ for $j=0,\cdots, 49$ (from violet to red respectively) below:

th =.; mu0 =.; x0 =.; x =.; t =.;
Clear[p]; 
p[x_, t_] :=   (x)^((-4 th + 1)/2 + mu0) (x0)^(1/2 - mu0) 2/(1 - th)
   1/(2 t) Exp[-(2/(1 - th)^2) (x^(2 - 2 th) + x0^(2 - 2 th))/(
    4 t)] BesselI[n, 2/(1 - th)^2 (x x0)^(1 - th)/(2 t)];
Clear[NN]; 
NN[t_] := (1 - 
   GammaRegularized[(-1 + 2 mu0)/(-2 + 2 th), x0^(2 - 2 th)/(
    2 t (-1 + th)^2)]);

(*Choose parameters .*)
{mu0, th, x0} = {1/3, -2, 1};
(*Plot*)
ts = Array[# &, 50, {1/50, 1}];
Clear[cols]; 
cols[m_] := Table[ColorData["Rainbow", i/(m - 1)], {i, 0, m - 1}];
xs = Array[# &, 500, {1/100, 3}];
pl1 = ListPlot[
  Table[Transpose[{xs, (p[#, ts[[i]]]/NN[ts[[i]]]) & /@ xs}], {i, 1, 
    Length[ts]}], PlotRange :> All, PlotStyle :> cols[Length[ts]], 
  AxesLabel -> {"x", "pdf[x]"}, 
  PlotLabel -> "mu0,th,x0=" <> ToString[N@{mu0, th, x0}], 
  ImageSize :> 400]

---

Clearly as $\theta \rightarrow 1_-$ the process in question goes into the geometric Brownian motion and as such the following limit must hold true:

\begin{equation} \lim\limits_{\theta \rightarrow 1_-} \rho_{X_t}(x,t) = \frac{1}{x} \frac{1}{\sqrt{2\pi t}} \cdot \exp\left( -\frac{1}{2 t} \left[ \log(\frac{x}{x_0}) - (\mu_0- \frac{1}{2}) \cdot t\right]^2\right) \tag{3} \end{equation}

My problem is that I am unable to verify this limit. I did the following calculations:

\begin{eqnarray} &&\lim\limits_{\theta \rightarrow 1_-} \rho_{X_t}(x,t) = \\ &&\frac{x_0^{\frac{1}{2}-\mu _0} x^{\mu _0+\frac{1}{2} (1-4 \theta )} e^{-\frac{x^{2-2 \theta }+x_0^{2-2 \theta }}{2 t (1-\theta )^2}} }{t (1-\theta )} \cdot % \frac{\exp\left(\frac{\left(x x_0\right){}^{1-\theta }}{t (1-\theta )^2}\right)}{\sqrt{2\pi \frac{\left(x x_0\right){}^{1-\theta }}{t (1-\theta )^2} }} =\\ && \left( \frac{x}{x_0}\right)^{\mu_0 -\frac{1}{2}} \cdot \frac{1}{\sqrt{t}} \cdot x^{1-2\theta} \cdot e^{-\frac{(x^{1- \theta }+x_0^{1- \theta })^2}{2 t (1-\theta )^2}} \cdot \frac{1}{\sqrt{2 \pi (x_0 x)^{1-\theta}}} \underbrace{=}_{\theta \rightarrow 1_-} \\ &&e^{(\mu_0-\frac{1}{2}) \log(\frac{x}{x_0})} \cdot \frac{1}{\sqrt{2 \pi t}} \cdot \frac{1}{x} \cdot \exp\left( -\frac{1}{2 t} \left[ \log(\frac{x}{x_0})\right]^2\right) =\\ && \frac{1}{\sqrt{2 \pi t}} \cdot \frac{1}{x} \cdot \exp\left( -\frac{1}{2 t} \left[ \log(\frac{x}{x_0}) - (\mu_0 - \frac{1}{2}) \cdot t \right]^2 \right) \cdot e^{\frac{1}{2} (\mu_0-\frac{1}{2})^2 t} \tag{4} \end{eqnarray}

In the second line I used the asymptotic form of the modified Bessel function from here, in the third line I merged the two exponentials together and then simplified the result, in the forth line I took the limit $\theta \rightarrow 1_-$ and then in the fifth line I simplified everything. Clearly the result is different from the right hand side of equation $(3)$. Where is the error?

Answer 1

Assume $0<\mu_0 <\frac{1}{2}$ is fixed. Note that $$ I_{\frac{{2\mu _0 - 1}}{{2\theta - 2}}} \!\!\left( {\frac{{(xx_0 )^{1 - \theta } }}{{t(1 - \theta )^2 }}} \right) = I_{-\left( {\frac{1}{2}-\mu _0} \right)\frac{1}{{\theta - 1}}} \!\!\left( {\left( {\frac{1}{2}-\mu _0} \right)^{ - 2} \frac{{(xx_0 )^{1 - \theta } }}{t}\left( {\frac{1}{2}-\mu _0} \right)^2\frac{1}{{(\theta - 1)^2 }}} \right). $$ Therefore, you need to use an asymptotics for $I_{-\nu} (a\nu ^2 )$ as $\nu \to +\infty$ with $a>0$. By the uniform asymptotic expansions of the modified Bessel functions (take $z=a\nu$) $$ I_{ - \nu } (a\nu ^2 ) = I_\nu (a\nu ^2 ) + \frac{2}{\pi }\sin (\pi \nu )K_\nu (a\nu ^2 ) \sim \frac{{\exp (a\nu ^2 - 1/(2a))}}{{\sqrt {2\pi a} \nu }} $$ as $\nu\to+\infty$ with any $a>0$. The asymptotics $$ I_\alpha(x) \sim \frac{\mathrm{e}^x}{\sqrt{2\pi x}} \quad (x\to +\infty), $$ that you used, is valid only if $\alpha^2=o(x)$. However, in your case $\alpha^2=\Theta(x)$. Hence the appearance of the extra term $-1/(2a)$ in the exponential.