Global Hopf bifurcation of a delayed phytoplankton-zooplankton system considering toxin producing effect and delay dependent coefficient

1.   Introduction

For a continuous risk outcome $0<y<1$, a model with a random effect has potentially a wide application in portfolio risk management, especially, for stress testing ^{[1,2,7,16,19]}, capital allocation, conditional expected shortfall estimation ^[3,11,17].

Given fixed effects ${x = (x}_{1}, {x}_{2}, \dots , {x}_{k})$, two widely used regression models to estimate the expected value $E\left(y\right|x)$ are: the fraction response model ^[10] and Beta regression model ^[4,6,8]. There are cases, however, where tail behaviours or severity levels of the risk outcome are relevant. In those cases, a regression model may no longer fit in for the requirements. In addition, a fraction response model of the form $E\left(y|x\right) = \mathrm{\Phi }\left({a}_{0}+{a}_{1}{x}_{1}+\dots +{a}_{k}{x}_{k}\right)$ may not be adequate when data exhibits significant heteroscedasticity, where $\mathrm{\Phi }$ is a map from ${R}^{1}$ to the open interval $\left(0, 1\right).$

In this paper, we assume that the risk outcome $y$ is driven by a model:

where $s$ is a random continuous variable following a known distribution, independent of fixed effects ${(x}_{1}, {x}_{2}, \dots , {x}_{k})$. Parameters ${a}_{0}, {a}_{1}, \dots , {a}_{k}$ are constant, while parameter $b$ can be chosen to be dependent on ${(x}_{1}, {x}_{2}, \dots , {x}_{k})$ when required, for example, for addressing data heteroscedasticity.

Given random effect model (1.1), the expected value $E\left(y\right|x)$ can be deduced accordingly. It is given by the integral ${\int }_{\mathit{\Omega} }\mathrm{\Phi }\left({a}_{0}+{a}_{1}{x}_{1}+\dots +{a}_{k}{x}_{k}+bs\right)f\left(s\right)ds$ over the domain $\mathit{\Omega}$ of $s$, where $f$ is the probability density of $s$. Given the routine QUAD implemented in SAS and Python, this integral can be evaluated as quickly as other function calls. Relative error tolerance for QUAD is 1.49e-8 in Python and is 1e-7 in SAS. But one can rescale the default tolerance to a desired level when necessary. This leads to an alternative regression tool to the fraction response model and Beta regression model.

We introduce a family of interval distributions based on variable transformations. Probability densities for these distributions are provided (Proposition 2.1). Parameters of model (1.1) can then be estimated by maximum likelihood approaches assuming an interval distribution. In some cases, these parameters get an analytical solution without the needs for a model fitting (Proposition 4.1). We call a model with a random effect, where parameters are estimated by maximum likelihood assuming an interval distribution, an interval distribution model.

In its simplest form, the interval distribution model $y = \mathrm{\Phi }\left(a+bs\right)$, where $a$ and $b,$ are constant, can be used to model the loss rate as a random distribution for a homogeneous portfolio. Let ${y}_{\alpha }$ and ${s}_{\alpha }$ denote the $\alpha$ -quantiles for $y$ and $s$ at level $\alpha$, $0<\alpha <1.$ Then ${y}_{\alpha } = \mathrm{\Phi }\left(a+b{s}_{\alpha }\right).$ The conditional expected shortfall for loss rate $y$, at level $\alpha$, can then be estimated as the integral $\frac{1}{1-\alpha }{\int }_{[{s}_{\alpha }, +\infty )}\mathrm{\Phi }\left(a+bs\right)f\left(s\right)ds$, where $f$ is the density of $s.$ Meanwhile, a stress testing loss estimate, derived from a model on a specific scenario, can be compared in loss rate to severity ${y}_{\alpha }( = \mathrm{\Phi }\left(a+b{s}_{\alpha }\right))$, to position its level of severity. A loss estimate may not have reached the desired, for example, 99% level yet, if it is far below ${y}_{0.99}$, and far below the maximum historical loss rate. In which case, further recalibrations for the model may be required.

The paper is organized as follows: in section 2, we introduce a family of interval distributions. A measure for tail fatness is defined. In section 3, we show examples of interval distributions and investigate their tail behaviours. We propose in section 4 an algorithm for estimating the parameters in model (1.1).

