Optimal investment problem with M-CEV model: closed form solution and applications to the algorithmic trading

TL;DR

This paper derives explicit solutions for an optimal investment problem under the M-CEV model using Laplace transforms, providing formulas for strategy analysis and extensions for algorithmic trading applications.

Contribution

It presents a closed-form solution for the M-CEV model's optimal investment problem and develops elementary and continued fraction approximations for practical analysis.

Findings

Explicit optimal strategy expressions using hypergeometric functions

Asymptotic and approximation formulas for parameter analysis

Extensions applicable to algorithmic trading strategies

Abstract

This paper studies an optimal investment problem under M-CEV with power utility function. Using Laplace transform we obtain explicit expression for optimal strategy in terms of confluent hypergeometric functions. For obtained representations we derive asymptotic and approximation formulas contains only elementary functions and continued fractions. These formulas allow to make analysis of impact of model's parameters and effects of parameters misspecification. In addition we propose some extensions of obtained results that can be applicable for algorithmic strategies.

Figures

Equations

$$d{B}_s$$

$$d{S}_s/S_s$$

$$\sigma (S_s,s)=a{S}^{\beta },\lambda (S,s)=b+c{\sigma }^2(S,s)=b+c{a}^2{S}^{2\beta },q(s)=q,r(s)=r,\alpha =r-q+b,$$

$$d{S}_s/S_s=\left[\alpha +c{a}^2{S}^{2\beta }\right]ds+a{S}^{\beta }d{W}_s.$$

$$d{X}_s=r(X_s-\pi S_s)ds+{\pi }_{s}d{S}_s.$$

$$J(X,S,t)=\underset{\pi }{\sup }\mathbb{E}[U(X_T)|X_t=X,S_t=S],$$

$$U(x)=\frac{x^{\gamma }}{\gamma }.$$

