Non-stationary Feller process with time-varying coefficients

We study the non-stationary Feller process with time varying coefficients. We obtain the exact probability distribution exemplified by its characteristic function and cumulants. In some particular cases we exactly invert the distribution and achieve the probability density function. We show that for sufficiently long times this density approaches a Gamma distribution with time-varying shape and scale parameters. Not far from the origin the process obeys a power law with an exponent dependent of time, thereby concluding that accessibility to the origin is not static but dynamic. We finally discuss some possible applications of the process.


I. INTRODUCTION
The Feller process is a one-dimensional diffusion with linear drift and linear diffusion coefficient vanishing at the origin. One of its most distinctive characteristics is that it never attains negative values. This, along with the fact that the process amounts to a continuous representation of a number of branching and birth-dead processes, have made Feller process an ideal candidate for modeling many phenomena in physical and social sciences [1,2]. Examples span from theoretical biology and neurobiology [3][4][5][6][7], population growth [8][9][10][11][12], radiation physics [13] and even economics and financial markets [14][15][16][17][18][19]. From a more formal point of view, the process is also valuable because some one-dimensional diffusions can be mathematically transformed into the Feller process, as it was proved by Capocelli and Ricciardi years ago [20].
In its original formulation [21] the parameters are time independent and the process is stationary. This means that it is invariant under time translations and that, as time progresses, the probability distribution tends towards a stationary value given by the Gamma distribution [1]. However, as countless situations show, stationary processes are a convenient idealization and many empirical studies indicate that real data susceptible to be modeled by diffusion processes are not fully stationary [22].
There are, however, very few attempts at addressing this question for the Feller process (as well as for other diffusions). One exception is, to our knowledge, the recent work by Gan and Waxman [2] where the authors, using the method of spectral decomposition, obtain an explicit solution for the distribution function of the process with time varying coefficients. Although they achieve the solution only in the special and singular case when the linear drift also vanishes at the origin. This is rather restricting because, as we will see below, the value of the drift at the origin determines where the process tends to * Electronic address: jaume.masoliver@ub.edu in the course of time. Fixing this value to zero limits the scope and applications of the solution (even though it is valuable in population dynamics when one is faced with possible extinctions [9,12]).
In this paper we address the study of the Feller process in the most possible general form. The main result is knowing the exact probability distribution materialized by the Laplace transform of the probability density function, that is, by its characteristic function. The transformed density is inverted in some particular but relevant cases (including that of the original process). In addition we also get some asymptotic expressions valid for large to moderate times showing the emergence of the Gamma distribution but with time-varying shape and scale parameters.
Another significant aspect is the behavior of the process near the origin which, in turn, determines whether or not the origin is accessible. In the stationary case accessibility depends on the value taken by a constant ratio between two parameters governing the deterministic and random components of the dynamics of the process [1,21]. We will see that for the non-stationary process the behavior of the probability density function near the origin is given by a power law with a time-dependent exponent which implies that the likelihood of attaining the origin is not fixed but changing with time.
The paper is organized as follows. In Sec. II we study the dynamics and the general properties of the nonstationary process. In Sec. III we obtain the exact probability distribution through the knowledge of its characteristic function and cumulants. Section IV is devoted to obtaining the exact form of the probability density function for some special cases. In Sec. V we get approximate expressions valid for moderate to long times. Section VI focusses on the behavior of the probability density near the origin and its attainability by the process. A short summary of the main results with few possible applications is in Sec. VII and more technical details are in several Appendices.

