Stochastic Square of the Brennan-Schwartz Diffusion Process: Statistical Computation and Application

Abstract

In this paper, we study a new one-dimensional homogeneous stochastic process, termed the Square of the Brennan-Schwartz model, which is used in various contexts. We first establish the probabilistic characteristics of the model, such as the analytical expression solution to Itô’s stochastic differential equation, after which we determine the trend functions (conditional and non-conditional) and the likelihood approach in order to estimate the parameters in the drift. Then, in the diffusion coefficient, we consider the problem of parameter estimation, doing so by a numerical approximation. Finally, we present an application to population growth by the use of real data, namely the growth of the total population aged 65 and over, resident in the Arab Maghreb, to illustrate the research methodology presented.

Appendices

Appendix A: Ergodicity and stationary distribution of the SBSDP

Here, we study the asymptotic behaviour of the process proposed in this paper, analyse the problem of ergodicity and the existence of the stationary distribution of the process and explicitly obtain its density function.

In general (see Nobile and Ricciardi (1984) and Nicolau (2005)), a diffusion process {x(t), t ≥ 0}, with state space I = (l, r) , is governed by the following SDE:

$$ \begin{array}{@{}rcl@{}} dx(t)=a(x(t))dt+b(x(t))dW_{t}, \quad x_{0}=x, \end{array} $$

where W_t is a standard Wiener process and x is either a constant value or a random value independent of W_t. We assume that a(x) and b(x) have continuous derivatives.

Let $\displaystyle { s(z)=\exp \left \{-{\int }_{z_{0}}^{z}\frac {2a(u)}{b^{2}(u)}du\right \}}$ be the scale density function (z₀ is an arbitrary point inside I).

The speed density function is: m(u) = (b²(u)s(u))^− 1. And we denote by:

$$ S[x, y]= {{\int}_{x}^{y}}s(u)du, \quad S(l, y]= \lim_{x\rightarrow l}{{\int}_{x}^{y}}s(u)du \quad \text{and} $$

$$ S[x, r)= \lim_{y\rightarrow r}{{\int}_{x}^{y}}s(u)du,$$

where, l < x < y < r, and then if:

$$ \begin{array}{@{}rcl@{}} S(l, x]=S[x, r)=\infty \quad \text{and} \quad {{\int}_{l}^{r}}m(u)du<\infty, \end{array} $$

the process {x(t), t ≥ 0} is ergodic and has an invariant (stationary) density function which is given by:

$$ f(x)=m(x)/{{\int}_{l}^{r}}m(u)du. $$

In our diffusion, the drift and diffusion coefficient are, respectively:

$$a(x)=\alpha x+\beta \sqrt{x} \quad \text{and} \quad b^{2}(x)=\sigma^{2} x^{2},$$

and I = (0, ∞). In this case, it follows that:

$$s(z)=k z^{-\frac{2\alpha}{\sigma^{2}}}e^{\frac{4\beta}{\sigma^{2}\sqrt{z}}}, \quad \text{with} \quad k=z_{0}^{\frac{2\alpha}{\sigma^{2}}}e^{-\frac{4\beta}{\sigma^{2}\sqrt{z_{0}}}}.$$

and we have, for 0 < x < y < ∞

$$ \begin{array}{@{}rcl@{}} S[x, y]={{\int}_{x}^{y}}s(u)du=k{{\int}_{x}^{y}}u^{-\frac{2\alpha}{\sigma^{2}}}e^{\frac{4\beta}{\sigma^{2}\sqrt{u}}}du, \end{array} $$

with the variable change v = u^− 1/2, the latter expression can be written as

$$ \begin{array}{@{}rcl@{}} S[x, y]= 2k{\int}_{1/\sqrt{y}}^{1/\sqrt{x}}v^{\frac{4\alpha}{\sigma^{2}}-3}e^{\frac{4\beta} {\sigma^{2}}v}dv. \end{array} $$

(20)

On the one hand, taking the limit as x tends to 0 in (20). We have, for β > 0, S(0, y] = ∞ (the case of β = 0 is excluded, because the process is lognormal, and this is not ergodic). On the other hand, taking the limit when y tends to ∞ in (20), we have, for $\alpha \leq \displaystyle {\frac {\sigma ^{2}}{2}}$, S[x, ∞) = ∞. And therefore, for $\alpha \leq \displaystyle {\frac {\sigma ^{2}}{2}}$ and β > 0,

