关于广义扩散过程概率密度函数θ→1⁻极限验证的错误排查问询
问题背景
给定$\mu_0 \in {\mathbb R}$且$\theta < 1$，考虑如下扩散过程：
$$d X_t = \mu_0 X_t^{2\theta- 1} \cdot dt + X_t^\theta dB_t$$
其中$B_t$是布朗运动，初始值$X_0 = x_0$。
蒙特卡洛模拟实现
下面是参数$\mu_0,x_0 = 1/3, 1$，$\theta = -2,-5/4,-1/2,1/4,1$（曲线颜色从紫到红）的蒙特卡洛模拟代码：
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]
福克-普朗克方程及解析解
过程的概率密度$\rho_{X_t} (x,t) := P\left( x \le X_t \le x+dx \right)/dx$满足如下福克-普朗克方程：
$$
\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}
$$
其解析解为：
$$
\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}
$$
其中$I\nu()$是修正贝塞尔函数。
解析解的验证代码
以下代码验证式(2)在$\theta < 1$和$\theta = 1$时均满足方程(1)：
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
概率密度随时间演化的绘制代码
下面是式(2)的概率密度随时间$t= 1/50 + j/50$（$j=0,\cdots, 49$，曲线颜色从紫到红）的绘制代码：
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]
极限验证的疑惑
当$\theta \rightarrow 1_-$时，该过程收敛到几何布朗运动，因此理论上应有如下极限：
$$
\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}
$$
但我自行推导验证这个极限时，得到了与式(3)不一致的结果，我的推导过程如下：
$$
\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}
$$
推导步骤说明：
- 第二行使用了修正贝塞尔函数的渐近形式
- 第三行合并指数项并进行代数化简
- 第四行取$\theta \rightarrow 1_-$的极限
- 第五行对结果做进一步整理
显然式(4)与式(3)的结果存在差异，请问我的推导过程中哪里出现了错误？
备注：内容来源于stack exchange，提问作者Przemo
