A Stochastic Differential Equation Inventory Model

Inventory for an item is being replenished at a constant rate whilst simultaneously being depleted by demand growing randomly and in relation to the inventory level. A stochastic differential equation is put forward to model this situation with solutions to it derived when analytically possible. Probabilities of reaching designated a priori inventory levels from some initial level are considered. Finally, the existence of stable inventory states is investigated by solving the Fokker–Planck equation for the diffusion process at the steady state. Investigation of the stability properties of the Fokker–Planck equation reveals that a judicious choice of control strategy allows the inventory level to remain in a stable regime.


Introduction
It has been recognized for some time that the demand for some items may be proportional to the inventory on display. Baker and Urban [1] argued that the demand rate of an item is of a polynomial functional form, dependent on the inventory level. A thorough review of such demand models has been carried out by Urban [2]. In his 2005 article Urban stated that two distinct functional models for the demand rate were dominant, Type I Models, where the demand rate advanced on the initial inventory alone, and Type II Models, where the demand rate was a continuous (predominantly a power) function of the inventory level. Tsoularis [3] used the type of demand function in this article to solve an optimal inventory control problem but no deep analysis of the stochastic differential equation was undertaken.
In this work we propose that the demand growth rate, d, is a continuous function of the inventory level, x, assuming the quadratic form: where d 1 and d 2 are both positive constants that establish the concavity (d (x) < 0) of the function (1). The growth in demand is at its most rapid for small inventory levels, governed primarily by the coefficient d 1 , rising gradually at a decreasing rate (d (x) > 0, d (x) < 0),  2 . This behaviour of the demand rate is a feature of the Type II Models reviewed by Urban [2]. However, when the inventory exceeds d 1 2d 2 , the growth in demand declines rapidly (d (x) < 0, d (x) < 0) and ceases altogether (d(x) 0) when x d 1 d 2 . This constitutes a departure from Type II Models and allows the realistic possibility of saturation in demand when the product inventory reaches a sufficiently high level.
We introduce next an element of stochasticity in the demand growth Eq. (1) by allowing the dominant growth parameter, d 1 , to be a random variable that evolves according to whered 1 is the mean value, η(t) is the continuous time Gaussian white noise and σ is the diffusion coefficient measuring the intensity of the disturbance. In discrete time, white noise is a sequence of independent uncorrelated random variables. In continuous time, however, the autocorrelation is the Dirac delta function, E[η(t)η(t + s)] δ(s). The white noise, although not an actual physical process, is a useful approximation to physical situations where noise is inherently present in the dynamics of a process [4]. The stochastic equivalent of the demand growth form (1) is The change in the actual demand induced by d(x) in an infinitesimal interval, dt, is where η(t) is symbolically written as the derivative of Brownian motion,w(t), as the Wiener process is nowhere differentiable. In a small time interval, dt, demand grows by d(x)dt and the inventory is replenished at a rate udt, so the infinitesimal inventory change is (u-d(x))dt. Using (4) we can formulate the following stochastic differential equation (SDE):