$$S[x, \infty)=S(0, y]=\infty.$$

The speed density in this case is

$$ \displaystyle{m(x)=\frac{1}{k\sigma^{2}}x^{\frac{2\alpha}{\sigma^{2}}-2}e^{\frac{-4\beta} {\sigma^{2} \sqrt{x}}}}, $$

and we have:

$$ \begin{array}{@{}rcl@{}} {\int}_{0}^{\infty}m(x)dx&=&\frac{1}{k\sigma^{2}}{\int}_{0}^{\infty}x^{\frac{2\alpha}{\sigma^{2}}-2}e^{\frac{-4\beta} {\sigma^{2} \sqrt{x}}}dx\\ &=&\frac{2}{k\sigma^{2}}{\int}_{0}^{\infty}v^{\frac{-4\alpha}{\sigma^{2}}+1}e^{-\frac{4\beta} {\sigma^{2}} v}dv, \end{array} $$

and according to Gradshteyn and Ryzhik (1979) 3.18. p.317, for ν > 0 and μ > 0,

$$ {\int}_{0}^{\infty}x^{\nu-1}e^{-\mu x}dx=\mu^{-\nu}{\Gamma}(\nu),$$

we have, for $\alpha < \displaystyle {\frac {\sigma ^{2}}{2}}$ and β > 0,

$$ \begin{array}{@{}rcl@{}} \displaystyle{{\int}_{0}^{\infty}m(x)dx=\frac{2}{k\sigma^{2}} \left( \frac{4\beta}{\sigma^{2}}\right)^{\left( \frac{4\alpha}{\sigma^{2}}-2\right)} {\Gamma}\left( 2-\frac{4\alpha}{\sigma^{2}}\right)<\infty}, \end{array} $$

Then, by combining the two conditions, we deduce that for $\alpha < \displaystyle {\frac {\sigma ^{2}}{2}}$ and β > 0, the process is ergodic and its stationary density function is given by the following expression:

$$ \begin{array}{@{}rcl@{}} \displaystyle{f(x)=m(x)/{\int}_{0}^{\infty}m(u)du=\frac{\mu^{\lambda}}{2{\Gamma}\left( \lambda\right)} x^{-\left( \frac{\lambda}{2}+1\right)}e^{-\frac{\mu}{\sqrt{x}}}} \end{array} $$

where $\lambda =2-\displaystyle {\frac {4\alpha }{\sigma ^{2}}}$ and $\mu =\displaystyle {\frac {4\beta }{\sigma ^{2}}}$, and the conditions of ergodicity in terms of λ and μ are equivalent to λ > 0 and μ > 0.

Appendix B: Approximate estimator of the diffusion coefficient of the SBSDP

To estimate the parameter σ, we used an extension of the method described by Chesney and Elliott (1995). This method has been used by Giovanis and Skiadas (1999) in a paper on a stochastic logistic innovation diffusion model studying electricity consumption in Greece and the United States, and in another on a stochastic Bass innovation diffusion model to study the growth of electricity consumption in Greece, see Skiadas and Giovanis (1997).

In our study, an approximate estimator of the σ parameter between two observations has the general form:

$$ \begin{array}{@{}rcl@{}} \hat{\sigma}_{(t-1,t)}=\frac{\mid x(t)-x(t-1)\mid}{\sqrt{x(t) x(t-1)} }, \end{array} $$

Then, we have

$$ \begin{array}{@{}rcl@{}} \hat{\sigma}_{(0,1)}=\frac{\mid x(1)-x(0)\mid}{\sqrt{x(1) x(0)} }, \end{array} $$

$$ \begin{array}{@{}rcl@{}} \hat{\sigma}_{(1,2)}=\frac{\mid x(2)-x(1)\mid}{\sqrt{x(2) x(1)} }, \end{array} $$

$$ \begin{array}{@{}rcl@{}} \vdots \end{array} $$

$$ \begin{array}{@{}rcl@{}} \hat{\sigma}_{(n-1,n)}=\frac{\mid x(n)-x(n-1)\mid}{\sqrt{x(n) x(n-1)} }. \end{array} $$