II. DYNAMICS AND GENERAL PROPERTIES
In its most general form the Feller process can be described by the following Langevin equation: where ξ(t) is the derivative of the standard Wiener process. That is, ξ(t) zero-mean Gaussian white noise All stochastic differential equations of this paper are interpreted in the sense of Itô.
Although from a mathematical point of view the parameters of the process are given by arbitrary functions of time, from a physical point of view it is more convenient to assume that α(t) > 0, β(t) ≥ 0 and k(t) > 0 are positive and smooth functions. This will be one of our basic assumptions for the rest of the paper.
Note that α(t) has dimensions of (time) −1 and sets a dimensionless time scale defined by where t 0 is any initial time. Observe that t = t 0 corresponds to τ = 0 and also that τ = τ (t) is an increasing function of the original time scale because dτ /dt = α(t) > 0 (t ≥ t 0 ). In other words, "future maps into future" (and past into past). For the rest of the paper we also assume that α(t) is such that That is to say, t → ∞ implies τ → ∞ and vice versa. This is our second and last restriction on the parameters of the process.
In terms of this new time scale we write the stochastic equation (1) as where t = t(τ ) is the inverse function implicitly defined in Eq. (2). Defining the rescaled parameters and the rescaled noise the Langevin equation for X(τ ) reads The rescaled noise η(τ ) is zero-mean Gaussian white noise: Indeed, that η(τ ) is Gaussian noise with zero mean is evident from Eq. (4). We next show that η(τ ) is delta correlated. From definition (4), we have Let us recall the standard property of the delta function, is any differentiable function and x 0 is the solution to g(x 0 ) = 0. Then, since t(τ 1 ) = t(τ 2 ) is equivalent to τ 1 = τ 2 [we assume that t(τ ) is an univariate function], we write where the prime refers to the derivative and t(τ ) and τ (t) are inverse functions for which t ′ = 1/τ = 1/α [cf. Eq.
(2)]. Hence which proves Eq. (6). In the representation given by Eq. (5) the process has a state and time dependent diffusion D(X, τ ) = 2D(τ )X, which for large values of X enhances the effects of noise (with an intensity depending of time) while as X approaches zero the effect of noise diminishes. Hence, when the process reaches the origin, the drift drags it toward the asymptotic value: Indeed, near the origin the process (5) approximates to the deterministic process dx/dτ = −[x − m(τ )] and from its solution we readily see that x(τ ) → m ∞ as τ → ∞. Consequently, if m(τ ) ≥ 0 the process, starting at some positive value, cannot reach the negative region and the Feller process is nonnegative [otherwise the noise term in Eq. (5) would become imaginary].
Moreover, the asymptotic value m ∞ coincides with the stationary average of the process, In effect, the formal solution of Eq. (5) such that X(0) = x 0 is As is well known, in the Itô interpretation the output process X(τ ) and the input white noise η(τ ) are uncorrelated. Hence which can be written as and the limit τ → ∞ proves Eq. (8). Let us observe that when m(τ ) = m is constant, m ∞ = m, and the process tends to m. This is the reason why in many setting (such as, for example, economics [14][15][16][17][18][19]) m is called the "normal level" of the process. We, therefore, see that for the non-stationary Feller process, and under rather general circumstances (i.e., m(τ ) ≥ 0), the origin is a singular boundary that the process cannot cross. A closely related problem is whether or not the origin is attainable, in other words, whether the value X = 0 can or cannot be reached. This is a key issue in many practical situations (for instance, in population dynamics where attaining the origin signifies annihilation). We will prove later on that the answer to this question is time dependent and determined by a varying parameter balancing deterministic motion and fluctuations. The problem of classifying the different types of boundaries appearing in diffusion processes was thoroughly studied by Feller himself in the early 1950s and we refer the reader to the literature for a more complete account on this topic [23][24][25].