The Stochastic Differential Inventory Equation
By making the notational substitutions, α d 1 and D d 1 d 2 , we may rewrite (5) as the SDE: so that the mean demand growth,d 1 x − d 2 x 2 , adopts the more familiar logistic form, αx 1 − x D . We shall assume an initial inventory, x 0 ∈ [0, D) and that the control variable, u, assumes a constant value throughout some interval [0,t]. Moreover (6) will be valid for x ∈ [0, D) only, that is, the inventory will obey (6) so long as x is bounded form above by D, and (6) will be no longer valid for x ≥ D. The mean value of the random parameter d 1 , α, is the demand growth rate per unit of inventory (the random parameter whilst the inventory is low, and D is the level of inventory that places an upper limit on demand growth, so no further growth in demand is possible beyond this value. If demand saturation is unlikely to occur at small values then D → ∞ and demand grows linearly with inventory, which in this case evolves thus: Strong and weak solutions to SDEs are subject to certain conditions. In strong solutions, the Brownian motions are given on a given probability space whereas in weak solutions the Brownian motion is chosen. The existence and uniqueness of a continuous path strong solution to SDE (6) is subject to the following two conditions [5]: (i) The Lipschitz condition: for some constant independent of time t, L, in some time interval, [0, T ] and x 1 , . This is essentially a smoothness condition which is fulfilled in the sufficient as the functions, u − αx 1 − x D and σ x, are continuously differentiable. (ii) The growth condition: in some finite time interval, [0, T ]. This condition is imposed to prevent the solution to (6) becoming infinite in [0, T ], which is always the case when u and σ are bounded from above by L.
Provided these existence and uniqueness conditions are met, for an initial condition, a strong solution to (6) exists as a continuous path. The Itô solution is a Markov process given by the integral equation The first integral in (8) is a standard Riemann integral and the second is an Itô integral, whose convergence is interpreted in the mean square sense (L 2 ), lim

Solution to SDE (6)
In this section a solution to the temporally homogeneous process (6) is presented for two distinct cases: (i) when the order rate, u, is a nonzero constant, and (ii) when u 0.

Solution to SDE (6) with u 0
First introduce the integrating factor: Define y(t) x(t)F(t), with y(0) x 0 , and proceed to derive the random differential equation This is a Riccati equation with a random coefficient, − α F D , of the nonlinear term, and a random forcing function, Fu. We shall not make any further attempt to solve Eq. (9) in this article.
Directly from (6) we obtain a differential equation for the mean inventory, E[x(t)]: To derive the equation for the variance of the solution to (6) we use Itô's formula: (10) and (11) are the first two differential equations in a recursive scheme of differential equations involving higher moments, E[x n (t)], n ≥ 3, obtained by repeated application of Itô's formula on d(x n (t)).

Solution to SDE (6) when the Noise is Small
In many practical applications, the diffusion parameter, σ , is small. In such cases, it is reasonable to assume that the solution to the SDE will be a stochastic perturbation of the deterministic solution as σ → 0. We assume a solution to (6) of the form where y(t) is the solution to the deterministic differential equation which can be obtained by separation of variables. The drift term, u − αx 1 − x D , can be expanded via a Taylor series around the solution y(t): The diffusion term, σ x, can also be expanded as a power series: We substitute next the expansions (13) and (14) in (6) and equate coefficients of like powers of σ to obtain an infinite set of stochastic differential equations. We write below only the first two which are often adequate in practice: SDE (16) is a time-dependent Ornstein-Uhlenbeck process whose solution is the first order linearization of (6) around the deterministic solution found directly from (15). The solution to (16) with the obvious initial condition, x 1 (0) 0, is The Itô integral solution (17) is a Gaussian random variable, and as such it can assume any value with finite probability. The validity of the series expansion is investigated in [6] where it is shown that it is asymptotic,

Solution to SDE (6) with u 0
When u 0, (6) is a Bernoulli equation which is transformed via the substitution z 1 y to the standard form and finally, by virtue of The solution is not a Gaussian process as the Wiener process, w(t), appears in the exponent.

Solution to SDE (7)
In this section a solution to the temporally homogeneous SDE (7) is presented for two distinct cases: (i) when the order rate, u, is a nonzero constant, and (ii) when u 0.

Solution to SDE (7) with u 0
As in "Solution to SDE (6) with u 0" section, first introduce the integrating factor, F(t) exp σ 2 t 2 − σ w(t) , then proceed along the same lines as before we arrive at the final solution: , for a normally distributed random variable Z, is applicable here. The expected value of x(t) is then given by the following formula, which is independent of σ : As t → ∞,, The average inventory will be increasing in time if The variance can be found by first solving the differential equation using the integrating factor, e 2α−σ 2 t . The solution is then Solution to SDE (7) with u 0 When u 0 (7) reduces to: a geometric Brownian motion with the well-known solution,

The Boundary at x 0
The inventory, x =0, is an intrinsic boundary of the diffusion process (6) because the diffusion coefficient, σ x, vanishes there. Its classification is dependent on the integrability of the scale function, S(ξ ) Dσ 2 ξ , and the value of the scale function is As the scale function is divergent at x 0, x 0 is a natural boundary according to the Russian literature classification scheme [7]. According to another classification scheme proposed by Feller [5,8], one classifies the point, x 0, based on the convergence of the integral The above integral is convergent, and according to Feller, x 0 is an entrance boundary. An entrance boundary cannot be reached from the interior of the state space, that is the inventory cannot vanish in finite time from some initial value, x 0 0. The inventory however, can start from x 0 0 and quickly build up to nonzero values.
In the absence of any orders, u 0, and The scale function converges in the vicinity of x 0. The boundary, x 0, is attracting and the inventory never attains the boundary zero in finite time.

First Passage Times and Probabilities of Exit Through Absorbing Barriers
In this section we look at how long the inventory, initially at x 0 at time t 0, remains in the interval (x l , x r ), which is assumed to contain x 0 , x l < x 0 < x r . By erecting artificial absorbing barriers at x l and x r we investigate the probability that the inventory crosses over either barrier, x l or x r , and the mean first passage time to x l or x r . Define π l probability of exit through x l , π r probability of exit through x r ,

Passage Probabilities for SDE (6) with D < ∞
The solution to the ordinary differential equation yields the probabilities of exit through either x l or x r , starting from x 0 , with the boundary conditions: π l (x l ) 1, π r (x r ) 0, if the exit is through, x l , π l (x l ) 0, π r (x r ) 1, if the exit is through, x r , π l (x 0 ) + π r (x 0 ) 1. The probabilities are given by The definite integrals in (25)

Mean Passage Times Through Boundaries with D < ∞
The solution to the following ordinary differential equation with the boundary conditions yields the mean first passage times through either x l or x r . The solution to (28) with integration constants k 1 , k 2 is

Passage Probabilities for SDE (7) with D ∞
The solution to the ordinary differential equation, with the same boundary conditions as in the last section, gives the exit probabilities π l (x 0 ) x r x 0 x which can then be integrated term by term. If 2u σ 2 is a large quantity however, on account of the order, u, being several orders of magnitude larger than σ 2 , the integrals in (31) behave like Laplace integrals [9] of the form,

Mean Passage Times Through Boundaries with D ∞
The differential equation is in this case with the boundary conditions

Stationary Solution of the Fokker-Planck Equation and Existence of Stable States
The Kolmogorov forward equation of Fokker-Planck equation governs the evolution of the transition probability density, f (x(t)|x 0 x(0)), henceforth denoted by f .

Stationary Solution to the Fokker-Planck Equation for SDE (6)
For (6) the Fokker-Planck equation reads If as t → ∞ the system attains a probability density, f * (x), independent of time, the system exhibits stationary behaviour. In this case (35) becomes an ordinary differential equation with solution where N is the integration (normalization) constant.
To qualify as a probability density, (36) must be normalizable, that is, Integrating (37) yields the definite integral The variable substitution, z 2u σ 2 x , in (38) leads to the following series representation: If u 0, it can be seen from (38) that N does not exist, and in this case, f * (x) δ(x), in [0, D). This is because the boundary, x 0, is an attracting boundary and the stationary probability mass will be concentrated entirely on zero.

Stationary Solution to the Fokker-Planck Equation for SDE (7)
The Fokker-Planck equation now reads and the stationary density reads where the normalization constant, N, is now furnished by the much simpler expression:

Extrema of Stationary Densities
The qualitative behaviour of the inventory process is determined by the extrema of the stationary density [10]. For (36) the extrema are supplied by the two roots of the quadratic equation: given by Both roots are real if u D ≤ (α+σ 2 ) 2 4α . For u D (α+σ 2 ) 2 4α , we have a double root, which is an inflection point at x * 1 x * 2 2u α+σ 2 . If the condition, u D < (α+σ 2 ) 2 4α , holds, the stable root, x * 1 , falls below D when either σ 2 ≤ α, or σ 2 > max α, u D . Differentiation of (43) with respect to x reveals that x * 1 is a relative maximum and x * 2 > x * 1 is a relative minimum. The inventory tends to move away from the relative minimum, x * 2 , towards the relative maximum, x * 1 , which represents the stable inventory state of the diffusion process.

Probabilistic Potentials
Use the exponent in (37) to define the following function: If u D < α 4 , then φ(x) possesses two extrema in [0,D) given by is uniformly decreasing so no global minimum exists, as the order rate, u, always exceeds demand growth (constantly positive drift).
The probabilistic potential, φ(x), is analogous to the potential (Lyapunov) function in Classical Mechanics. The root ξ 1 is a stable minimum and the second root, ξ 2 > ξ 1 , is an unstable maximum. The inventory will tend to drift towards levels that minimize φ(x) and maximize f * (x). But the maximum, x * 1 , of f * (x) does not in general coincide with the minimum, ξ 1 , of φ(x), unless the diffusion term is just an additive constant, independent of x. In practice, stable inventory values will be those that fall within the valley of φ(x) and the peak of f * (x). If φ(x) does not have a minimum but f * (x) still has a maximum, then there is a non-negligible probability that the inventory will fall anywhere in the range, [0, D), which is clearly an undesirable consequence when D is large in relation to the existing stock.
The density (41) has a unique maximum at and the associated probabilistic potential function, φ(x) 2u σ 2 x + 2α ln x σ 2 , a unique minimum at

Approximation of the Stationary Density by a Normal Density in the Vicinity of Its Extremum
The probability density function (36), We can approximate f * (x) by a Gaussian density, g(x), in the neighbourhood of x * 1 [10]. The Gaussian density must have the form so that g (x * 1 ) 0 and g (x * 1 ) −c f * (x * 1 ) ( f * ) (x * 1 ). The area, A, under the Gaussian function, g(x), is given by where erf(x) is the well known error function, erf(x) 2 √ π x 0 e −t 2 dt [11]. The effective width, ε, of the peak of the Gaussian function, g(x), is the width of the rectangle that has the same height as its peak, f * (x * 1 ), and the same area, A. So The probability density (41), can also be approximated by the normal density The effective width for (51) is similarly given by The effective width is useful in practice as it designates the range of approximately stable inventory levels located within ± ε 2 the theoretically obtained stable steady states, x * 1 (D < ∞) and x * (D ∞).
The Gaussian density, g(x), is a credible approximation to f * (x) when both φ(x) and f * (x) possess extreme values in [0,D), so that the effective width for g(x) represents the basin of stability.

A Numerical Example
We close the paper by a simple numerical demonstration of the key findings. Let x 0 5, u 5, D 400, α 0.2, σ 0.3. Figure 1 below shows 5 sample path realizations and the evolution of the mean inventory level. Figure 2 below displays the stationary probability density, f * (x), and its Gaussian approximation, g(x). The maximum density inventory is x * 1 ≈ 18 from (44) and the effective width is ε ≈ 18 from (50), hence the range of stable inventory values is approximately [9,27]. Finally, Fig. 3 illustrates the probabilistic potential, φ(x), maximized at ξ 1 ≈ 27 from (46).
Finally, suppose the inventory planner wants to estimate for instance, the probability that the inventory, starting from x 0 5 will either double to x r 10 or drop to x l 4, when the replenishment rate is u 5 and α 0.2, σ 0.3, D 400. From (25), (26) and (27) we obtain the probability estimates, π(x l 4) 0.0177 and π(x r 10) 0.9823. If the replenishment rate drops to u 2 for instance, the probabilities become roughly equal, π(x l 4) 0.5089 and π(x r 10) 0.4911.

Discussion
We have presented in this work a continuous time mathematical model for a randomly evolving inventory of an item for which demand is growing at a gradually slowing rate in relation to the inventory's availability, whilst being simultaneously satisfied at a constant order rate. The model put forward is a temporally homogeneous stochastic differential equation described by (6) and in a simpler reduced version by (7), with Itô solutions provided when analytically possible. To assess the direction the inventory is likely to take, theoretical estimates of the probabilities of attaining arbitrary inventory levels from some current inventory state under a fixed replenishment scheme are explicitly determined in "First Passage Times and Probabilities of Exit Through Absorbing Barriers" section. Finally, in "Stationary Solution of the Fokker-Planck Equation and Existence of Stable States" section the issue of long term stable stock levels is thoroughly addressed and the constraint on the replenishment rate, u, for a stable inventory regime is explicitly obtained. The size of the stable regime derived from the approximation of the probability density to a Gaussian density, will depend on the magnitude of the diffusion parameter, σ . The probabilities of reaching prescribed inventory levels and the determination of the stability regime are useful practical parameters for the inventory planner.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.