The method described by Chesney and Elliot was then applied to each time interval, and the average of these estimators for n observations of a sample path of the process is the approximate estimator. It takes the following form:

$$ \begin{array}{@{}rcl@{}} \hat{\sigma}=\frac{1}{n-1}\sum\limits_{t=1}^{n}\frac{\mid x(t)-x(t-1)\mid}{x(t)\sqrt{x(t-1)} }. \end{array} $$

Other approximate estimators of σ can be obtained by the same procedure, for example:

• Using Itô’s lemma to the transformationy = ln(x(t)) in (1) as follows:

$$ \begin{array}{@{}rcl@{}} \displaystyle{d(lnx(t))=\frac{dx(t)}{x(t)} -\frac{\sigma^{2}}{2} dt}, \end{array} $$

By substituting:

$$ (d(ln x(t)))^{2}= \sigma^{2} dt,$$

Considering that

$$\left( d(ln x(t))\right)\simeq ln(x(t))- ln(x(t-1)),$$

an approximate estimator of σ is:

$$ \begin{array}{@{}rcl@{}} \hat{\sigma}=\mid ln(x(t))-ln(x(t-1))\mid, \end{array} $$

Then, for n observations of a sample path of the process, the resulting approximate estimator has the following expression:

$$ \begin{array}{@{}rcl@{}} \hat{\sigma}=\frac{1}{n-1}\sum\limits_{t=1}^{n}\mid ln(x(t))-ln(x(t-1))\mid. \end{array} $$

• Alternatively, from the stochastic differential equation (1) of the variablex(t), we obtain:

$$ \begin{array}{@{}rcl@{}} \displaystyle{\left( \frac{dx(t)}{x(t)}\right)= \left( \alpha+ \frac{\beta}{\sqrt{x(t)}}\right) dt+ \sigma dw(t)}, \end{array} $$

and then:

$$ \begin{array}{@{}rcl@{}} \displaystyle{\left( \frac{dx(t)}{x(t)}\right)^{2}= \sigma^{2} dt}, \end{array} $$

Considering that

$$d(x(t))\simeq x(t)-x(t-1),$$

An approximate value of σ is:

$$ \begin{array}{@{}rcl@{}} \hat{\sigma}= \frac{\mid x(t)-x(t-1)\mid}{x(t)} , \end{array} $$

Therefore, for n observations of a sample path of the process, the resulting approximate estimator has the following expression:

$$ \begin{array}{@{}rcl@{}} \hat{\sigma}=\frac{1}{n-1}\sum\limits_{t=1}^{n}\frac{\mid x(t)-x(t-1)\mid}{x(t)}. \end{array} $$

By applying the three expressions of $\hat {\sigma }$ to our real data, we obtain the following results:

Table 9 Values obtained using the three expressions of $\hat {\sigma }$

Similar values are obtained by all three methods.

Appendix C: Probabilistic characteristics and statistical inference of the SLDP

The SLDP from the SBSDP

In the equation (1) when β = 0, the homogeneous lognormal diffusion process is obtained as a particular case in which the infinitesimal moments are given by:

$$ \begin{array}{@{}rcl@{}} A_{1}(x)=\alpha x \hspace{0.5cm} , \hspace{0.5cm} A_{2}(x)=\sigma^{2} x^{2},\hspace{0.5cm} \end{array} $$

This satisfies the following Itô’s SDE:

$$ dx(t)=\alpha x(t)dt+ \sigma x(t) dw(t) \hspace{0.5cm},\hspace{0.5cm} x(t_{0})=x_{t_{0}}\hspace{0.5cm}, $$

where σ > 0 and α are real parameters, W_t is a standard Wiener process and $x_{t_{0}}$ is fixed in $\mathbb {R}^{*}_{+}$.

• The analytical expression of the SLDP

  By taking β = 0 in equation (3) of the Section 2.2, the previous SDE has a unique solution which is given analytically by the following expression:

  $$ \begin{array}{@{}rcl@{}} x(t)&=&x_{t_{0}} \left( \exp\left[\left( \alpha-\frac{\sigma^{2}}{2}\right) \left( t-t_{0}\right)+\sigma\left( w(t)-w(t_{0})\right)\right]\right) \end{array} $$

