# Transition probability density function for a non-trivial diffusion process.

Let $$\mu$$ and $$\sigma > 0$$ and $$\beta_1 \ge 0$$ and $$\beta_2 \ge 0$$ be real numbers. Consider a stochastic process $$X_t$$ that satisfies the following stochastic differential equation:

$$$$d X_t = \mu X_t^{\beta_1} dt + \sigma X_t^{\beta_2} dB_t \tag{1}$$$$ where $$B_t$$ is the Brownian motion.

My objective is to derive the transition probability density function $$P_x\left( X_t \in dz\right) := P\left( X_t \in (z,z+dz) | X_0 = x \right)$$ for the process in question. Here we will be using a result known from literature (see section 4.11 pages 149-161 in K Ito& H P McKean Jr, "Diffusion Processes and their Sample Paths") as below.

##### Theorem:

The Laplace transform (with respect to time) of the transition probability density in question reads:

$$$$\int\limits_0^\infty e^{-\lambda t} P_x\left( X_t \in dz\right) dt = G_\lambda(z,x) \cdot m(dz) \tag{2}$$$$

Here $$m(z):= \left( 1/2 \sigma^2(z) S(z) \right)^{-1}$$ is the speed measure with $$S(y) := \exp\left(- \int\limits_1^y 2 \mu(\xi)/\sigma^2(\xi) d\xi \right)$$ being the scale measure.

Then $$G_\lambda​(z,x)$$ is the Green's function of the operator $$\lambda - {\mathfrak G}^{*}$$ with $${\mathfrak G}^{*}$$ being the infinitesimal generator, i.e. such an operator that its action on a trial function $$f$$ reads $$\left({\mathfrak G}^{*} f\right)(x) = \lim\limits_{\epsilon \rightarrow 0} E \left[ f(X_\epsilon) - f(x) \right]/\epsilon$$ subject to $$X_0=x$$.