III. PROBABILITY DISTRIBUTION AND CUMULANTS
Our main objective is knowing the probability distribution of the non-stationary Feller process (1). In this section we will determine the distribution of probability through the characteristic function (CF) as well as the cumulants.
In the time scale defined by Eq. (2), the probability density function (PDF) of the process X(τ ), with initial condition Recall that x = 0 is a singular boundary that the process cannot cross. A sufficient condition for this to happen is that the probability flux through x = 0 is zero. We thus look for solutions of the initial-value problem (9) and (10) that meet this condition: where is the flux of probability through x at time τ . Note that in terms of the flux the FPE (9) can be written as We will now proceed to obtain the exact probability distribution via the characteristic function of the process, the latter defined as the Laplace transform of the PDF with respect to x (bear in mind that x ≥ 0): The Laplace transform of Eq. (13) under condition (11) reads and substituting for Eq. (12) we have We, therefore, obtain the following partial differential equation of first order with initial condition [cf. Eq. (10)] The exact solution to this initial-value problem can be obtained by the method of characteristics [26]. This is detailed in Appendix A with the result and the function φ(τ ) by We next write the probability distribution in terms of the original time scale t. Note that obtaining the CFp(σ, t|x 0 , t 0 ) certainly amounts to substituting in Eq. (17) the dimensionless time scale τ by the function τ (t) given in Eq. (2). However, we also have to write θ τ (ξ) [cf. Eq. (18)] as a function of t which is rather involved.
In the Appendix B we show that where, and t(ξ) is defined by φ[t|t(ξ)] = ξ. That is, Equation (21) is our main result and completely determines, through the characteristic function, the probability distribution of the non-stationary Feller process. However, the exact analytical inversion of this equation, thus obtaining the PDF, seems to be beyond reach except for few special cases to be detailed in the next section. It is, nonetheless, possible to get various approximate forms (see Secs. V and VI) as well as obtaining exact expression for the cumulants of any order.

IV. SOME SPECIAL CASES
We will now obtain exact expressions of the probability density function (PDF) in three particular instances.
(i) Suppose first that the drift of the processes also vanishes at the origin. In other words, β(t) = 0 which is equivalent to m(τ ) = 0 [cf. Eq. (3)]. This is a rather peculiar situation because the normal level coincides with the singular boundary meaning that the process tends to that singularity as time progresses (see discussion in Sec. II). This is the case treated in Ref. [2]. Now the parameter θ τ (ξ) defined in Eq. (23) also vanishes and the CF (21) readŝ where φ(t, t 0 ) is given by Eq. (22) and In the Appendix C we show that the Laplace inversion of Eq. (27) is where I 1 (z) is a modified Bessel function. Equation (29) agrees with the PDF obtained by Gan and Waxman who have stressed the singular character of the origin [2] something absent in Feller's original formulation [21]. Indeed, the term involving the Dirac function reflects the singularity of the origin. For, as we readily see from Eq. (29) p(x, t|0, t 0 ) = δ(x), which implies that if the process starts at x 0 = 0 remains there forever.
(ii) Suppose now that β(t) and k 2 (t), are proportional to each other for all t ≥ t 0 . In this case the ratio θ t (ξ) defined in Eq. (23) is constant: where θ > 0 is a positive constant. Note incidentally that in the dimensionless time scale τ this case corresponds to having the normal level m(τ ) proportional to the diffusion coefficient D(τ ) for all τ ≥ 0.
The characteristic function, Eq. (21), now readŝ The integral in the exponential is easily evaluated and we havê p(σ, t|x 0 , t 0 ) = 1 The Laplace inversion of this expression is sketched in the Appendix D and the exact PDF is where I θ−1 (z) is a modified Bessel function of order θ −1.
(iii) Our last special case for which the PDF can be obtained exactly is when all parameters are time inde-pendent: which corresponds to the original Feller process [21]. Now and since θ is constant this case is a particular case of (ii) above. Therefore, the PDF of the stationary process will be given by Eq. (33) with [cf. Eq. (20)] Note that since the process is now stationary we can set t 0 = 0 without loss of generality.