$$\left\{{\begin{array}[]{l}dS_{s}/S_{s}=\mu(V_{s},s)ds+\sigma(V_{s},s)dW_{s}^{1},\quad S_{t}=S,\\ dV_{s}=b(V_{s},s)ds+a(V_{s},t)dW_{s}^{2},\quad V_{t}=v,\end{array}}\right.$$

$$J(X,S,v,t)=\frac{X^{\gamma }}{\gamma }{f}^{1/\delta }(v,t),\delta =1+{\rho }^2\frac{\gamma }{1-\gamma }.$$

$$\left\{{\begin{array}[]{l}f_{t}+\frac{1}{2}a^{2}(v,t)f_{vv}+\left[b(v,t)+\rho\frac{\gamma(\mu(v,t)-r(t))a(v,t)}{(1-\gamma)\sigma(v,t)}\right]f_{v}+\frac{\gamma\delta}{1-\gamma}\left[\frac{(\mu(v,t)-r(t))^{2}}{2\sigma^{2}(v,t)}+(1-\gamma)r\right]f=0,\\ f(v,T)=1,\end{array}}\right.$$

$${\pi }^{\ast }(X,S,v,t)=\frac{X}{S(1-\gamma )}\left(\frac{\mu (v,t)-r(t)}{{\sigma }^2(v,t)}+\frac{\rho }{\delta }\frac{a(v,t)f_v(v,t)}{\sigma (v,t)f(v,t)}\right).$$

$$S=v,\rho =1,a(S,s)=S\sigma (S,s),b(S,s)=S\mu (S,s).$$

$$J(X,S,t)=\frac{X^{\gamma }}{\gamma }{f}^{1/\delta }(S,t),\delta =\frac{1}{1-\gamma }.$$

$$\left\{{\begin{array}[]{l}\mathcal{L}f\equiv f_{t}+\frac{a^{2}S^{2\beta+2}}{2}f_{vv}+\delta S[\alpha-\gamma r+ca^{2}S^{2\beta}]f_{v}+\frac{\delta(1-\delta)}{2a^{2}}\left[\left(\alpha-r\right)S^{-\beta}+cS^{\beta}\right]^{2}f+r\gamma\delta f=0,\\ f(v,T)=1,\end{array}}\right.$$

$${\pi }^{\ast }(X,S,t)=X\left(\delta \frac{\alpha -r+c{a}^2{S}^{2\beta }}{a^2{S}^{2\beta +1}}+\frac{f_S(v,t)}{f(v,t)}\right).$$

$$f(S,t)={\int }_0^{\infty }{f}_G(S,t;\xi )g(\xi )d\xi .$$

$$z=\frac{\Lambda }{S^{2\beta }},\tau =a^2{\beta }^2\Lambda (T-t),\Lambda =\frac{\sqrt{\delta }}{a^2|\beta |}\sqrt{{\alpha }^2-\gamma r^2},$$

$$f(S,t)={\int }_0^{\infty }{F}_G(z,\tau ;\xi )g\left({\left(\frac{\Lambda }{\xi }\right)}^{1/2\beta }\right)d\xi .$$

$$F_G(z,\tau ;\xi )=\frac{1}{2}\exp \left\{R\tau +Q(z-\xi )-\frac{(z+\xi )}{2}\coth (\tau )\right\}{\left(\frac{z}{\xi }\right)}^{\lambda +1/2}\frac{1}{\sinh (\tau )}{I}_{2\eta }\left(\frac{\sqrt{z\xi }}{\sinh (\tau )}\right),$$

$$\lambda =-\frac{1}{2}-\frac{1}{2\beta }\left(\frac{1}{2}-\delta c\right),\eta =\sqrt{{\left(\lambda +\frac{1}{2}\right)}^2+\frac{\delta (1-\delta )c^2}{4a^4{\beta }^2}},$$

$$Q=\frac{\delta (\alpha -\gamma r)}{\Lambda \beta a^2},R=\frac{r\delta }{a^2{\beta }^2\Lambda }-2Q\lambda -\frac{\delta (1-\delta )(\alpha -r)c}{\Lambda a^4{\beta }^2}.$$

$$f(S,t)={\int }_0^{\infty }{F}_G(z,\tau ;\xi )d\xi .$$

$${\int }_0^{\infty }{x}^{\mu -\frac{1}{2}}{e}^{-\alpha x}{I}_{2\nu }\left(2\beta \sqrt{x}\right)dx=\frac{\Gamma \left(\mu +\nu +\frac{1}{2}\right)}{\Gamma (2\nu +1)}{\beta }^{-1}{e}^{\frac{{\beta }^2}{2\alpha }}{\alpha }^{-\mu }{M}_{-\mu ,\nu }\left(\frac{{\beta }^2}{\alpha }\right),$$

$$\mathrm{Re}\left(\mu +\nu +\frac{1}{2}\right)>0.$$

$$f(S,t)=e^{R\tau +zB(\tau )}{D}^{\lambda }(\tau )\frac{\Gamma (\eta -\lambda +1/2)}{\Gamma (1+2\eta )}{e}^{-\frac{z}{2}A(\tau )}{\left(zA(\tau )\right)}^{\lambda }{M}_{\lambda ,\eta }\left(zA(\tau )\right),$$

$$A(\tau )=\frac{1}{2{\sinh }^2(\tau )[\coth (\tau )+Q]},B(\tau )=\frac{Q^2-1}{2[\coth (\tau )+Q]},D(\tau )={\sinh }^2(\tau )[\coth (\tau )+Q{]}^2.$$

$${\left(z\frac{d}{dz}z\right)}^n\left(e^{-z/2}{z}^{k-1}{M}_{k,\mu }(z)\right)=\frac{\Gamma (\mu +k+n+1/2)}{\Gamma (\mu +k+1/2)}{e}^{-z/2}{z}^{k+n-1}{M}_{k+n,\mu }(z).$$

$${\pi }^{\ast }(X,S,t)=X\left(\delta \frac{\alpha -r+c{a}^2{S}^{2\beta }}{a^2{S}^{2\beta +1}}+\left[B(\tau )+\frac{\lambda +\eta +1/2}{z}\frac{M_{\lambda +1,\eta }\left(A(\tau )z\right)}{M_{\lambda ,\eta }\left(A(\tau )z\right)}\right]\frac{dz}{dS}\right).$$

$$M_{\lambda ,\eta }(x)=e^{-x/2}{x}^{1/2+\eta }\Psi (\theta ,\omega ,x),\theta =1/2+\eta -\lambda ,\omega =1+2\eta$$

$${\pi }^{\ast }(X,S,t)=\frac{X}{S}\left[\frac{\delta (\alpha -r)/a^2-2\beta \Lambda B(\tau )}{S^{2\beta }}+\delta c+2\beta (\theta -\omega )\frac{\Psi \left(\theta -1,\omega ,\Lambda A(\tau )S^{-2\beta }\right)}{\Psi \left(\theta ,\omega ,\Lambda A(\tau )S^{-2\beta }\right)}\right].$$

$$\frac{\Psi (\theta -1,\omega ,x)}{\Psi (\theta ,\omega ,x)},\frac{M_{\lambda +1,\eta }\left(A(\tau )z\right)}{M_{\lambda ,\eta }\left(A(\tau )z\right)}.$$

Full text

1 Introduction

1.1 Motivation

Many academic papers about optimal investment problems assume that the asset price follows geometric Brownian Motion(GBM). However, there are a lot of empirical studies showing this simple model does not properly fit to real market data. Known drawbacks are the following: GBM model does not capture volatility smile/skew effects; they ignore the probability of the underlying’s default; the constant coefficients do not allow calibration of this model to the real term structure of interest rates and dividend yields etc. Our motivation is to extend the results of GBM models to a more realistic model. In order to obtain a more realistic fit to the market data we can use more sophisticated models based, for example, on Levy processes or on fractional Brownian motion. But although their dynamics are more realistic, these complicated models are not usually analytically tractable. Hence quantitative analysis is complicated and any qualitative analysis is impossible. We must try compromise between realistic modelling and the availability of analytical or quasi-analytical expressions.

In this paper we solve an optimization problem assuming the Modified Constant Elasticity of Variance (i.e. M-CEV) model for the asset’s price and a power utility over the final wealth for a finite horizon agent. This model was introduced in Heath and Platen (2002) and is a natural extension of the famous CEV model (see Cox(1975)). We choose this model for the following reasons: this model captures the volatility smile effect; allows non-zero probability of the underlying’s default (M-CEV process can touch zero while GBM is always positive); and it is analytically tractable. Also this model is applicable to algorithmic trading strategies because the M-CEV process has a mean-reversion property for some of the model’s parameters. Let us mention that the time-dependent extension of this model can be found in Linetsky and Carr (2006). For the M-CEV model we obtain a closed-form solution in terms of confluent hypergeometric functions. Despite the availability of many numerical solvers (i.e. PDE solvers or Monte-Carlo) explicit formulas are still relevant. There are several reasons to pursue a closed form solution: first, they show dependencies between model parameters and optimal policy, therefore we can obtain some non-trivial qualitative effects. Second, properly programmed closed-form solutions give faster and more efficient code than a lot of available numerical solvers(PDE solvers or Monte-Carlo). In addition, simple tractable models can serve as a benchmark in practical situations. Quite often, practitioners prefer to introduce ad hoc corrections to a simple model than to use a more involved model with a large number of parameters.

Another important point is the utility choice. There are some popular utility functions considered in the literature: logarithmic, power and exponential. Obviously, each utility gives a different optimal strategy that maximizes expected utility over the terminal wealth. It is well known that the optimal strategy in the case of a logarithmic utility does not depend on the time to the end of the investing period and the trading rules of an exponential utility investor is not sensitive to the current wealth (see Merton (1990)). In order to capture time and wealth dependencies we choose a power utility.

1.2 Previous research

There are a lot of papers about similar problems: T. Zariphopoulou (2001) considered the problem for stochastic volatility models and derived optimal policy as the solution of the parabolic PDE. Some closed-form solutions and asymptotic expansions for various models can be found in Kraft (2004), Chacko and Viceira (2005), Boguslavskaya and Muravey (2015). A detailed review of papers about closed-form solutions and asymptotics can be found in Chan and Sircar (2015).

Applications of the utility maximization problems to algorithmic trading were discussed in Boguslavsky and Boguslavskaya (2004), Liu and Longstaff (2000).

1.3 The main results and structure of the paper

The main result of this paper is the closed form solution for the expected utility maximization in the finite horizon with power utility and M-CEV model. We derive asymptotic and approximation formulas containing only elementary functions and continued fractions. The structure of this paper is as follows: first we define the problem. Then we present a closed form solution for the M-CEV model. This is followed by the algorithm of numeric implementation and an analysis of parameter misspecification. Applications of the obtained results to algorithmic trading strategies are then given. All proofs are in Appendix A.

2 Problem definition

2.1 Model setup

Consider a simple market consisting of a risk-free bond $B_{t}$ and a risky asset (i.e. stock) $S_{t}$. The bond and stock prices are driven by SDE:

[TABLE]

where $W_{s}$ is a standard Wiener process, $r(s)\geq 0$, $q(s)\geq 0$, $\sigma(S,s)>0$ and $\lambda(S,t)\geq 0$ are the time-dependent risk-free interest rate, the time-dependent dividend yield, the time- and state- dependent instantaneous stock volatility, and the time- and state- dependent default intensity, respectively.The M-CEV model has the following specifications:

[TABLE]

and defined by this corresponded SDE

[TABLE]

Let us mention that Heath and Platen considered model (2.3) with $c=1$. The case of $c\neq 1$ is not extension of original M-CEV model because this case can be reduced to the original model by a simple change of measure. We will use specification (2.3) with $c\neq 1$ to analyze the impact of parameter $c$ directly. The optimal investment problem can be treated in the general portfolio optimization framework. Assuming no market frictions and an absence of transaction costs, the wealth dynamics for a control $\pi_{s}$ is given by

[TABLE]

Here $\pi_{s}$ is the investor position in stock(i.e. the number of units of the asset held). We assume that there are no restrictions on $\pi_{s}$, so short selling is allowed and there are no marginal requirements on wealth $X_{s}$. We solve the expected terminal utility maximization problem for an agent with a prespecified time horizon $T$ and initial wealth $X_{0}>0$. The value function $J(X,S,t)$ is the expectation of the terminal utility conditional on the information available at time $t$ ($S_{t}=S$, $X_{t}=X$).

[TABLE]

where $U(X)$ is the power utility function

[TABLE]

2.2 Known results

In this section we provide some known results used later in this paper. The first result is about a reduction of the original problem (2.5) with power utility (2.6) to the Parabolic partial differential equation (PDE).

Theorem 2.1** (Zariphopoulou).**

Assume that the asset price process $S_{t}$ follows SDE

[TABLE]

where $W_{s}^{1}$ and $W_{s}^{2}$ are correlated Wiener processes with coefficient $\rho$ and the investor has power utility function (2.6). In these assumptions the value function (2.5) can be represented as(i.e. distortion transformation)

[TABLE]

Function $f$ is a solution of the linear parabolic PDE boundary problem

[TABLE]

Optimal policy $\pi^{*}(S_{t},X_{t},v_{t},t)$ is given in the feedback form

[TABLE]

It is easy to show that the T.Zariphopoulou result can be applied to the M-CEV model (2.3) by substitution

[TABLE]

Proposition 2.1**.**

For M-CEV model the value function $J(X,S,t)$ is given by

[TABLE]

Function $f$ solves Cauchy problem

[TABLE]

and optimal policy $\pi^{*}(X,S,t)$ is given by

[TABLE]

The main difficulty is to solve boundary problem (2.13). In the next section we present a closed-form solution of (2.13) in terms of confluent hypergeometric functions.

3 Main results

Consider the Cauchy problem (2.13) with arbitrary initial function $f(S,T)=g(S)$. It is known that its solution can be represented as a convolution product with Green function $f_{G}(S,t;\xi)$

[TABLE]

Using the Laplace transform method we obtain the explicit representation for Green function $f_{G}(S,t;\xi)$ in terms of Modified Bessel function $I_{\nu}(z)$ (for definition see Abramovitz and Stegun (1971)). Hence the solution of problem (2.13) can be easily obtained by application of formula (3.1) with initial function $g(S)\equiv 1$. For convenience we will use scaled space and inverse time variables $z$ and $\tau$:

[TABLE]

for function $f(S,t)$ we have the following representation

[TABLE]

In the next theorem we introduce explicit formulas for Green function $F_{G}(z,\tau;\xi)$.

Theorem 3.1**.**

Green function $F_{G}(z,\tau;\xi)$ is given by

[TABLE]

where $\lambda$, $\eta$, $R$ and $Q$ are constants

[TABLE]

[TABLE]

Hence the solution of boundary problem (2.13) can be represented as

[TABLE]

We can perform these integrations explicitly by using the following relation (see Gradshteyn and Ryzhik (1980), formula 6.643.2) between Modified Bessel function $I_{\nu}(z)$ and Whittaker function $M_{\lambda,\eta}$(z) (see Abramowitz and Stegun (1973))

[TABLE]

In the result we have the following formula for function $f$

[TABLE]

where $\Gamma(x)$ is the Euler gamma function and functions $A(\tau)$, $B(\tau)$ and $D(\tau)$ are given by

[TABLE]

The expression for $f_{S}/f$ is obtained by using differential rules for Whittaker functions (see Abramowitz and Stegun (1973))

[TABLE]

Hence the optimal policy $\pi^{*}(X,S,t)$ is

[TABLE]

Using the following relation between Whittaker function and Kummer function

[TABLE]

and compute derivative $dz/dS$ we obtain alternative formulas for $\pi^{*}(X,S,t)$:

[TABLE]

4 Numerics

4.1 Numerical algorithm

If we want to build any quantitative trading strategy based on the obtained results we should have a numerical algorithm to compute expression 3.10 (or 3.12) for any parameters. It consists of only elementary functions except the term

[TABLE]

Obviously, computation of these special functions is not a problem for packages such as MATLAB or Mathematica. However, production codes are mostly written in C++ and we can not use these packages. In this context we must provide fast and efficient computation for this non-elementary term in a C++ environment. There are libraries containing numerical algorithms for special functions (e.g. C++ GSL package has numerics for the Kummer confluent hyper-geometric function used in (4.1)). Hence we can compute (4.1) by the following scheme: if we have singularity in (4.1) we use asymptotic formulas(4.1), in other situations we use GSL. However, evaluations of Kummer functions can significantly slow down the computational speed of the algorithm and this approach is not suitable if speed is critical. In this section we provide a fast numerical scheme based on asymptotic expansions of term (4.1). The main idea is very simple. We construct two series expansions directly for term (4.1):

[TABLE]

and use first or second series depending on value of variable $x$. We compute approximation of series (4.2) recursively with the following stopping criteria: we stop evaluations if the difference between $N+1$ and $N$ truncated series is sufficiently small. In the next theorem we provide explicit formulas for coefficients $c_{s}$ and $d_{s}$.

Theorem 4.1**.**

Coefficients $c_{s}$ and $d_{s}$ in expansions (4.2) are defined by recursive formulas

[TABLE]

Figure (1) illustrates the convergence rate with fixed parameters $\theta$, $\omega$ and variable $x$ (left sub-figure) and accuracy of approximations (right sub-figure). For these tests we set $\theta=5.24$ and $\omega=1.42$. We illustrate convergence rate at the point $x=10$ and for accuracy illustration we set $N=80$ in expansions for small argument and set $N=8$ for large.

Remark 4.1**.**

If we truncate series (4.2) by 1 or 2 terms we obtain the following asymptotic for (4.1)

[TABLE]

Moreover, recursive application of the following relation between Kummer functions (see Abramowitz and Stegun (1973))

[TABLE]

turns out to the continued fraction representation

[TABLE]

where

[TABLE]

4.2 Computational speed test

In this section we present a computational speed benchmark. We perform tests on the standard laptop with an Intel Core i7-3537U processor and a GCC 6.3.1 C++ compiler. Both algorithms compute (4.1) for any $z$ by $10^{4}$ times. Our algorithm performs computations with accuracy $\epsilon=10^{-10}$. Parameters are set to $\theta=5.24$ and $\omega=1.42$. GSL routines have a predefined accuracy and we can-not change it. Let us also mention that our algorithm can compute (4.1) for large values (e.g. $x>732$) while GSL routines have overflow errors. Hence we need in some modifications of GSL routines (e.g. we can compute $e^{-x}\Psi(\theta,\omega,x)$ to avoid overflow) in case of large values of function argument. Figure (2) illustrates a comparison between our method based on formulas (4.2 - 4.3) and direct computation of the numerator and denominator in (4.1) using GSL routines. Parameters are set to The left sub-figure illustrates the computational speed of algorithms based on our formula for small arguments (red line) and the GSL algorithm (blue line). For $0<x<1$ our algorithm is faster than GSL, but for $1<x<4$ GSL is faster. Let us note that if we change $\epsilon$ we will have other results. The right sub-figure illustrates speed’s comparison in case of large argument $x$. In this case our solution is faster than GSL at whole segment $40<x<732$. For $732<x<1000$ GSL routines can not evaluate function value. We suggest the low speed of GSL routines may be caused by exponential grow of the Kummer functions for large arguments (it also can cause overflow errors). The source C++ codes can be found at GitHub repository (see link in the references).

4.3 Parameters misspecification

This section contains several numerical examples that illustrate optimal strategy and the effects of parameters’ misspecification. Figure 3 demonstrates the wealth dynamics for 3 different investment strategies.The first strategy (PnL is colored by red) consists of only bond investments. We have invested all of the initial wealth $X_{0}=100$ in bond $B_{t}$ with initial value $B_{0}=X_{0}$ and interest rate $r=0.04$. The second strategy (blue line) consists of only stock investments. The stock process has initial value $S_{0}=100$, average return $\alpha=0.045$, volatility $a=0.4$, default intensity $c=0.8$ and skewness $\beta=-0.4$. In these strategies we do not have any portfolio re-balancing during the whole investing period $T=1$. Positions in the third strategy (yellow color) are defined by formulas 3.10. The investor’s risk aversion is $\gamma=-4$. The second figure 4 illustrates the terminal wealth distribution with true and misspecified parameters. For these tests we have $10^{6}$ simulations to compute terminal wealth distributions. These examples show that calibration errors in average return $\alpha$ and skewness $\beta$ are more critical than errors in volatility level $a$ and default intensity $c$.

5 Applications to the algorithmic trading

In this section we propose a statistical arbitrage strategy based on our obtained results. Consider an arbitrageur trading a mean-reverting asset. Suppose that the trader knows the ’fair’ mean price of the asset (i.e. long term mean) and he knows that price will be return to this mean price. Generally, in this framework a trader can make profit by take a long position when the asset is below its long-term mean and a short when it is above. The question is in the size of the trader’s position and how the position should be optimally managed depending on the price process parameters and trader’s current wealth. This optimal trading problem can also be treated in the general portfolio optimization framework and it corresponds to the zero interest rates case i.e. we must set $r(t)=0$ in all formulas. Therefore the wealth for a control $\pi$ is given by

[TABLE]

here $\pi_{s}$ is the trader position in the mean-reverting asset. It is well known that the original M-CEV price process $S_{t}$ can be mean-reverting if $\alpha<0$. Without loss of generality, we consider the case of a square-root diffusion process which corresponds to parameters

[TABLE]

This leads to the following mean reverting process

[TABLE]

Parameter $\kappa$ is the reversion speed, $\bar{S}$ is long-term mean and $a$ is the volatility level.

Proposition 5.1**.**

For the process (5.3) the value function and optimal control allows representation (3.8) and (3.10). The parameters have the following representation

[TABLE]

Functions $A(\tau)$, $B(\tau)$ and $D(\tau)$ are given by

[TABLE]

Optimal position is

[TABLE]

If we intend to perform trading strategies based on optimal control of a square-root process we must properly construct the mean-reverting asset. One of the standard approaches to mean-reversion trading is pair trading. In this case we construct a mean-reverting asset as a difference(i.e. spread) between two co-integrated assets. In almost all cases this spread has zero long term mean. Therefore we can not consider a square-root process for pair trading because it has only positive values if $2\kappa\bar{S}>a^{2}$. This leads us to change the difference to another mean-reverting asset. We propose to make a strategy for FX rates. They are always positive and can be mean-reverting. Hence we can model it using a square-root process. Figure (5) illustrates a trading strategy based on USD/CAD historical data. We consider daily data over 6.5-year time period from 01/01/2011 to 26/06/2017. We calibrate parameters of square-root process $\kappa$ ,$\bar{S}$ and $a$ on 01/01/2011-01/07/2016 daily rates. We obtain following values:

[TABLE]

We test our strategy on the 01/07/2016-26/06/2017 time period, $T=0.9961$. The investor’s risk aversion is set to $\gamma=-7$, initial wealth is $X_{0}=1000$. Strategy based on formula (5.6) has $4.33\%$ return, $0.6464$ Sharpe ratio and $-6.54\%$ maximum drawdown, while buy and hold (blue line) strategy has $2.61\%$, $0.3911$ and $-7.03\%$.

Appendix A Appendix: Proofs

A.1 Theorem 3.1

It is easy to show that unknown Green function $F_{G}(z,\tau;\xi)$ can be represented as

[TABLE]

where $z$ and $\tau$ are defined in (3.2) and function $h(z,\tau)$ solves Cauchy problem ($\delta(z)$ is the Dirac delta function)

[TABLE]

Let $G(z;\zeta)$ be a Laplace transform of the function $h(z,\tau)$:

[TABLE]

It turns out to the following ODE for function $G$

[TABLE]

The homogeneous equation in (A.3) is called Whittaker equation and have two linearly independent solutions, namely $M_{-\zeta/2,\eta}(z)$ and $W_{-\zeta/2,\eta}(z)$ (see Abramowitz and Stegun (1973)).It is easy to show that the solution of non-homogeneous problem (A.3) can be represented as

[TABLE]

Using this relation between Whittaker functions and modified Bessel function(see Gradshteyn and Ryzhik (1980), formula 6.669.4)

[TABLE]

we obtain new formula for $G(z,\zeta)$

[TABLE]

Next we introduce new integration variable $\nu$

[TABLE]

In the result we have

[TABLE]

Inverting the Laplace transform, we recover the formula for $h(z,t)$

[TABLE]

where $N$ is a number such that all residues of the integrand are to the right of it. Using the well-known representation of Dirac function

[TABLE]

and changing the order of integration in (A.5), we get

[TABLE]

Note, that $\tau\geq 0$. Thus, we can complement the range of integration in (A.6) to the whole line, and, using the definition of Dirac’s function, namely $\int_{-\infty}^{\infty}\delta(\zeta-z)u(\zeta)d\zeta=u(z)$ for any continuous $u$, we get the main formula (3.4) for $F_{G}$.

A.2 Theorem 4.1

Consider the quotinent of two series

[TABLE]

It is equivalent to

[TABLE]

or

[TABLE]

Hence the coefficients $c_{s}$ solve the following linear system:

[TABLE]

For $c_{s}$ in 4.3 we use the following definition of Kummer function (see Abramovitz and Stegun (1972))

[TABLE]

and for $d_{s}$ we use the asymptotic of Kummer function for large argument (see Abramovitz and Stegun (1972))

[TABLE]