Note that the Green's function is easily determined by the eigen-functions and it reads $$G_\lambda(z,x) = S(x)/{\mathfrak W}_\lambda(x) \left( G_\lambda(x) F_\lambda(z) 1_{z < x} + F_\lambda(x) G_\lambda(z) 1_{z > x} \right)$$ where $$F_\lambda(z)$$ and $$G_\lambda(z)$$ is a strictly increasing and a strictly decreasing eigen-function of the infinitesimal generator to an eigenvalue $$\lambda$$ and $${\mathfrak W}_\lambda(x) := F_\lambda^{'}(x) G_\lambda(x) - G_\lambda^{'}(x) F_\lambda(x)$$ is the Wronskian of the generator in question.

##### What have I managed to achieve?:

I have found the result for $$(\beta_1,\beta_2) = (1, \beta)$$ with $$\beta >1$$ and I present the derivation below.

Here the infinitesimal generator reads $${\mathfrak G}_z := \mu z \cdot d/dz + 1/2 \sigma^2 z^{2 \beta} d^2/ d z^2$$ and its eigenfunctions (compare my previous question on a similar topic) read:

$$\begin{eqnarray} F(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{3a} \\ G(z) &=& F_{1,1}\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{3b} \end{eqnarray}$$ where $$U$$ is the confluent hypergeometric function.

Note that $$F^{'}(z) = \frac{\lambda z^{1-2 \beta } U\left(\frac{\lambda }{2 (\beta -1) \mu }+1,2+\frac{1}{2 (\beta -1)},\frac{z^{2-2 \beta } \mu }{s^2 (\beta -1)}\right)}{(\beta -1) s^2} >0$$ and $$G^{'}(z) = -\frac{2 \lambda z^{1-2 \beta } \, _1F_1\left(\frac{\lambda }{2 (\beta -1) \mu }+1;2+\frac{1}{2 (\beta -1)};\frac{z^{2-2 \beta } \mu }{s^2 (\beta -1)}\right)}{(2 \beta -1) s^2} < 0$$ (note that from the integral representations we see that both $$U$$ and $$F_{1,1}$$ are positive --well up to a multiplicative constant in the second case) and as such $$F(z)$$ and $$G(z)$$ are strictly increasing and strictly decreasing as it should be.

Now the Wronskian is computed as follows:

$$\begin{eqnarray} {\mathfrak W}_\lambda(x) &=& {\mathfrak W}_\lambda(1) \cdot \exp\left( -\int\limits_1^x \frac{2\mu(\xi)}{\sigma^2(\xi)} d\xi \right) \tag{4a} \\ &=& {\mathfrak W}_\lambda(1) \cdot \exp\left( - \frac{2 \mu}{\sigma^2} \frac{x^{2-2\beta} -1}{2-2\beta}\right) \tag{4b} \\ &=& -\frac{2 (\beta -1) \Gamma \left(1+\frac{1}{2 (\beta -1)}\right) e^{\frac{\mu }{(\beta -1) \sigma ^2}} \left(\frac{\mu }{(\beta -1) \sigma ^2}\right)^{\frac{1}{2-2 \beta }}}{\Gamma \left(\frac{\lambda }{(2 \beta -2) \mu }\right)} \cdot \exp\left( - \frac{2 \mu}{\sigma^2} \frac{x^{2-2\beta} -1}{2-2\beta}\right) \tag{4c} \end{eqnarray}$$ The first line above simply follows form the first order differential equation being satisfied by the Wronskian. The second line follows form the fact that $$\mu(\xi), \sigma(\xi) = \mu \cdot \xi, \sigma \cdot \xi^\beta$$. The step to the final line is the most laborious, but not really that mind boggling. Here one employs the integral representations in order to express the Wronskian at unity as a double integral in which one changes variables to obtain the neat final result.

Now we compute the inverse Laplace transform as the Bromwich integral and we use Cauchy theorem (being applied to the usual contour in the negative complex half-plane-- the second & the third quadrants) to evaluate that integral. As seen from $$(4c)$$ the Laplace transform has infinitely many simple poles at $$\lambda_n = (2\beta-2) \mu (-n)$$ for $$n=0,1,2,\cdots$$. Then the final result reads:

$$\begin{eqnarray} &&\left. P_x \left( X_t \in dz \right)/dz = \frac{2 \mu \left(\frac{\mu }{(\beta -1) \sigma ^2}\right)^{\frac{1}{2 (\beta -1)}}}{\sigma ^2 \Gamma \left(1+\frac{1}{2 (\beta -1)}\right)} \cdot z^{-2 \beta } e^{-\frac{\mu z^{2-2 \beta }}{(\beta -1) \sigma ^2}} \cdot \right.\\ && \left. \sum\limits_{n=0}^\infty \frac{(-1)^n}{n!} % \, _1F_1\left(-n;1+\frac{1}{2 (\beta -1)};\frac{\mu (x \vee z)^{2-2 \beta }}{(\beta -1) \sigma ^2}\right) U\left(-n,1+\frac{1}{2 (\beta -1)},\frac{\mu (x \wedge z)^{2-2 \beta }}{(\beta -1) \sigma ^2}\right) % \cdot e^{-2(-1+\beta) \cdot \mu \cdot n \cdot t} \right. \tag{5} \end{eqnarray}$$

Now we have checked numerically that the result above matches with the solution known from literature (see equations (4) and (5) in V Linetsky, R Mendoza, The Constant Elasticity of Variance Model). Here we go:

In[1991]:= (*Verify whether our solution matches that known from \
literature?*)
Clear[mu, s, b, z, lmb, phi, v, W, W1, W0, x];
{mu, s, b} = {1/2, 1/3 + 3/10, 5/2};
{x, z} = RandomReal[{1, 3}, 2, WorkingPrecision -> 50];
t = 3/2;
(*Our solution:*)
res1 = (2  mu ((mu/((-1 + b) s^2))^((1/(2 (-1 + b))))) )/(
s^2 Gamma[1 + 1/(2 (-1 + b))])
E^(-((mu z^(2 - 2 b))/((-1 + b) s^2))) z^(-2 b)
Table[(-1)^n/
n! HypergeometricU[-n, 1 + 1/(2 (-1 + b)),
mu/((-1 + b) s^2) Min[x, z]^(2 - 2 b)] Hypergeometric1F1[-n,
1 + 1/(2 (-1 + b)),
mu/((-1 + b) s^2) Max[x, z]^(2 - 2 b)] E^(-
2 (-1 + b) mu n t), {n, 0, M}];
Total[res1] // N

(*Solution known from literature: Equations (4) & (5) page 2 in V \
Linetsky & R Mendoza The Constant Elasticity of Variance Model.Here \
the solution was obtained,using different techniques,meaning by \
re-scaling the current price and by changing time in the solution of \
a simpler problem,i.e.a problem with the drift parameter being set to \
zero.*)

Clear[tau]; tau[t_] := (Exp[2 mu (b - 1) t] - 1)/(2 mu (b - 1));
res2 = Exp[-mu t] ((Exp[-mu t] z)^(-2 b + 1/2) x^(1/2))/(
s^2 Abs[b - 1] tau[t])
BesselI[1/(2 Abs[ b - 1]), (x^(-b + 1) (Exp[-mu t] z)^(-b + 1))/(
s^2 (b - 1)^2 tau[
t])] Exp[-((x^(-2 b + 2) + (Exp[-mu t] z)^(-2 b + 2))/(
2 s^2 (b - 1)^2 tau[t]))];
res2 // N

Out[1996]= 0.0481819

Out[1999]= 0.0481819


Having said all this my question is how do we find the solution in the generic case, i.e. when $$\beta_1 \neq 1$$? Another obvious question would be to provide a literature reference, if it exists, on that problem solution.

Here we will only focus on finding the eigenfunctions of the infinitesimal generator. In other words we want to solve the following ODE below:

$$$$\left(\mu z^{\beta_1} \frac{d}{d z} + \frac{\sigma^2}{2} z^{2 \beta_2} \frac{d^2}{d z^2} - \lambda\right) f(z) =0 \tag{1}$$$$

Note that even this problem is very hard in general. Indeed, if we do the usual trick, i.e. change the ordinate in such a way in order to annihilate the coefficient at the first derivative by writing $$f(z):= m(z) v(z)$$ where , in this case, $$m(z) := e^{-\frac{\mu z^{\beta _1-2 \beta _2+1}}{\left(\beta _1-2 \beta _2+1\right) \sigma ^2}}$$ then the resulting equation for the function $$v(z)$$ takes the form: $$\frac{v(z) z^{-2 \beta _2} \left(\mu z^{\beta _1} \left(-\frac{\mu z^{\beta _1-2 \beta _2}}{\sigma ^2}-\frac{\beta _1-2 \beta _2}{z}\right)-2 \lambda\right)}{\sigma ^2}+v''(z)=0$$ which leaves us hopeless because in the relevant section in

Polyanin, A. D.; Zaitsev, Valentin F., Handbook of exact solutions for ordinary differential equations., Boca Raton, FL: CRC Press. xxvi, 787 p. (2003). ZBL1015.34001.

, meaning in section 2.1.2-1 there is only one entry (entry 10) that is vaguely reminiscent of what we have above (it is only in very special cases that the above reduces to that entry).

In view of this we have to analyze the problem from scratch. The first observation we make is that the solution to $$(1)$$ can be always written as a double infinite series as follows $$f(z) := \sum\limits_{n=0}^\infty \sum\limits_{m=0}^\infty a_{n,m} z^{\alpha + (1-\beta_1) n + (1-\beta_2) m}$$ for some coefficients $$( a_{n,m} )_{(n,m) \in {\mathbb N}^2}$$. It is fairly easy to write down the recurrence relations for those coefficients but, as it turns out, it is not easy to solve those relations (except numerically but this is not that interesting because the solution cannot be used to to construct our transition densities). However the functional form of the conjectured solution suggests that things might simplify if either $$(1-\beta_1) = \theta \cdot (1-\beta_2)$$ or $$(1-\beta_2) = \theta \cdot (1-\beta_1)$$ where $$\theta \in {\mathbb Z}$$. This is what we are going to analyze now.

##### First case:

We assume that $$(1-\beta_1) = \theta \cdot (1-\beta_2)$$ where $$\theta \in {\mathbb Z}$$ and we substitute for $$y = A z^{1-\beta_2}$$ in $$(1)$$. This leads to a following ODE below:

$$$$\frac{1}{2} A^2 \left(\beta _2-1\right){}^2 \sigma ^2 f''(y)+\frac{1}{2} \left(\beta _2-1\right) A^{-\frac{2}{\beta _2-1}} y^{-\theta -1} f'(y) \left(\beta _2 \sigma ^2 A^{\frac{2 \beta _2}{\beta _2-1}} y^{\theta }-2 \mu y^2 A^{\frac{2}{\beta _2-1}+\theta }\right)-\lambda f(y)=0 \tag{2}$$$$

Note that the ODE above is simpler because it has polynomial coefficients in the abscissa .

See code snippet for the derivation of $$(2)$$:

In[874]:= Clear[mu, s, lmb, b1, b2, y, z, xi, f, th, A];
(*Exact eigenfunctions of the infinitesimal generator.*)
(*First case: (1-b1) = th(1-b2)  with th in integers.*)
zz[y_] = z /. First@Solve[y == A z^(1 - b2), z];
u[y_] := D[zz[y], y];
ex = Nest[(mu z^b1 D[#, z] + s^2/2 z^(2 b2) D[#, {z, 2}] - lmb #) &,
f[z], 1];
ex = ex /. Derivative[1][f][z] :> 1/u[y] D[f[y], y] /.
Derivative[2][f][z] :> 1/u[y] D[ 1/u[y] D[f[y], y], y] /.
f[z] :> f[y] /. z :>  zz[y];
ex = PowerExpand[ex] /. b1 :> th (b2 - 1) + 1;
ex = Collect[ex, Derivative[n_][f][y], Simplify[#] &];
ex = Expand[ex] /. y^n_ :> y^Simplify[n];
(*The later is then exactly solved in case th=0,1,2.*)
ex = Collect[ex, Derivative[n_][f][y], Simplify[PowerExpand[#]] &]

During evaluation of In[874]:= Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information.

Out[882]= -lmb f[y] +
1/2 A^(-(2/(-1 + b2))) (-1 + b2) y^(-1 -
th) (-2 A^(2/(-1 + b2) + th) mu y^2 +
A^((2 b2)/(-1 + b2)) b2 s^2 y^th) Derivative[1][f][y] +
1/2 A^2 (-1 + b2)^2 s^2 (f^\[Prime]\[Prime])[y]


Now, we used Mathematica find a clue as to how those solutions look like. We found the following three exactly solvable cases:

1. Here $$\theta = 0$$:

$$\begin{eqnarray} &&f[y]:=\\ &&C_1 \cdot U\left(\frac{\lambda }{2 \mu \left(\beta _2-1\right)},1+\frac{1}{2 \left(\beta _2-1\right)},\frac{y^2 \mu }{A^2 \sigma ^2 \left(\beta _2-1\right)}\right) + \\ && C_2 \cdot \,_1F_1\left(\frac{\lambda }{2 \mu \left(\beta _2-1\right)};1+\frac{1}{2 \left(\beta _2-1\right)};\frac{y^2 \mu }{A^2 \sigma ^2 \left(\beta _2-1\right)}\right) \tag{2a} \end{eqnarray}$$

1. Here $$\theta = 1$$:

$$\begin{eqnarray} &&f[y]:=\\ && C_1 \cdot e^{\frac{y \left(\mu -\sqrt{2 \lambda \sigma ^2+\mu ^2}\right)}{A \left(\beta _2-1\right) \sigma ^2}} U\left(\frac{\left(1-\frac{\mu }{\sqrt{\mu ^2+2 \lambda \sigma ^2}}\right) \beta _2}{2 \left(\beta _2-1\right)},\frac{\beta _2}{\beta _2-1},\frac{2 y \sqrt{\mu ^2+2 \lambda \sigma ^2}}{A \sigma ^2 \left(\beta _2-1\right)}\right) + \\ && C_2 e^{\frac{y \left(\mu -\sqrt{2 \lambda \sigma ^2+\mu ^2}\right)}{A \left(\beta _2-1\right) \sigma ^2}} \,_1F_1\left(\frac{\left(1-\frac{\mu }{\sqrt{\mu ^2+2 \lambda \sigma ^2}}\right) \beta _2}{2 \left(\beta _2-1\right)};\frac{\beta _2}{\beta _2-1};\frac{2 y \sqrt{\mu ^2+2 \lambda \sigma ^2}}{A \sigma ^2 \left(\beta _2-1\right)}\right) \tag{2b} \end{eqnarray}$$

1. Here $$\theta = 2$$:

$$\begin{eqnarray} &&f[y]:=\\ && C_1 \cdot y^{-\frac{\sigma ^2-2 \mu }{2 \left(\beta _2-1\right) \sigma ^2}} J_{\frac{\sigma ^2-2 \mu }{2 \sigma ^2 \left(\beta _2-1\right)}}\left(\frac{i \sqrt{2} y \sqrt{\lambda }}{A \sigma -A \sigma \beta _2}\right) + \\ && C_2 \cdot y^{-\frac{\sigma ^2-2 \mu }{2 \left(\beta _2-1\right) \sigma ^2}} Y_{\frac{\sigma ^2-2 \mu }{2 \sigma ^2 \left(\beta _2-1\right)}}\left(\frac{i \sqrt{2} y \sqrt{\lambda }}{A \sigma -A \sigma \beta _2}\right) \tag{2c} \end{eqnarray}$$

In here $$U$$ is the confluent hypergeometric function and $$J$$ and $$Y$$ are the Bessel functions of the first and the second kind. As always, the code below verifies numerically that the solutions are correct:

In[889]:= Clear[mu, s, lmb, b1, b2, y, z, f, th, A];
th = 0; f[
y_] := {HypergeometricU[lmb/(2 (-1 + b2) mu), 1 + 1/(2 (-1 + b2)), (
mu y^2)/(A^2 (-1 + b2) s^2)],
Hypergeometric1F1[lmb/(2 (-1 + b2) mu), 1 + 1/(2 (-1 + b2)), (
mu y^2)/(A^2 (-1 + b2) s^2)]};
{mu, s, b2, lmb, A, y} = RandomReal[{1, 5}, 6, WorkingPrecision -> 50];
ex

Clear[mu, s, lmb, b1, b2, y, z, f, th, A];
th = 1;
f[y_] := E^(((mu - Sqrt[mu^2 + 2 lmb s^2]) y)/(
A (-1 + b2) s^2)) {HypergeometricU[(
b2 (1 - mu/Sqrt[mu^2 + 2 lmb s^2]))/(2 (-1 + b2)),
b2/(-1 + b2), (2 Sqrt[mu^2 + 2 lmb s^2] y)/(A (-1 + b2) s^2)],
Hypergeometric1F1[(b2 (1 - mu/Sqrt[mu^2 + 2 lmb s^2]))/(
2 (-1 + b2)), b2/(-1 + b2), (2 Sqrt[mu^2 + 2 lmb s^2] y)/(
A (-1 + b2) s^2)]};
{mu, s, b2, lmb, A, y} = RandomReal[{1, 5}, 6, WorkingPrecision -> 50];
ex

Clear[mu, s, lmb, b1, b2, y, z, f, th, A];
th = 2;
f[y_] := y^(-((-2 mu + s^2)/(
2 (-1 + b2) s^2))) {BesselJ[(-2 mu + s^2)/(2 (-1 + b2) s^2), (
I Sqrt[2] Sqrt[lmb] y)/(A s - A b2 s)],
BesselY[(-2 mu + s^2)/(2 (-1 + b2) s^2), (I Sqrt[2] Sqrt[lmb] y)/(
A s - A b2 s)]};
{mu, s, b2, lmb, A, y} = RandomReal[{1, 5}, 6, WorkingPrecision -> 50];
ex

Out[892]= {0.*10^-48 + 0.*10^-49 I, 0.*10^-46}

Out[897]= {0.*10^-46 + 0.*10^-46 I, 0.*10^-48}

Out[902]= {0.*10^-47 + 0.*10^-47 I, 0.*10^-46 + 0.*10^-47 I}

1. Here $$\theta=-2$$.

We demonstrate below that the eigenfunctions are expressed through the bi-confluent Heun function, in the independent variable $$y^2$$, with following parameters $$-\gamma = (1-2 \beta_2)/(2-2 \beta_2), \delta = 0, -\epsilon=-1$$ and $$\alpha = 0, -q= - \lambda \sigma^2/2 \cdot 1/\sqrt{\mu} ((\beta_2-1) \sigma^2)^{-3/2}$$ (we used the notation from https://dlmf.nist.gov/31.12). Here we go:

In[312]:= (*First case.*)
Clear[mu, s, lmb, b1, b2, y, z, xi, f, th,
A, \[Mu], \[Sigma], \[Beta], \[Lambda], \[Theta]];
th = -2;
ex = -lmb f[y] +
1/2 A^(-(2/(-1 + b2))) (-1 + b2) y^(-1 -
th) (-2 A^(2/(-1 + b2) + th) mu y^2 +
A^((2 b2)/(-1 + b2)) b2 s^2 y^th) Derivative[1][f][y] +
1/2 A^2 (-1 + b2)^2 s^2 (f^\[Prime]\[Prime])[y];

ex1 = ex /. Derivative[1][f][y] :> 2 Sqrt[v] D[f[v], v] /.
Derivative[2][f][y] :> 2 Sqrt[v] D[2 Sqrt[v] D[f[v], v], v] /.
f[y] :> f[v] /. y :> Sqrt[v];
ex1 = ex1/Coefficient[ex1, Derivative[2][f][v]];
(*Maps onto the bi-confluent Heun equation with -\[Gamma] = (1-2 \
b2)/(2 -2 b2 ), \[Delta]=0, -\[Epsilon]=-1, and \[Alpha] = 0, -q = \
-((lmb (s^2) )/(2 Sqrt[mu] ((-1+b2) s^2)^(3/2))) see \
https://dlmf.nist.gov/31.12.*)
Collect[ex1 /. A -> mu^(1/4)/((-1 + b2) s^2)^(1/4),
Derivative[n_][f][v], FullSimplify]

Out[317]= -((lmb s^2 f[v])/(
2 Sqrt[mu] ((-1 + b2) s^2)^(3/2)
v)) + (-v + (1 - 2 b2)/(2 v - 2 b2 v)) Derivative[1][f][v] + (
f^\[Prime]\[Prime])[v]

##### Second case:

Here we assume that $$(1-\beta_2) = \theta \cdot (1-\beta_1)$$ where $$\theta \in {\mathbb Z}$$ and we substitute for $$y = A z^{1-\beta_1}$$ in $$(1)$$. This leads to a following ODE below:

$$$$\frac{1}{2} \left(\beta _1-1\right){}^2 \sigma ^2 A^{2 \theta } y^{2-2 \theta } f''(y)+\frac{1}{2} \left(\beta _1-1\right) y^{-2 \theta } f'(y) \left(\beta _1 \sigma ^2 y A^{2 \theta }-2 A \mu y^{2 \theta }\right)-\lambda f(y)=0 \tag{3}$$$$

Again, see code snippet for the derivation of $$(3)$$:

In[922]:= (*Second case: (1-b2) = th(1-b1)  with th in integers.*)
Clear[mu, s, lmb, b1, b2, y, z, f, th, A];
zz[y_] = z /. First@Solve[y == A z^(1 - b1), z];
u[y_] := D[zz[y], y];
ex = Nest[(mu z^b1 D[#, z] + s^2/2 z^(2 b2) D[#, {z, 2}] - lmb #) &,
f[z], 1];
ex = ex /. Derivative[1][f][z] :> 1/u[y] D[f[y], y] /.
Derivative[2][f][z] :> 1/u[y] D[ 1/u[y] D[f[y], y], y] /.
f[z] :> f[y] /. z :>  zz[y];
ex = PowerExpand[ex] /. b2 :> th (b1 - 1) + 1;
ex = Collect[ex, Derivative[n_][f][y], Simplify[#] &];
ex = Expand[ex] /. y^n_ :> y^Simplify[n];
(*The later is then exactly solved in case th=0.*)
ex = Collect[ex, Derivative[n_][f][y], Simplify[PowerExpand[#]] &]

During evaluation of In[922]:= Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information.

Out[930]= -lmb f[y] +
1/2 (-1 + b1) y^(-2 th) (A^(2 th) b1 s^2 y -
2 A mu y^(2 th)) Derivative[1][f][y] +
1/2 A^(2 th) (-1 + b1)^2 s^2 y^(2 - 2 th) (f^\[Prime]\[Prime])[y]


Now, again, with the help of Mathematica we have found only one exactly solvable case which is different from the ones listed above.

1. Here $$\theta = 0$$.

$$\begin{eqnarray} &&f[y]:=\\ && C_1 \cdot y^{\frac{\sigma -\sqrt{8 \lambda +\sigma ^2}}{2 \sigma -2 \beta _1 \sigma }} \, _1F_1\left(-\frac{\sigma -\sqrt{\sigma ^2+8 \lambda }}{2 \sigma -2 \sigma \beta _1};\frac{\sqrt{\sigma ^2+8 \lambda }}{\sigma -\sigma \beta _1}+1;-\frac{2 A \mu }{y \sigma ^2 \left(\beta _1-1\right)}\right) + \\ && C_2 \cdot y^{\frac{\sqrt{8 \lambda +\sigma ^2}+\sigma }{2 \sigma -2 \beta _1 \sigma }} \, _1F_1\left(-\frac{\sigma +\sqrt{\sigma ^2+8 \lambda }}{2 \sigma -2 \sigma \beta _1};\frac{\sqrt{\sigma ^2+8 \lambda }}{\sigma \left(\beta _1-1\right)}+1;-\frac{2 A \mu }{y \sigma ^2 \left(\beta _1-1\right)}\right) \tag{3a} \end{eqnarray}$$

Again, the code below verifies the solutions numerically. Here we go:

In[944]:= Clear[mu, s, lmb, b1, b2, y, z, f, th, A];
th = 0; f[
y_] := {y^((s - Sqrt[8 lmb + s^2])/(2 s - 2 b1 s))
Hypergeometric1F1[-((s - Sqrt[8 lmb + s^2])/(2 s - 2 b1 s)),
1 + Sqrt[8 lmb + s^2]/(s - b1 s), -((2 A mu)/((-1 + b1) s^2 y))],
y^((s + Sqrt[8 lmb + s^2])/(2 s - 2 b1 s))
Hypergeometric1F1[-((s + Sqrt[8 lmb + s^2])/(2 s - 2 b1 s)),
1 + Sqrt[8 lmb + s^2]/((-1 + b1) s), -((
2 A mu)/((-1 + b1) s^2 y))]};
{mu, s, b1, lmb, A, y} = RandomReal[{1, 5}, 6, WorkingPrecision -> 50];
ex

Out[947]= {0.*10^-47, 0.*10^-49}

##### Overall conclusion:

It is definitely true that solutions exist in both case for an arbitrary integer value of $$\theta$$. However, most likely, those solutions are not expressed through the standard functions being available in the symbolic calculation software. As a matter of fact some of those, unknown, solutions might be expressed through the Heun functions or the confluent versions of those. This is an open question,which needs to be investigated, when this is the case . Unfortunately, according to my knowledge, Mathematica does not know how to operate with Heun functions (but Maple does) and as such it is very hard to verify any possible results. Again, it would be interesting to analyze this problem deeper at some point in the future.

We provide an additional solution, i.e. the eigenfunctions of the infinitesimal generator, when $$\beta_1 - 1 = \theta (\beta_2-1)$$ and $$\theta=2 {\tilde \theta}$$ where $$-{\tilde \theta} \in {\mathbb N}_+$$. This is a generalization of the first case in my first answer above.

First of all a substitution $$y = A z^{1-\beta_2}$$ transforms the ODE $$(1)$$ into the following, simpler, ODE with polynomial coefficients. Here we go:

$$$$f^{''}(y) + \left( {\mathfrak D} \frac{1}{y} + y^{1-\theta} \right) f^{'}(y) - {\mathfrak C} \lambda f(y) = 0 \tag{I}$$$$

where $$A:= 2^{-\frac{1}{\theta -2}} \left(\frac{\left(1-\beta _2\right) \sigma ^2}{\mu }\right){}^{\frac{1}{\theta -2}}$$ then $${\mathfrak D}:= \beta_2/(-1+\beta_2)$$ and $${\mathfrak C}:= \frac{2}{A^2 \left(\beta _2-1\right){}^2 \sigma ^2}$$.

Now we use the Frobenius series expansion method to construct two linearly independent solutions of the ODE $$(I)$$. In other words we postulate $$f(y) = \sum\limits_{n=0}^\infty a_n y^{n+\alpha}$$. Inserting this into the ODE $$(I)$$ we get (from the left to the right) the following equation:

$$\begin{eqnarray} \underline{ a_n \cdot (n+\alpha)_{(2)} \cdot y^{n+\alpha-2}+ a_n \cdot {\mathfrak D} (n+\alpha)_{(1)} \cdot y^{n+\alpha-2} }+ a_n \cdot (n+\alpha)_{(1)} y^{n+\alpha-\theta} + -{\mathfrak C} \cdot \lambda a_n y^{n+\alpha} =0 \tag{Ia} \end{eqnarray}$$ It is clear that out of all those different powers of $$y$$ it is the underlined terms that produce the lowest power exponents. By annihilating those terms, at the lowest value of $$n$$, we obtain the indicial equation for $$\alpha$$. We have two possibilities.

1. $$a_0 \neq 0$$ and $$a_1=0$$ and $$\alpha_{(2)} + {\mathfrak D} \alpha_{(1)}=0$$. This gives $$\alpha = 0$$. The recurrence reads: $$$$a_{n+2} = \frac{{\mathfrak C} \lambda a_n - (n+\theta) a_{n+\theta}}{(n+2)(n+1 + {\mathfrak D})}$$$$

2. $$a_0 = 0$$ and $$a_1 \neq 0$$ and $$(1+\alpha)_{(2)} + {\mathfrak D} (1+\alpha)_{(1)}=0$$. This gives $$\alpha = -1$$. The recurrence reads: $$$$a_{n+2} = \frac{{\mathfrak C} \lambda a_n - (n-1+\theta) a_{n+\theta}}{(n+1)(n+0 + {\mathfrak D})}$$$$ In both cases we have a non-trivial recurrence that involves three different indices $$a_{n+2} = f_n \cdot a_n + g_n \cdot a_{n+\theta}$$ subject to either $$a_0,a_1=1,0$$ or $$a_0,a_1=0,1$$. Now, it is not that hard to see that the general solution to this recurrence can be given in terms of of multiple sum over strictly increasing integer sequences as visualized in the following picture below:

$$\begin{eqnarray} &2&, \cdots, 2, &2-\theta&, 2, \cdots, 2, &2-\theta&, 2, \cdots, 2, &2-\theta&, 2 \cdots, &2 \\ &1& \le &j_1 & < &j_2 & < \cdots < &j_l & \le &k \end{eqnarray}$$

subject to $$n+2- \zeta = 2 k- l \theta$$ (the difference between the current index and the initial index has to be equal to the $$L_1$$ norm of the sequence of the steps above. Therefore the solution to the recurrence in question reads:

$$\begin{eqnarray} a_{n+2} = a_\zeta \cdot \sum\limits_{k=1}^{\lfloor \frac{n+2-\zeta}{2} \rfloor} \sum\limits_{l=0}^k \delta_{n+2-\zeta,2 k- l \theta} \sum\limits_{1 \le j_1 < j_2 < \cdots < j_l \le k} \prod\limits_{\eta=0}^l \left( \prod\limits_{\xi=0}^{j_{\eta+1}-j_\eta-2} f_{n - 2 j_\eta + \eta \cdot \theta-2 \xi} \right) \left( g_{n - 2 j_{\eta+1} + \eta \cdot \theta+2} 1_{\eta < l} + 1_{\eta=l} \right) \tag{Ib} \end{eqnarray}$$ where $$\zeta = 0,1$$. Here $$\eta=0,\cdots, l$$ labels the positions of the $$g$$ terms (the jumps by $$2-\theta$$) whereas $$\xi=0,\cdots,j_{\eta+1}-j_\eta-2$$ labels the positions labels the positions of the $$f$$ terms (the jumps by $$2$$) preceding a particular $$g$$ term. One sanity check on this result is the following. For a fixed $$\eta$$ the subscript always jumps down by a factor of two including the jump from the last $$f$$ to the subsequent $$g$$. However the the jumps in the index from the $$g$$ factor to the first $$f$$ factor after it is always equal to $$2-\theta$$.

Now we replace $$\theta \rightarrow 2 {\tilde \theta}$$ and then we can take out a common factor from the term in the multiple sums below:

$$\begin{eqnarray} a_{n+2} = a_\zeta \cdot \sum\limits_{k=1}^{\lfloor \frac{n+2-\zeta}{2} \rfloor} \sum\limits_{l=0}^k \delta_{n+2-\zeta,2 k- l 2 {\tilde \theta}} \left( \prod\limits_{\xi=0}^{k-1-l {\tilde \theta}} f_{n-2 \xi} \right) \sum\limits_{1 \le j_1 < j_2 < \cdots < j_l \le k} \prod\limits_{\eta=0}^{l-1} \frac{g_{n-2 j_{\eta+1} +\eta 2 {\tilde \theta} + 2}} {\prod\limits_{\xi=1}^{1-\theta} f_{n-2 j_{\eta+1} +\eta 2 {\tilde \theta} + 4-2\xi}} \tag{Ic} \end{eqnarray}$$

Now we insert $$(f_n,g_n) = ({\mathfrak C} \lambda, -(n-\zeta+2 {\tilde \theta}))/((n-\zeta+2)(n-\zeta+1+{\mathfrak D}) )$$ and then the solution reads:

$$\begin{eqnarray} a_{n+2}^{(\zeta)} &=&a_\zeta \cdot\sum\limits_{k=1}^{\lfloor \frac{n+2-\zeta}{2} \rfloor} \sum\limits_{l=0}^k \delta_{n+2-\zeta, 2k-2 l {\tilde \theta}} \cdot ({\mathfrak C} \lambda)^{k-l} \cdot 2^{(-2 {\tilde \theta}+1) l} \cdot \prod\limits_{\xi=0}^{k-1-l {\tilde \theta}} \frac{1}{(n-\zeta-2 \xi+2)(n-\zeta-2 \xi+1+{\mathfrak D})}\\ && \underline{ \sum\limits_{1 \le j_1 < j_2 < \cdots < j_l \le k} \prod\limits_{\eta=0}^{l-1} \left( j_{\eta+1} - \eta {\tilde \theta} + \frac{\zeta-n-2}{2}\right)^{(-{\tilde \theta}+1)} \cdot \left( j_{\eta+1} - \eta {\tilde \theta} + \frac{\zeta-n-1-{\mathfrak D}}{2}\right)^{(-{\tilde \theta})} } \\ \tag{II} % \end{eqnarray}$$ The amazing thing is that the multiple sum above can be carried out in closed form by using the results from here.

Now we we took $$-{\tilde \theta}= 0,1,\cdots,5$$ and for each of those values we construct the two independent solutions , i.e. $$\zeta=0,1$$, and then verify that both of them satisfy the ODE $$(I)$$. The output is composed of sequences $${\tilde \theta}, \zeta$$ then the difference between the coefficients being computed numerically (by solving the recurrence relation) and "analytically" (by using $$(II)$$) and finally the residual, i.e. the result of the action of the differential operator in $$(I)$$ on the series expansion where all the terms bigger then a prescribed threshold are being neglected. Here we go:

(*Now replace \[Theta] by 2 \[Theta] in the above. We have:*)
Clear[a, aa, \[Theta], \[Lambda], \[Eta], \[Theta], DD, CC, th, f, g, \
j]; M = 20;

Do[(*Loop over \[Theta] *)
Do[a[-p] = 0, {p, 1, -2 \[Theta]}];
Do[(*Loop over \[Zeta]*)
If[\[Zeta] == 1, {a[0], a[1]} = {0, 1}, {a[0], a[1]} = {1, 0}];

f[n_] := (CC \[Lambda] )/((n - \[Zeta] + 2) (n - \[Zeta] + 1 + DD));
g[n_] := - (n - \[Zeta] + 2 \[Theta])/((n - \[Zeta] +
2) (n - \[Zeta] + 1 + DD));

(*Numerically*)
Do[
a[n + 2] = f[n] a[n] + g[n] a[n + 2 \[Theta]];
, {n, 0, M}];

(*"Closed form solution"*)
Do[
S = 0;
Do[
Do[
If[! (n + 2 - \[Zeta] - (2 k - l 2 \[Theta]) == 0), Continue[]];
S0 = Which[l == 0, 1, l > 0,
Sum[

Product[Pochhammer[
j[\[Eta] + 1] - \[Eta] \[Theta] + (\[Zeta] - n - 2)/
2, -\[Theta] + 1] Pochhammer[
j[\[Eta] + 1] - \[Eta] \[Theta] + (\[Zeta] - n - 1 - DD)/
2, -\[Theta]], {\[Eta], 0, l - 1}]
,
Evaluate[
Sequence @@
Table[{j[\[Eta]], If[\[Eta] == 1, 0, j[\[Eta] - 1]] + 1,
k}, {\[Eta], 1, l}]]]
];
S0 =
S0 *Product[
1/((n - \[Zeta] - 2 \[Xi] + 2) (n - \[Zeta] - 2 \[Xi] + 1 +
DD)), {\[Xi], 0, k - 1 - l \[Theta]}] (CC \[Lambda])^(
k - l) 2^((-2 \[Theta] + 1) l);
S += S0;
, {l, 0, k}];
, {k, 1, Floor[(n + 2 - \[Zeta])/2]}];
aa[n + 2] = S;
, {n, 0, M}];
(*Use sequence above to construct a general solution to the ODE in \
question.*)
res = (D[#, {y, 2}] + ( DD 1/y + y^(1 - 2 \[Theta])) D[#, y] -
CC \[Lambda] #) & /@ Sum[a[n] y^(n - \[Zeta]), {n, 0, M + 2}];
res = Simplify[(Expand[res] /. y^n_ :> 0 /; n >= M)];
Print["\[Theta],\[Zeta],a[n]-aa[n],residual=", {\[Theta], \[Zeta],
Simplify[Table[(a[n + 2] - aa[n + 2]), {n, 0, M}]], res}];
, {\[Zeta], 0, 1}];
, {\[Theta], -0, -5, -1}];


Here we provide a new set of solutions, i.e. the local solutions to the eigen-equation of the infinitesimal generator, in the case when $$\beta_1 = 1 + \theta (\beta_1-1)$$ with $$\theta=1,2,3,\cdots$$. This is a generalization of the second case of my answer above.

A substitution $$y = A z^{1-\beta_1}$$ transforms the ODE $$(1)$$ into the following simpler ODE with integer power law coefficients. We have:

$$$$f^{''}(y) + \left({\mathfrak D} \frac{1}{y} + y^{2 \theta-2} \right) f^{'}(y) - {\mathfrak C} \lambda y^{2 \theta-2} f(y) \tag{II}$$$$

where $$A= (\frac{(1-\beta_1) \sigma^2}{(2 \mu)})^{\frac{1}{1-2 \theta}}$$ then $${\mathfrak D} = \beta_1/(-1+\beta_1)$$ and $${\mathfrak C} := \frac{2}{(A^{2\theta}(-1+\beta_1)^2 \sigma^2)}$$.

Now we proceed in the same way as before, i.e. we use the Frobenius method to construct power series expansion solutions $$f(y):= \sum\limits_{n=0}^\infty a_n y^{n+\alpha}$$ . Inserting this into $$(II)$$ gives(from the left to the right) the following:

$$$$\underline{a_n (n+\alpha)_{(2)} y^{n+\alpha-2} + a_n {\mathfrak D} (n+\alpha)_{(1)} y^{n+\alpha-2}} + a_n(n+\alpha)_{(1)} y^{n+\alpha+2 \theta-3} - {\mathfrak C} \lambda a_n y^{n+\alpha+2 \theta-2} =0 \tag{IIa}$$$$

Because $$\theta \ge 1$$ it is clear that the underlined terms will produce the lowest power exponents. Annihilating those terms at the lowest values of $$n$$ gives us the indicial equation. Pretty much in the same way as before we have the following two possibilities:

1. $$a_0\neq 0$$ and $$a_1=0$$ and $$\alpha_{(2)} + {\mathfrak D} \alpha_{(1)} = 0$$ which gives $$\alpha =0$$.

2. $$a_0= 0$$ and $$a_1 \neq0$$ and $$(1+\alpha)_{(2)} + {\mathfrak D} (1+\alpha)_{(1)} = 0$$ which gives $$\alpha =-1$$.

The recurrence relation reads, in both cases, as follows:

$$$$a_{n+2} = \frac{{\mathfrak C} \lambda a_{n-\theta_1} - (n+\alpha -\theta_2) a_{n-\theta_2} }{(n+2+\alpha)(n+1+\alpha +{\mathfrak D})}$$$$ where $$(\theta_1,\theta_2) = (2\theta-2, 2\theta-3)$$. The recurrence involves, again, three different points and is of the form $$a_{n+2}= f_n a_{n-\theta_1} + g_n a_{n-\theta_2}$$ subject to either $$(a_0,a_1) = (1,0)$$ or $$(a_0,a_1) = (0,1)$$. By similar considerations as in my previous answer above the general solution to that recurrence reads:

$$\begin{eqnarray} a_{n+2} = a_\zeta \sum\limits_{k=1}^{\lfloor \frac{n+2-\zeta}{2+ \theta_1 \wedge \theta_2} \rfloor} \sum\limits_{l=0}^k \delta_{n+2-\zeta,(2+\theta_1)k + l (\theta_2-\theta_1)} % \sum\limits_{1 \le j_1 < j_2 < \cdots < j_l \le k} \prod\limits_{\eta=0}^l \left( g_{n-(2+\theta_1)(j_{\eta+1}-\eta-1) - \eta (2+\theta_2)} 1_{\eta where $$\zeta=0,1$$.

Now we take $$(\theta_1,\theta_2) = (2\theta-2,2\theta-3)$$ and $$(f_n,g_n) = ({\mathfrak C} \lambda,-(n-\zeta-\theta_2))/((n-\zeta+2)(n-\zeta+1+{\mathfrak D}))$$ and inserting those into $$(IIb)$$ we get the final formula for the coefficients. Here we go:

$$\begin{eqnarray} &&a_{n+2} = a_\zeta \cdot \sum\limits_{k=1}^{\lfloor \frac{n+2-\zeta}{2 \theta-1} \rfloor} \sum\limits_{l=0}^k \delta_{k, \frac{n+l+2-\zeta}{2 \theta}} \cdot \left( \frac{{\mathfrak C} \lambda}{(2\theta)^2} \right)^{k-l} \left( \frac{1}{2 \theta} \right)^l \cdot \sum\limits_{1 \le j_1 < j_2 , \cdots < j_l \le k} \\ && \prod\limits_{\eta=0}^l \frac{1}{ (\frac{n+\eta-\zeta+2}{-2\theta} + j_\eta)^{(j_{\eta+1} - j_\eta - \delta_{eta,l} )}} \cdot \frac{1}{ (\frac{n+\eta-\zeta+1+{\mathfrak D}}{-2\theta} + j_\eta)^{(j_{\eta+1} - j_\eta - \delta_{eta,l} )}} \cdot \prod\limits_{\eta=1}^{l} \left( \frac{n+\eta-\zeta+2}{-2 \theta} + j_\eta \right) \end{eqnarray}$$ for $$n=0,1,2,\cdots$$. Now the solution to the ODE $$(II)$$ reads $$f(y) = \sum\limits_{n=0}^\infty a_n y^{n-\zeta}$$ as the code snippet below demonstrates:

(*Eigenfunctions of the infinitesimal generator: Second case.*)
(*The recurrence to be solved is the one above at \[Theta]1= 2 \
\[Theta]-2 and \[Theta]2 = 2 \[Theta]-3.*)
up[a_, n_] := Pochhammer[a, n];
\[Theta] = RandomInteger[{2, 5}];
\[Zeta] = RandomInteger[{0, 1}]; M = 30;

Do[
Do[
{DD, CC, \[Lambda]} =
RandomReal[{0, 1}, 3, WorkingPrecision -> 50];
{\[Theta]1, \[Theta]2} = {2 \[Theta] - 2, 2 \[Theta] - 3};
Do[a[-p] = 0, {p, 1, 2 \[Theta] - 2}];
If[\[Zeta] == 1, {a[0], a[1]} = {0, 1}, {a[0], a[1]} = {1, 0}];

(*Use sequence above to construct a general solution to the ODE \
in question.*)
f[n_] := (
CC \[Lambda] )/((n - \[Zeta] + 2) (n - \[Zeta] + 1 + DD));
g[n_] := - (n - \[Zeta] - \[Theta]2)/((n - \[Zeta] +
2) (n - \[Zeta] + 1 + DD));

(*Numerically*)
Do[
a[n + 2] = f[n] a[n - \[Theta]1] + g[n] a[n - \[Theta]2];
, {n, 0, M}];
(*"Closed form solution"
*)
Do[
S = 0;
Do[
Do[
If[! ((n + l + 2 - \[Zeta])/(2 \[Theta]) == k), Continue[]];
S0 = Which[
l == 0, (-1)^k/(k! up[(DD - 1)/(-2 \[Theta]) - k, k]),
l > 0,
Sum[
Product[

1/up[(n + \[Eta] - \[Zeta] +
2)/(-2 \[Theta]) + (If[\[Eta] == 0, 0, j[\[Eta]]]),
If[\[Eta] == l, k, j[\[Eta] + 1]] -
If[\[Eta] == 0, 0, j[\[Eta]]]] 1/
up[(n + \[Eta] - \[Zeta] + 1 +
DD)/(-2 \[Theta]) + (If[\[Eta] == 0, 0, j[\[Eta]]]),
If[\[Eta] == l, k, j[\[Eta] + 1]] -
If[\[Eta] == 0, 0, j[\[Eta]]]] , {\[Eta], 0,

l}] Product[((
n + \[Eta] - \[Zeta] +
2)/(-2 \[Theta]) + (j[\[Eta]])), {\[Eta], 1, l}]
,
Evaluate[
Sequence @@
Table[{j[\[Eta]], If[\[Eta] == 1, 0, j[\[Eta] - 1]] + 1,
k}, {\[Eta], 1, l}]]]
];
S +=
S0 ((CC \[Lambda])/(2 \[Theta])^2)^(k - l) (1/(2 \[Theta]))^l;
, {l, 0, k}];
, {k, 1, Floor[(n + 2 - \[Zeta])/(2 \[Theta] - 1)]}];
aa[n + 2] = S;
, {n, 0, M}];

diffas = N[Table[(a[n + 2] - aa[n + 2]), {n, 0, M}], 50];

(*Use sequence above to construct a general solution to the ODE \
in question.*)

Clear[CC, DD, \[Lambda]];
f[n_] := (
CC \[Lambda] )/((n - \[Zeta] + 2) (n - \[Zeta] + 1 + DD));
g[n_] := - (n - \[Zeta] - \[Theta]2)/((n - \[Zeta] +
2) (n - \[Zeta] + 1 + DD));

(*Numerically*)
Do[
a[n + 2] = f[n] a[n - \[Theta]1] + g[n] a[n - \[Theta]2];
, {n, 0, M}];

res = (D[#, {y, 2}] + ( DD 1/y + y^\[Theta]1) D[#, y] -
CC \[Lambda] y^\[Theta]1 #) & /@ {Sum[
a[n] y^(n - \[Zeta]), {n, 0, M + 2}]};
res =
Collect[(Expand[res] /. y^n_ :> 0 /; n >= M), y, Together[#] &];
Print[
"\[Theta],\[Zeta],a-aa,residual=", {\[Theta], \[Zeta], diffas,
res}];
, {\[Zeta], 0, 1}];
, {\[Theta], 2, 5, +1}];


We give the eigenfunctions of the infinitesimal generator for another much broader family of processes that includes the first case family as a special case. To be precise we take: $$\beta_1-1 = \theta_1/\theta_2 (\beta_2-1)$$ where $$\theta_{1,2}$$ are both integers and $$\theta_1 \ge 1$$ and $$2 \theta_2-\theta_1-2 \ge 0$$.

Here the substitution $$y = A z^{\frac{1-\beta_2}{\theta_2}}$$ leads to the following ODE with integer power law coefficients. We have:

$$$$f^{''}(y) + \left( {\mathfrak D} \frac{1}{y} + y^{2 \theta_2-\theta_1-1}\right) f^{'}(y) - {\mathfrak C} \cdot \lambda f(y) = 0 \tag{III}$$$$

where $$A = (\frac{\sigma^2 (1-\beta_2)}{2 \mu \theta_2})^{1/(\theta_1-2 \theta_2)}$$ then $${\mathfrak D} = 1 + \frac{\theta_2}{-1+\beta_2}$$ and $${\mathfrak C} := A^{-2 \theta_2} \theta_2^2 \frac{2}{(-1+\beta_2)^2 \sigma^2}$$.

Now we set $$f(y) := \sum\limits_{n=0}^\infty a_n y^{n+\alpha}$$. Inserting this conjecture into $$(III)$$ gives the following:

$$$$\underline{ a_n (n+\alpha)_{(2)} \cdot y^{n+\alpha-2} + a_n {\mathfrak D} (n+\alpha)_{(1)} \cdot y^{n+\alpha-2}} + a_n (n+\alpha)^{(1)} y^{n+\alpha-\theta_1+2(\theta_2-1)} + -{\mathfrak C} \cdot \lambda \cdot a_n \cdot y^{n+\alpha} \tag{IIIa}$$$$

Now, due to our constraints on $$\theta_{1},\theta_2$$ it is the underlined terms that yield the smallest power law exponents at $$n=0,1$$. The indicial equation is therefore determined by annihilating those terms. We have two possibilities:

1. $$a_0\neq 0$$, $$a_1=0$$ and $$\alpha_{(2)} + {\mathfrak D} \alpha_{(1)} =0$$ which gives $$\alpha=0$$.

2. $$a_0= 0$$, $$a_1\neq 0$$ and $$(1+\alpha)_{(2)} + {\mathfrak D} (1+\alpha)_{(1)} =0$$ which gives $$\alpha=-{\mathfrak D}$$.

The recurrence relation reads in both cases as follows:

$$\begin{eqnarray} a_{n+2} &=& \frac{{\mathfrak C} \lambda a_n - a_{n+2-(2 \theta_2-\theta_1)}}{(n+2+\alpha)(n+1+\alpha+{\mathfrak D})} \tag{IIIb} \end{eqnarray}$$ where $$n=0,1,2,\cdots$$ and the initial conditions being either $$a_0,a_1,\alpha=1,0,0$$ or $$a_0,a_1,\alpha=0,1,-{\mathfrak D}$$.

Now by using the general solution to a three point recurrence relation (see $$(IIb)$$ in my second answer to this question) we write down the following "closed form " expression for the coefficients:

$$\begin{eqnarray} &&a_{n+2} = a_\zeta \sum\limits_{k=1}^{\lfloor \frac{n+2-\zeta}{2+(0 \wedge (-2+2 \theta_2-\theta_1))}\rfloor} \sum\limits_{l=0}^k \delta_{n+2-\zeta,2 k+l(-2+2 \theta_2-\theta_1)} \cdot (\frac{{\mathfrak C} \lambda}{2^2})^{k-l} 2^{-l} \cdot \left( (-1-\frac{\theta_1}{2}+\theta_2)! (-0-\frac{\theta_1}{2}+\theta_2)! \right)^l \cdot \\ && \frac{\Gamma(\frac{n+\epsilon+\alpha+2}{-2})}{\Gamma(\frac{n-l(-2+2\theta_2-\theta_1)+\epsilon+\alpha+2}{-2} + k)} \cdot \frac{\Gamma(\frac{n+\epsilon+\alpha+1+{\mathfrak D}}{-2})}{\Gamma(\frac{n-l(-2+2\theta_2-\theta_1)+\epsilon+\alpha+1+{\mathfrak D}}{-2} + k)} \cdot \\ && \sum\limits_{1 \le j_1 < j_2 < \cdots < j_l \le k} \prod\limits_{\eta=1}^l \binom{\frac{1}{2}(-2-n-\alpha-\epsilon-\eta(2+\theta_1-2 \theta_2)+2 j_\eta)}{-0-\frac{\theta_1}{2}+\theta_2} \binom{\frac{1}{2}(-3-{\mathfrak D}-n-\alpha-\epsilon-\eta(2+\theta_1-2 \theta_2)+2 j_\eta)}{-1-\frac{\theta_1}{2}+\theta_2} \tag{IIIc} \end{eqnarray}$$ where $$\zeta=0,1$$ and $$\alpha = 0$$ or $$\alpha=-{\mathfrak D}$$ when $$\zeta=0$$ or $$\zeta=1$$ respectively. Again, when $$\theta_1$$ is even then the multiple sum can be done in closed form by using the results from here. The eigenfunctions are then constructed from the ansatz as follows $$f(y)= \sum\limits_{n=0}^\infty a_n y^{n+\alpha}$$.

The code snippet below verifies the solution in question. Here we go:

(*Eigenfunctions of the infinitesimal generator: Third case.\[Beta]1 \
= 1+\[Theta]1/\[Theta]2 (\[Beta]2-1);*)
(* We have:*)
Clear[a, aa, \[Theta]1, \[Theta]2, \[Lambda], \[Epsilon], DD, CC, f, \
g, z, y]; M = 10; mmax = 5;
(*Assumption: 2-2 \[Theta]2+\[Theta]1 <=0 and Element[\[Theta]2|\
\[Theta]1, Integers]*)
up[a_, n_] := Pochhammer[a, n];
Do[Do[
Do[
Do[a[-p] = 0, {p, 1, 2 Max[\[Theta]1, \[Theta]2]}];
If[\[Zeta] ==
0, {a[0], a[1], \[Alpha]} = {1, 0, 0}, {a[0],
a[1], \[Alpha]} = {0, 1, -DD}];

f[n_] := (
CC  \[Lambda] )/((n + \[Epsilon] + \[Alpha] +
2) (n + \[Epsilon] + \[Alpha] + 1 + DD));
g[n_] := - (n + \[Epsilon] + \[Alpha] + 2 -
2 \[Theta]2 + \[Theta]1)/((n + \[Epsilon] + \[Alpha] +
2) (n + \[Epsilon] + \[Alpha] + 1 + DD));

(*Numerically*)
Do[
a[n + 2] = f[n] a[n] + g[n] a[n + 2 - 2 \[Theta]2 + \[Theta]1];
, {n, 0, M}];

(*"Closed form solution"*)
If[False,
(*Insert (*Some old code.*) from above if needed.*)
,
(*Here \[Theta]1 is odd.*)
Do[
S = 0;
Do[
Do[

If[! (n +
2 - \[Zeta] - ((2 + 0) k +
l (-2 + 2 \[Theta]2 - \[Theta]1 - 0)) == 0),
Continue[]];
S0 = Which[
l == 0, 1,
l > 0,
Sum[

Product[Binomial[
1/2  (-2 -
n - \[Alpha] - \[Epsilon] - \[Eta]  (2 + \
\[Theta]1 - 2  \[Theta]2) + 2  j[\[Eta]]), -(\[Theta]1/
2) + \[Theta]2]  Binomial[
1/2  (-3 - DD -
n - \[Alpha] - \[Epsilon] - \[Eta]  (2 + \
\[Theta]1 - 2  \[Theta]2) + 2  j[\[Eta]]), -1 - \[Theta]1/
2 + \[Theta]2] , {\[Eta], 1, l}]
,
Evaluate[
Sequence @@
Table[{j[\[Eta]],
If[\[Eta] == 1, 0, j[\[Eta] - 1]] + 1, k}, {\[Eta],
1, l}]]]
];

S += S0 ((CC \[Lambda])/2^2 )^(
k - l) (2)^-l Gamma[(n + \[Epsilon] + \[Alpha] + 2)/-2]/
Gamma[(n - (l) (-2 +

2 \[Theta]2 - \[Theta]1) + \[Epsilon] + \[Alpha] \
+ 2)/-2 + k] (
Gamma[(n + \[Epsilon] + \[Alpha] + 1 +
DD)/-2] ((-1 - \[Theta]1/
2 + \[Theta]2)! (-0 - \[Theta]1/2 + \[Theta]2)!)^
l)/Gamma[(
n + \[Epsilon] - (l) (-2 +
2 \[Theta]2 - \[Theta]1) + \[Alpha] + 1 + DD)/-2 +
k];
, {l, 0, k}];
, {k, 1,
Floor[(n + 2 - \[Zeta])/(
2 + Min[0, -2 + 2 \[Theta]2 - \[Theta]1])]}];
aa[n + 2] = S;
, {n, 0, M}];
];

amaa = FullSimplify[Table[(a[n + 2] - aa[n + 2]), {n, 0, M}]];

(*Use sequence above to construct a general solution to the ODE \
in question.*)
res = (D[#, {y, 2}] + (
DD 1/y + y^(2 \[Theta]2 - \[Theta]1 - 1)) D[#, y] -
CC \[Lambda] #) & /@ {Sum[
a[n] y^(n + \[Alpha]), {n, 0, M + 2}]};
res = res /. \[Epsilon] :> 0;
res1 =
Collect[(Expand[res] /. y^n_ :> 0 /; n + DD >= M || n >= M), y,
Together];
Print[
"\[Theta]1,\[Theta]2,\[Zeta],a[]-aa[],residual=", {\[Theta]1, \
\[Theta]2, \[Zeta], amaa, res1}];

, {\[Zeta], 0, 1}];
, {\[Theta]2, Ceiling[(\[Theta]1 + 2)/2], 2 mmax}];
, {\[Theta]1, 1, mmax}];