V. LONG-TIME ASYMPTOTICS
In this section we get some interesting approximations to the PDF which are effective for sufficiently long times.
Let us first remember that one of our basic assumptions, as expressed in Sec. II, is that the time-dependent parameter α(t) > 0 is such that if t is long so is the new time scale τ defined in Eq. (2). In other words, t → ∞ implies τ → ∞.
We thus start off with the CF in the form given by Eq. (17): Since τ is large, the definition of θ τ (ξ) given in Eq. (18) allows for the following approximation (τ ≫ 1), which enables us to estimate the integral appearing in the CF as Hence, the CF is approximated bŷ (τ ≫ 1). Note that this approximate CF has de same form than that of Eq. (32). Therefore, the approximate PDF would be given by Eq. (33) with the constant parameter θ replaced by the time-varying θ(τ ).
We also observe that, within the same degree of approximation for which Eq. (36) holds, we may write [cf. Eq. (20)] which, after neglecting the exponentially small term e −τ , yields If, as we did in obtaining Eq. (39), we also neglect exponentially small terms in Eq. (38), we get (τ ≫ 1), where, since time is large, the dependence on the initial value x 0 has vanished. On the other hand,

and using the Laplace inversion formula [27]
(γ > 0) we obtain for τ ≫ 1 the following approximation The PDF in real time, p(x, t), will be given by Eq. (41) after replacing τ by the function τ (t) given in Eq. (28). Moreover, from Eqs. (3) and (37) we obtain and Therefore (t ≫ t 0 ). We thus see that for long times the PDF is approximated by a Gamma distribution with shape and scale parameters, θ(t) and D(t) respectively, dependent of time.
Let us remark that the asymptotic approximation given either by Eq. (41) or Eq. (44) does not apply to the singular case β(t) = 0 [which implies θ(t) = 0]. Indeed, in such a case Eq. (40) readsp(σ, t) ≃ 1 which after Laplace inversion yields Reflecting the fact that the homogeneous drift, −α(t)x, drags the process towards the origin as t → ∞. Let us incidentally note that the same result is indeed achieved from the exact solution (29) in the limit t → ∞. We end this section by getting the asymptotic expression of cumulants and moments. At first sight the most direct way of obtaining the long-time form of cumulants is starting from the general expression given in Eq. (26) and taking there the limit t → ∞. However, it turns out to be much simpler to start off from the asymptotic expression of the CF given by Eq. (40) and apply the definition (25). The asymptotic cumulants thus read (n = 1, 2, 3, . . . ). The connection between cumulants and moments becomes very involved as n increases [25]. It is again simpler the direct calculation of moments based on Eq. (40). Indeed, since and long-time moments are only determined by the normal level of the process. In the opposite case when θ(t) ≪ 1, long-time moments are determined by the diffusion coefficient: Note that when θ(t) → θ tends to a finite constant [note that case (ii) of previous the section is a particular instance] also leads to moments dominated by the diffusion coefficient:

VI. NEAR THE ORIGIN
Another situation where it is possible to know an approximate PDF is when the process is near the origin. The behavior of the process for small values of x is rather significant in many fields such as, for instance, the firing of neurons [5] and also in population dynamics where, as we mentioned before, attaining x = 0 means extinction [8,9,12].
The problem is closely related to the question of whether or not the non-stationary Feller process can access the singular boundary x = 0. For the stationary process this problem was addressed by Feller himself who, after analyzing the PDF near x = 0, concluded that if the constant ratio θ = 2β/k 2 ≤ 1 the origin is an accessible boundary, while if θ > 1 it is not [21] (see also Ref. [1]). In the stationary process the accessible character of x = 0 is time independent and fixed forever once we know the parameters of the model. We will now show that for the non-stationary process attaining the origin is not static but time dependent.
Let us incidentally note that the question of reaching the origin would be more elegantly addressed by evaluating, instead of the PDF, the hitting probability of first reaching x = 0, as we recently did for the stationary process [1]. However, obtaining hitting probabilities needs the knowledge of first-passage time densities, something rather involved for non-stationary processes. We will, therefore, analyze the behavior of the PDF near the origin as was done by Feller [21] .
The starting point of our asymptotic analysis is the characteristic function in the form given by Eq. (21) and we will use this exact expression for knowing the behavior of the PDF near the origin. The later inference is achieved after using Tauberian theorems which show that the small x behavior of p(x, t|x 0 , t 0 ) is determined by the large σ behavior of its Laplace transform [28].
Following this procedure we prove in the Appendix E that for large values of σ we havê where A(t|x 0 , t 0 ) is defined in Eq. (E7) of Appendix E and θ(t) in Eq. (42). The CF is thus approximated bŷ p(σ, t|x 0 , t 0 ) ≃ A(t|x 0 , t 0 ) σ θ(t) as σ → ∞. Using Tauberian theorems [28] we see that the behavior of the PDF near the origin will be given by the Laplace inversion of this expression: We, therefore, see that that the PDF near the origin follows a power law with a time-varying exponent and that at x = 0 the PDF is different from zero only if θ(t) ≤ 1: Following Feller [21] (see also [1]) we conclude that if θ(t) ≤ 1 the origin is an accessible boundary while if θ(t) > 1 it is not. Therefore, attaining the origin is not fixed but varies with time.