2.   Interval distributions generated by transformations

Interval distributions introduced in this section are defined for a risk outcome over a finite open interval $\left({c}_{0}, {c}_{1}\right),$ where ${c}_{0}<$ ${c}_{1}$ are finite numbers. These interval distributions can potentially be used for modeling a risk outcome over an arbitrary finite interval, including interval (0, 1), by maximum likelihood approaches.

Let $D = ({d}_{0}, {d}_{1})$, ${d}_{0}< {d}_{1}$, be an open interval, where ${d}_{0}$ can be finite or $-\infty$ and ${d}_{1}$ can be finite or $+\infty$.

Let

be a transformation with continuous and positive derivatives ${\mathrm{\Phi }}^{\mathrm{'}}\left(\mathrm{x}\right) = {\rm{ \mathsf{ ϕ} }}\left(x\right)$. A special example is $\left({c}_{0}, {c}_{1}\right) = \left(0, 1\right)$, and $\mathrm{\Phi }: D\to (0, 1)$ is the cumulative distribution function (CDF) of a random variable with a continuous and positive density.

Given a continuous random variable $s,$ let $f$ and $F$ be respectively its density and CDF. For constants $a$ and $b>0$, let

where we assume that the range of variable $(a+bs)$ is in the domain $D$ of $\mathrm{\Phi }.$ Let $g\left(y, a, b\right)$ and $G\left(y, a, b\right)$ denote respectively the density and CDF of $y$ in (2.2).

Proposition 2.1. Given ${\mathrm{\Phi }}^{-1}\left(y\right),$ functions $g\left(y, a, b\right)$ and $G\left(y, a, b\right)$ are given as:

where

Proof. A proof for the case when $\left({c}_{0}, {c}_{1}\right) = \left(0, 1\right)$ can be found in ^[18]. The proof here is similar. Since $G\left(y, a, b\right)$ is the CDF of $y$, it follows:

By chain rule and the relationship $\mathrm{\Phi }\left[{\mathrm{\Phi }}^{-1}\left(y\right)\right] = y,$ the derivative of ${\mathrm{\Phi }}^{-1}\left(y\right)$ with respect to $y$ is

Taking the derivative of $G\left(y, a, b\right)$ with respect to $y$, we have

One can explore into these interval distributions for their shapes, including skewness and modality. For stress testing purposes, we are more interested in tail risk behaviours for these distributions.

Recall that, for a variable X over (− $\infty , +\infty ),$ we say that the distribution of X has a fat right tail if there is a positive exponent $\alpha >0$, called tailed index, such that $P\left(X>x\right) \sim {x}^{-\alpha }.$ The relation $\sim$ refers to the asymptotic equivalence of functions, meaning that their ratio tends to a positive constant. Note that, when the density is a continuous function, it tends to 0 when $x\to +\infty .$ Hence, by L’Hospital’s rule, the existence of tailed index is equivalent to saying that the density decays like a power law, whenever the density is a continuous function.

For a risk outcome over a finite interval $\left({c}_{0}, {c}_{1}\right)$, ${c}_{0}, < {c}_{1},$ however, its density can be $+\infty$ when approaching boundaries ${c}_{0}$ and ${c}_{1.}$ Let ${y}_{0}$ be the largest lower bound for all values of $y$ under (2.2), and ${y}_{1}$ the smallest upper bound. We assume ${y}_{0} = {c}_{0}$ and ${y}_{1} = {c}_{1}.$

We say that an interval distribution has a fat right tail if the limit $\mathit{li}{m}_{y⤍{y}_{1}^{-}} \ \ g \left(y, a, b\right) = +\infty ,$ and a fat left tail if $\mathit{li}{m}_{y⤍{y}_{0}^{+}} \ \ g \left(y, a, b\right) = +\infty$, where $y⤍{y}_{0}^{+}$ and $y⤍{y}_{1}^{-}$ denote respectively $y$ approaching ${y}_{0}$ from the right-hand-side, and ${y}_{1}$ from the left-hand-side. For simplicity, we write $y⤍{y}_{0}$ for $y⤍{y}_{0}^{+},$ and $y⤍{y}_{1}$ for $y⤍{y}_{1}^{-}.$

Given $\alpha >0,$ we say that an interval distribution has a fat right tail with tailed index $\alpha$ if $\mathit{li}{m}_{y⤍{y}_{1}} \ \ g \left(y, a, b\right){\left({y}_{1}-y\right)}^{\beta } = +\infty$ whenever $0<\beta <\alpha ,$ and $li{m}_{y⤍{y}_{1}} \ \ g \left(y, a, b\right){\left({y}_{1}-y\right)}^{\beta } = 0$ for $\beta >\alpha .$ Similarly, an interval distribution has a fat left tail with tailed index $\alpha$ if $li{m}_{y⤍{y}_{0}} \ \ g \left(y, a, b\right){\left({y-y}_{0}\right)}^{\beta } = +\infty$ whenever $0<\beta <\alpha ,$ and $li{m}_{y⤍{y}_{0}} \ \ g \left(y, a, b\right)({y-{y}_{0})}^{\beta } = 0$ for $\beta >\alpha .$ Here the status at $\beta = \alpha$ is left open. There are examples (Remark 3.4), where an interval distribution has a fat right tail with tailed index $\alpha$, but the limit $li{m}_{y⤍{y}_{1}} \ \ g \left(y, a, b\right){\left({y}_{1}-y\right)}^{\alpha }$ can either be $+\infty$ or 0. Under this definition, a tailed index of an interval distribution with a continuous density is always larger than 0 and less or equal to 1, if it exists.

Recall that, for a Beta distribution with parameters $\alpha >0$ and $\beta >0,$ its density is given by $f\left(x\right) = \frac{{x}^{\alpha -1}{\left(1-x\right)}^{\beta -1}}{B(\alpha , \beta )},$ where $B(\alpha , \beta )$ is the Beta function $.$ Under the above definition, Beta distribution has a fat right tail with tailed index $(1-\beta )$ when $0<\beta <1$, and a fat left tail with tailed index $(1-\alpha )$ when $0<\alpha <1.$

Next, because the derivative of $\mathrm{\Phi }$ is assumed to be continuous and positive, it is strictly monotonic. Hence ${\mathrm{\Phi }}^{-1}\left(y\right)$ is defined. Let

Then $li{m}_{y⤍{y}_{0}}z$ exists (can be $-\infty )$, and the same for $li{m}_{y⤍{y}_{1}}z$ (can be $+\infty )$. Let $li{m}_{y⤍{y}_{0}} \ \ z = {z}_{0}$, and $li{m}_{y⤍{y}_{1}} \ \ z = {z}_{1}$. Rewrite $g\left(y, a, b\right)$ as $g\left(\mathrm{\Phi }\left(z\right), a, b\right)$ by (2.7). Let ${\partial \left[g\left(\mathrm{\Phi }\left(\mathrm{z}\right), a, b\right)\right]}^{-\frac{1}{\beta }}/\partial z$ denote the derivative of ${\left[g\left(\mathrm{\Phi }\left(\mathrm{z}\right), a, b\right)\right]}^{-1/\beta }$ with respect to $z$.

Lemma 2.2. Given $\beta >0$, the following statements hold:

(ⅰ) $\mathit{li}{m}_{y⤍{y}_{0}} \ \ g \left(y, a, b\right){\left({y-y}_{0}\right)}^{\beta } = \mathit{li}{m}_{z⤍{z}_{0}} \ \ g \left(\mathrm{\Phi }\left(\mathrm{z}\right), a, b\right){\left(\mathrm{\Phi }\left(\mathrm{z}\right)-{\mathrm{y}}_{0}\right)}^{\beta }$ and $li{m}_{y⤍{y}_{1}} \ \ g \left(y, a, b\right){\left({y}_{1}-y\right)}^{\beta } = li{m}_{z⤍{z}_{1}} \ \ g \left(\mathrm{\Phi }\left(z\right), a, b\right){\left({y}_{1}-\mathrm{\Phi }\left(z\right)\right)}^{\beta }$.

(ⅱ) If $li{m}_{y⤍{y}_{0}} \ \ g \left(y, a, b\right) = +\infty$ and $li{m}_{z⤍{z}_{0}}{\{\partial \left[g\left(\mathrm{\Phi }\left(\mathrm{z}\right), a, b\right)\right]}^{-\frac{1}{\beta }}/\partial z\}/{\rm{ \mathsf{ ϕ} }}(\mathrm{z})$ is 0 (resp. $+\infty )$, then $li{m}_{y⤍{y}_{0}} \ \ g \left(y, a, b\right){(y-{y}_{0})}^{\beta } = +\infty$ (resp. 0).

(ⅲ) If $li{m}_{y⤍{y}_{1}} \ \ g \left(y, a, b\right) = +\infty$ and $li{m}_{z⤍{z}_{1}}-{\{\partial \left[g\left(\mathrm{\Phi }\left(\mathrm{z}\right), a, b\right)\right]}^{-\frac{1}{\beta }}/\partial z\}/{\rm{ \mathsf{ ϕ} }}(\mathrm{z}\left)\right)$ is 0 (resp. $+\infty )$, then $li{m}_{y⤍{y}_{1}} \ \ g \left(y, a, b\right){\left({y}_{1}-y\right)}^{\beta } = +\infty$ (resp. 0).

Proof. The first statement follows from the relationship $y = \mathrm{\Phi }(\mathrm{z}$). For statements (ⅱ) and (ⅲ), we show only (ⅲ). The proof for (ⅱ) is similar. Notice that

By L’Hospital’s rule and taking the derivatives of the numerator and the denominator of (2.8) with respect to $z$, we have ${li{m}_{y⤍{y}_{1}}\left[g\left(y, a, b\right){\left({y}_{1}-y\right)}^{\beta }\right]}^{-1/\beta } = 0$ (resp. $+\infty$) if $li{m}_{z⤍{z}_{0}}{-\{\partial \left[g\left(\mathrm{\Phi }\left(z\right), a, b\right)\right]}^{-1/\beta }/\partial z\}/{\rm{ \mathsf{ ϕ} }}(z)$ is 0 (resp. $+\infty )$. Hence $li{m}_{y⤍{y}_{1}} \ \ g \left(y, a, b\right){\left({y}_{1}-y\right)}^{\beta } = +\infty$ (resp. 0).

For tail convexity, we say that the right tail of an interval distribution is convex if $g\left(y, a, b\right)$ is convex for ${y}_{1}-є<y< {y}_{1}$ for sufficiently small $є>0.$ Similarly, the left tail is convex if $g\left(y, a, b\right)$ is convex for ${y}_{0}<y< {y}_{0}+є$ for sufficiently small $є>0.$ One sufficient condition for convexity for the right (resp. left) tail is ${g}_{yy}^{''}(y, a, b)\ge 0$ when $y$ is sufficiently close to ${y}_{1}$ (resp. ${y}_{0})$.

Again, write $g\left(y, a, b\right) = g\left(\mathrm{\Phi }\left(z\right), a, b\right)$. Let

where $\mathrm{l}\mathrm{o}\mathrm{g}\left(x\right)$ denotes the natural logarithmic function. Then

By (2.9), (2.10), using (2.6) and the relationship ${z = \mathrm{\Phi }}^{-1}\left(y\right)$, we have

The following lemma is useful for checking tail convexity, it follows from (2.11).

Lemma 2.3. Suppose ${\rm{ \mathsf{ ϕ} }}\left(\mathrm{z}\right)>0$, and derivatives ${h}_{\mathrm{z}}^{\mathrm{'}}\left(\mathrm{z}\right), {h}_{\mathrm{z}}^{\mathrm{'}\mathrm{'}}\left(\mathrm{z}\right),$ and ${{\rm{ \mathsf{ ϕ} }}}_{\mathrm{z}}^{\mathrm{'}}\left(\mathrm{z}\right)$, with respect to $z,$ all exist. If ${h}_{zz}^{''}\left(z\right)\ge 0$ and ${h}_{z}^{'}\left(z\right){{\rm{ \mathsf{ ϕ} }}}_{\mathrm{z}}^{'}\left(z\right)\le 0,$ then ${g}_{yy}^{''}(y, a, b)\ge 0.$

3.   Examples of interval distributions and their tail behaviours

In this section, we focus on the case where $\left({c}_{0}, {c}_{1}\right) = \left(0, 1\right),$ and $\mathrm{\Phi }: D\to (0, 1)$ in (2.2) is the CDF of a continuous distribution $.$ This includes, for example, the CDFs for standard normal and standard logistic distributions.

One can explore into a wide list of densities with different choices for $\mathrm{\Phi }$ and $s$ under (2.2). We consider here only the following four interval distributions:

A. $s \sim N\left(\mathrm{0, 1}\right)$ and $\mathrm{\Phi }$ is the CDF for the standard normal distribution.

B. $s$ follows the standard logistic distribution and $\mathrm{\Phi }$ is the CDF for the standard normal distribution.

C. $s$ follows the standard logistic distribution and $\mathrm{\Phi }$ is its CDF.

D.D. $s \sim N\left(\mathrm{0, 1}\right)$ and $\mathrm{\Phi }$ is the CDF for standard logistic distribution.

Densities for cases A, B, C, and D are given respectively in (3.3) (section 3.1), (A.1), (A.3), and (A5) (Appendix A). Tail behaviour study is summarized in Propositions 3.3, 3.5, and Remark 3.6. Sketches of density plots are provided in Appendix B for distributions A, B, and C.

3.1.   Case A: the vasicek distribution and its tail behaviours

Using the notations of section 2, we have ${\rm{ \mathsf{ ϕ} }} = f$ and $\mathrm{\Phi } = F$. We claim that $y = \mathrm{\Phi }\left(a+bs\right)$ under (2.2) follows the Vasicek distribution ^[13,14].

By (2.5), we have

Therefore, we have

Again, using the notations of section 2, we have ${y}_{0} = 0$ and ${y}_{1} = 1$. With $z = {\mathrm{\Phi }}^{-1}\left(y\right),$ we have $li{m}_{y⤍0} \ \ z = -\infty$ and $li{m}_{y⤍1} \ \ z = +\infty .$ Recall that a variable $0<y<1$ follows a Vasicek distribution ^[13,14] if its density has the form:

where $p$ is the mean of y $,$ and $\rho$ is a parameter called asset correlation.

Proposition 3.1. Density (3.3) is equivalent to (3.4) under the relationships:

Proof. A similar proof can be found in ^[19]. By (3.4), we have

The following relationships are implied by (3.5):

Remark 3.2. The mode of $g\left(y, p, \rho \right)$ in (3.4) is given in ^[14] as $\mathrm{\Phi }\left(\frac{\sqrt{1-\rho }}{1-2\rho }{\mathrm{\Phi }}^{-1}\left(p\right)\right)$. We claim this is the same as $\mathrm{\Phi }\left(\frac{a}{1-{b}^{2}}\right)$. By (3.6), $1-2\rho = \frac{1-{b}^{2}}{1+{b}^{2}}$ and $\sqrt{1-\rho } = \frac{1}{\sqrt{1+{b}^{2}}}.$ Therefore, we have

This means $\mathrm{\Phi }\left(\frac{\sqrt{1-\rho }}{1-2\rho }{\mathrm{\Phi }}^{-1}\left(p\right)\right) = \mathrm{\Phi }\left(\frac{a}{1-{b}^{2}}\right).$

Proposition 3.3. The following statements hold for $g(y, a, b)$ given in (3.3):

(ⅰ) $g(y, a, b)$ is unimodal if $0<b<1$ with mode given by $\mathrm{\Phi }\left(\frac{a}{1-{b}^{2}}\right)$, and is in U-shape if $b>1$.

(ⅱ) $\mathrm{I}\mathrm{f}\ b>1, \mathrm{t}\mathrm{h}\mathrm{e}\mathrm{n} \ \ g(y, a, b)$ has a fat left tail and a fat right tail with tailed index $(1-1/{b}^{2})$.

(ⅲ) If $b>1,$ both tails of $g(y, a, b)$ are convex $,$ and is globally convex if in addition $a = 0.$

Proof. For statement (ⅰ), we have $-\left(1-{b}^{2}\right)<0$ when $0<b<1$. Therefore by (3.2) function $\mathrm{log}\left(\frac{{U}_{1}}{{U}_{2}}\right)$ reaches its unique maximum at $z = \frac{a}{1-{b}^{2}}$, resulting in a value for the mode at $\mathrm{\Phi }\left(\frac{a}{1-{b}^{2}}\right).$ If $b>1$, then $-\left(1-{b}^{2}\right)>0,$ thus by (3.2), $g(y, a, b)$ is first decreasing and then increasing when $y$ varying from 0 to 1. This means $(y, a, b$) is in U-shape.

Consider statement (ⅱ). First by (3.3), if $b>1,$ then $li{m}_{y⤍1}\ \ g\left(y, a, b\right) = +\infty$ and $li{m}_{y⤍0} \ \ g \left(y, a, b\right) = +\infty$. Thus $g\left(\mathrm{y}, a, b\right)$ has a fat right and a fat left tail. Next for tailed index, we use Lemma 2.2 (ⅱ) and (ⅲ). By (3.1),

By taking the derivative of (3.8) with respect to $z$ and noting that ${\rm{ \mathsf{ ϕ} }}\left(\mathrm{z}\right) = \frac{1}{\sqrt{2\pi }}\mathrm{exp}\left(-\frac{{z}^{2}}{2}\right),$ we have

Thus $li{m}_{z⤍+\infty }-\left\{\partial {\left[g\left(\mathrm{\Phi }\left(\mathrm{z}\right), a, b\right)\right]}^{-\frac{1}{\beta }}/\partial z\right\}/{\rm{ \mathsf{ ϕ} }}\left(\mathrm{z}\right)$ is 0 if $\frac{{b}^{2}-1}{\beta {b}^{2}}>1$, and is $+\infty$ if $\frac{{b}^{2}-1}{\beta {b}^{2}}<1.$ Hence by Lemma 2.2 (ⅲ), $g\left(y, a, b\right)$ has a fat right tail with tailed index $(1-1/{b}^{2})$. Similarly, for the left tail, we have by (3.9)

Thus $li{m}_{z⤍-\infty }\left\{\partial {\left[g\left(\mathrm{\Phi }\left(\mathrm{z}\right), a, b\right)\right]}^{-\frac{1}{\beta }}/\partial z\right\}/{\rm{ \mathsf{ ϕ} }}\left(\mathrm{z}\right)$ is 0 if $\frac{{b}^{2}-1}{\beta {b}^{2}}>1$, and is $+\infty$ if $\frac{{b}^{2}-1}{\beta {b}^{2}}<1.$ Hence $g\left(y, a, b\right)$ has a fat left tail with tailed index $(1-1/{b}^{2})$ by Lemma 2.2 (ⅱ).

For statement (ⅲ), we use Lemma 2.3. By (2.9) and using (3.2), we have

When $b>1,$ it is not difficult to check out that ${h}_{zz}^{''}\left(z\right)\ge 0$ and ${h}_{z}^{'}\left(z\right){{\rm{ \mathsf{ ϕ} }}}_{\mathrm{z}}^{'}\left(z\right)\le 0$ when $z⤍\pm \infty$ or when $a = 0$.

Remark 3.4. Assume $\beta = (1-1/{b}^{2})$ and $b>1.$ By (3.9), we see

is $+\infty$ for $a = 0$, and is 0 for $a>0.$ Hence for this $\beta ,$ the limit $li{m}_{y⤍1} \ \ g \left(y, a, b\right){\left(1-y\right)}^{\beta }$ can be either 0 or $+\infty ,$ depending on the value of $a$.

3.2.   Tail behaviours for interval distributions for cases B-D

For these distributions, we again focus on their tail behaviours. A proof for the next proposition can be found in Appendix A.

Proposition 3.5. The following statements hold:

(a) Density $g\left(y, a, b\right)$ has a fat left tail and a fat right tail for case B for all $b>0$, and for case C if $b>1.$ For case D, it does not have a fat right tail nor a fat left tail for any $b>0.$

(b) The tailed index of $g\left(y, a, b\right)$ for both right and left tails is 1 for case B for all $b>0$, and is $(1-\frac{1}{b})$ for case C for B for $b>1$.

Remark 3.6. Among distributions A, B, C, and Beta distribution, distribution B gets the highest tailed index of 1, independent of the choices of $b>0$.

4.   Algorithms for fitting interval distribution models

In this section, we assume that $\mathrm{\Phi }$ in (2.2) is a function from ${R}^{1}$ to $(0, 1)$ with positive continuous derivatives. We focus on parameter estimation algorithms for model (1.1).

First, we consider a simple case, where risk outcome $y$ is driven by a model:

where $b>0$ is a constant, $v = {a}_{0}+{a}_{1}{x}_{1}+\dots +{a}_{k}{x}_{k}$, and $s \sim N\left(0, 1\right),$ independent of fixed effects ${x = (x}_{1}, {x}_{2}, \dots , {x}_{k})$. The function $\mathrm{\Phi }$ does not have to be the standard normal CDF. But when $\mathrm{\Phi }$ is the standard normal CDF, the expected value $E\left(y\right|x)$ can be evaluated by the formula ${E}_{S}\left[\mathrm{\Phi }\left(a+bs\right)\right] = \mathrm{\Phi }\left(\frac{a}{\sqrt{1+{b}^{2}}}\right)$ ^[12].

Given a sample ${\left\{({x}_{1i}, {x}_{2i}, \dots , {x}_{ki}, {y}_{i})\right\}}_{i = 1}^{n},$ where $({x}_{1i}, {x}_{2i}, \dots , {x}_{ki}, {y}_{i})$ denotes the ${i}^{th}$ data point of the sample, let ${z}_{i} = {\mathrm{\Phi }}^{-1}\left({y}_{i}\right).$ and ${v}_{i} = {a}_{0}+{a}_{1}{x}_{1i}+\dots +{a}_{k}{x}_{ki}.$ By (2.3), the log-likelihood function for model (4.1) is:

where $f$ is the density of $s.$ The part of $\sum _{i = 1}^{n}\mathrm{l}\mathrm{o}\mathrm{g}{\rm{ \mathsf{ ϕ} }}\left({z}_{i}\right)$ is constant, which can be dropped off from the maximization.

Recall that the least squares estimators of ${{a}_{0}, a}_{1}, \dots , {a}_{k}$, as a row vector, that minimize the sum squares

has a closed form solution given by the transpose of ${\left({\mathrm{X}}^{T}\mathrm{X}\right)}^{-1}{\mathrm{X}}^{\mathrm{T}}\mathrm{Z}$ ^[5,9] whenever the design matrix $\mathrm{X}$ has a rank of $k,$ where

The next proposition shows there exists an analytical solution for the parameters of model (4.1).

Proposition 4.1. Given a sample ${\left\{({x}_{1i}, {x}_{2i}, \dots , {x}_{ki}, {y}_{i})\right\}}_{i = 1}^{n}$, assume that the design matrix has a rank of $k.$ If $s \sim N(0, 1),$ then the maximum likelihood estimates of parameters ${(a}_{0}, {a}_{1}, \dots , a)$, as a row vector, and parameter $b$ are respectively given by the transpose of ${\left({\mathrm{X}}^{T}\mathrm{X}\right)}^{-1}{\mathrm{X}}^{\mathrm{T}}\mathrm{Z},$ and ${b}^{2} = \frac{1}{n}\sum _{i = 1}^{n}{({z}_{i}-{v}_{i})}^{2}.$ In absence of fixed effects ${\{x}_{1}, {x}_{2}, \dots , {x}_{k}\}$, parameters ${a}_{0}$ and ${b}^{2}$ degenerate respectively to the sample mean and variance of ${z}_{1}$, ${z}_{2}, \dots , {z}_{n}.$

Proof. Dropping off the constant term from (4.2) and noting $f\left(z\right) = \frac{1}{\sqrt{2\pi }}\mathrm{exp}\left(-\frac{{z}^{2}}{2}\right)$, we have

Hence the maximum likelihood estimates $({a}_{0}, {a}_{1}, \dots , {a}_{k})$ are the same as least squares estimators of (4.3), which are given by the transpose of ${\left({\mathrm{X}}^{T}\mathrm{X}\right)}^{-1}{\mathrm{X}}^{\mathrm{T}}\mathrm{Z}.$ By taking the derivative of (4.4) with respect to $b$ and setting it to zero, we have ${b}^{2} = \frac{1}{n}\sum _{i = 1}^{n}{({z}_{i}-{v}_{i})}^{2}.$

Next, we consider the general case of model (1.1), where the risk outcome $y$ is driven by a model:

where parameter $w$ is formulated as $w = \mathrm{exp}\left(u\right),$ and $u = {b}_{0}+{b}_{1}{x}_{1}+ \dots +{b}_{k}{x}_{k}.$ We focus on the following two cases:

(a) $s \sim N\left(0, 1\right),$

(b) $s$ is standard logistic.

Given a sample ${\left\{({x}_{1i}, {x}_{2i}, \dots , {x}_{ki}, {y}_{i})\right\}}_{i = 1}^{n},$ let ${w}_{i} = \mathrm{e}\mathrm{x}\mathrm{p}({b}_{0}+{b}_{1}{x}_{1i}+ \dots +{b}_{k}{x}_{ki})$ and ${{u}_{i} = b}_{0}+{b}_{1}{x}_{1i}+ \dots +{b}_{k}{x}_{ki}.$ The log-likelihood functions for model (4.5), dropping off the constant part $\mathrm{log}\left({U}_{2}\right),$ for cases (a) and (b) are given respectively by (4.6) and (4.7):

Recall that a function is log-concave if its logarithm is concave. If a function is concave, a local maximum is a global maximum, and the function is unimodal. This property is useful for searching maximum likelihood estimates.

Proposition 4.2. The functions (4.6) and (4.7) are concave as a function of ${(a}_{0}, {a}_{1}, \dots , {a}_{k})$. As a function of ${(b}_{0}, {b}_{1}, \dots , {b}_{k}),$ (4.6) is concave.

Proof. It is well-known that, if $f(x)$ is log-concave, then so is $f(Az + b),$ where $Az + b$ : ${R^m} \to {R^1}$ is any affine transformation from the m-dimensional Euclidean space to the 1-dimensional Euclidean space. For (4.6), the function $f\left(x\right) = -({z-v)}^{2}\mathrm{exp}(-2u)$ is concave as a function of $v,$ thus function (4.6) is concave as a function of ${(a}_{0}, {a}_{1}, \dots , {a}_{k})$. Similarly, this function $f\left(x\right)$ is concave as a function of $u,$ so (4.6) is concave as a function of ${(b}_{0}, {b}_{1}, \dots , {b}_{k}).$

For (4.7), the linear part $-\left({z}_{i}-{v}_{i}\right)\mathrm{e}\mathrm{x}\mathrm{p}(-{u}_{i})$, as a function of ${(a}_{0}, {a}_{1}, \dots , {a}_{k})$, in (4.7) is ignored. For the second part in (4.7), we know $-\mathrm{log}\{1+\mathrm{e}\mathrm{x}\mathrm{p}[-(z-v)/\mathrm{e}\mathrm{x}\mathrm{p}\left(u\right)\left]\right\}$, as a function of $v$, is the logarithm of the CDF of a logistic distribution. It is well-known that the CDF for a logistic distribution is log-concave. Thus (4.7) is concave with respect to ${(a}_{0}, {a}_{1}, \dots , {a}_{k})$.

In general, parameters ${(a}_{0}, {a}_{1}, \dots , {a}_{k})$ and ${(b}_{0}, {b}_{1}, \dots , {b}_{k})$ in model (4.5) can be estimated by the algorithm below.

Algorithm 4.3. Follow the steps below to estimate parameters of model (4.5):

(a) Given ${(b}_{0}, {b}_{1}, \dots , {b}_{k})$, estimate ${(a}_{0}, {a}_{1}, \dots , {a}_{k})$ by maximizing the log-likelihood function;

(b) Given ${(a}_{0}, {a}_{1}, \dots , {a}_{k})$, estimate ${(b}_{0}, {b}_{1}, \dots , {b}_{k})$ by maximizing the log-likelihood function;

(c) Iterate (a) and (b) until a convergence is reached.

5.   Conclusions

With the interval distributions introduced in this paper, models with a random effect can be fitted for a continuous risk outcome by maximum likelihood approaches assuming an interval distribution. These models provide an alternative regression tool to the Beta regression model and fraction response model, and a tool for tail risk assessment as well.

Acknowledgments

Authors are very grateful to the third reviewer for many constructive comments. The first author is grateful to Biao Wu for many valuable conversations. Thanks also go to Clovis Sukam for his critical reading for the manuscript.

We would like to thank you for following the instructions above very closely in advance. It will definitely save us lot of time and expedite the process of your paper's publication.

Conflict of interest

The views expressed in this article are not necessarily those of Royal Bank of Canada and Scotiabank or any of their affiliates. Please direct any comments to Bill Huajian Yang at h_y02@yahoo.ca.