• The trend functions of the SLDP

  In the same way, by taking β = 0 in equations (4) and (5 of the Section 2.3, the conditional trend function of the process is:

  $$ \begin{array}{@{}rcl@{}} \mathbb{E}(x(t)\mid x(s)=x_{s})=x_{s} e^{\alpha(t-s)}. \end{array} $$

  (21)

  and by assuming the initial condition $\mathrm {P}\left (x(t_{0})=x_{t_{0}}\right )=1$, the trend function of the process is:

  $$ \begin{array}{@{}rcl@{}} \mathbb{E}(x(t))=x_{t_{0}} e^{\alpha(t-t_{0})}. \end{array} $$

  (22)

Parameter estimation in the SLDP

We now determine the estimator of the parameter α of the SLDP using the method described. The estimator of the drift parameters α is obtained by the maximum likelihood method, with continuous sampling.

• Likelihood estimation of the drift parameter:

  The vector form of the SDE of the SLDP can be written as:

  $$ \begin{array}{@{}rcl@{}} A_{t}(x(t))= x(t), \quad \theta^{*}= \alpha, \quad B_{t}(x(t))= \sigma x(t), \end{array} $$

  The corresponding vector H_T in this case is one-dimensional and is given by:

  $$ \begin{array}{@{}rcl@{}} H^{*}_{T}= \displaystyle{\frac{1}{\sigma^{2}}} {\int}_{t_{0}}^{T}\frac{dx(t)}{x(t)}, \end{array} $$

  and S_T has the following form:

  $$ \begin{array}{@{}rcl@{}} S_{T}= \frac{T-t_{0}}{\sigma^{2}} \end{array} $$

  Then the expression of the estimator is

  $$ \begin{array}{@{}rcl@{}} \displaystyle{\hat{\alpha}=\frac{ {\int}_{t_{0}}^{T}\frac{dx(t)}{x(t)}} {T-t_{0}}}, \end{array} $$

  The stochastic integral in the latter expression can be transformed into Riemann integrals by using Itô’s formula and thus:

  $$ \begin{array}{@{}rcl@{}} {\int}_{t_{0}}^{T}\frac{dx(t)}{x(t)}&=& \log(x_{T})-\log(x_{t_{0}}) +\frac{\sigma^{2}}{2}\left( T-t_{0}\right). \end{array} $$

  Therefore, the expression of the Maximum Likelihood estimator is:

  $$ \displaystyle{\hat{\alpha}=\frac{\left( \log(x_{T})- \log(x_{t_{0}}\right) +\frac{\sigma^{2}}{2}\left( T-t_{0}\right)} {T-t_{0}}}. $$

• Approximate estimator of the diffusion coefficient

The estimator of the coefficient diffusion parameter can be approximated using a method similar to that described Section 3.2. By following the same steps and for n observations of a sample path of the process, the resulting approximate estimator is:

$$ \begin{array}{@{}rcl@{}} \hat{\sigma}=\frac{1}{n-1}\sum\limits_{t=1}^{n}\frac{\mid x(t)-x(t-1)\mid}{\sqrt{ x(t) x(t-1)} }. \end{array} $$

Computational aspects in SLDP

• Approximated likelihood estimators

  As for the SBSDP, in order to use the above expression to estimate the parameter α, we must have continuous observations. Therefore, we use an alternative estimation procedure based on continuous time maximum likelihood estimators with suitable approximations of the integrals that appear in the expression. The Riemann-Stieljes integrals are approximated by means of the trapezoidal formula.

• Estimated trend functions

  By applying Zehna’s theorem and by taking β = 0 in the equations (21) and (22), the estimated trend function (ETF) and estimated conditional trend function (ECTF) of the process are obtained as:

  $$ \begin{array}{@{}rcl@{}} \widehat{\mathbb{E}}(x(t)/ x(s)=x_{s})=x_{s} e^{\hat{\alpha}(t-s)}, \end{array} $$
  $$ \begin{array}{@{}rcl@{}} \widehat{\mathbb{E}}(x(t))=x_{t_{0}} e^{\hat{\alpha}(t-t_{0})}. \end{array} $$