VII. CONCLUDING REMARKS
We briefly summarize the main findings of the paper. In this work we have concentrated on the study of the non stationary Feller process with time-varying coefficients: where W (t) is the standard Wiener process. We have assumed: (i) α(t) > 0, β(t) ≥ 0, and k(t) > 0 are positive and smooth functions of t ≥ t 0 where t 0 is an arbitrary initial time, and (ii) α(t) is such that Under these rather general conditions on the coefficients, which ensure a broad applicability of the process, we have been able to obtain the exact probability distribution exemplified by the characteristic function of the process, as well as cumulants of any order [cf. Eqs. (21) and (26)].
We have been able to exactly invert the CF to get the PDF in a few special cases and obtain some relevant approximations as well. Thus, for sufficiently long times the PDF approaches the Gamma distribution, (t ≫ t 0 ), with the shape of the distribution determined by and the scale by both parameters depending of time. Another situation where it has been possible to get an approximate PDF is when the process is not far from the origin. Thus, for small values of the state variable x, the PDF obeys a power-law with a time-varying exponent, where θ(t) > 0 is the ratio between the time-varying deterministic normal level β(t) and the strength of fluctuations k 2 (t)/2. Finally, the behavior near x = 0 also determines the question, rather significant in some settings, of whether the origin is accessible or not . We have proved that for strong fluctuations such that k 2 (t)/2 ≥ β(t) the origin is accessible while for milder fluctuations [i.e., k 2 (t)/2 < β(t)] it is nor. The special character of the origin is now dynamical rather than static (as was the situation of the stationary case). We finish this work by mentioning few possible physical applications of the non-stationary process. One field is econophysics where the Feller process is one of the most popular models for the volatility [16,17]. As a first approximation the normal level of the volatility -i.e., the parameter m(t) ≡ β(t)/α(t) [cf. Eq. (3)]-is assumed to be constant. However, there are evidences of seasonality in the behavior of volatility [29,30] which could be modeled assuming a periodic function of time as normal level.
A similar situation appears in the price of commodities [31,32].
We can find other potential applications in anomalous diffusion problems. In the broad and complex field of fractional dynamics [33,34], the non-stationary Feller process may represent a complementary and rather simple approximation to the problem that is not based on the fractional Brownian motion as driving noise. Thus, for instance, we have shown in Sec. V that when θ(t) ≫ 1 (t → ∞) the asymptotic mean-square value of the process is X 2 (t) ≃ m 2 (t) [cf. Eq. (47)]. If we further assume that as t → ∞ the normal level follows a power law m(t) ∼ t γ/2 (γ > 0) then the process shows the anomalous diffusion behavior: It is quite interesting that, under the assumption θ(t) ≫ 1, the anomalous behavior depends only on the normal level m(t) = β(t)/α(t) and it is independent of the noise intensity k 2 (t). Other different situations may indeed appear depending on the log-time behavior of the parameters α(t), β(t) and k(t) as we can see from the discussion at the end of Sec. V [cf. Eqs. (48) and (49)].