Alexander Meister
email: meistear@mathematik.uni-stuttgart.de

Optimal convergence rates for density estimation from grouped data

Keywords: aggregated data, cross-validation, deconvolution, minimax convergence rates, nonparametric estimation

We derive the optimal convergence rates for density estimation based on aggregated observations under common smoothness conditions for symmetric densities.

We study a procedure for data-driven bandwidth selection and give an extension to skew densities.

A c c e p t e d m a n u s c r i p t 1. Introduction

The statistical problem of estimating a density function when only aggregated data are observed has received considerable attention, where research is mainly motivated by real applications in the field of econometrics.

In the mathematical model, we define the i.i.d. random variables X ij , integer i, j, having the density function f X , which we aim to estimate. The only empirical access is

given by the data Y 1 , . . . , Y n where

Y i = m j=1 X ij , (1)
where m is the fixed size of the data groups. Note that our problem may also be seen as a missing data problem in time series analysis. Consider the moving average process (Z k ) k with Z k = m+k-1 j=k X j with i.i.d. X 1 , X 2 , . . . having the density f X . Then, assume that only the data Z 1+mk are observed for integer k ≥ 0.

The current note intends to advance the understanding of this problem by deriving the optimal rates of convergence under common smoothness conditions on f X with respect to the mean integrated squared error (MISE). In Section 2, we study the optimal rates for symmetric densities. In Section 3, we describe a procedure for data-driven bandwidth selection; in Section 4 we give an extension to skew densities. The proofs are deferred to Section 5.

To give a survey on related problems of indirect density estimation, we mention problems of reconstructing a density from measurement error; that topic has become famous as density deconvolution [START_REF] Stefanski | Deconvoluting kernel density estimators[END_REF], [START_REF] Carroll | Optimal rates of convergence for deconvolving a density[END_REF], [START_REF] Zhang | Fourier methods for estimating mixing densities and distributions[END_REF], [START_REF] Horowitz | Semiparametric estimation of regression models for panel data[END_REF], [START_REF] Fan | On the optimal rates of convergence for nonparametric deconvolution problems[END_REF] among many others). See the recent note of [START_REF] Machado | A note on identification with averaged data[END_REF] for a study on parameter identification from averaged data. A nonparametric approach to density estimation from aggregated observations is given by [START_REF] Linton | Nonparametric estimation with aggregated data[END_REF], where the authors consider a more complicated model with an additional independent error component for each data group. Pointwise asymptotic normality of the authors' estimator is studied at certain convergence rates. Another contribution is given by [START_REF] Schick | Root n consistent density estimation for sums of independent random variables[END_REF] where estimation of the density of the sums of independent random variables is studied; hence, somehow, our consideration can be seen as the corresponding inverse problem. [START_REF] Belomestny | Rates of convergence for constrained deconvolution problem[END_REF] (2005).

Minimax rates

We focus on those f X which are symmetric around zero and, hence, have a real-valued Fourier transform f f t X (t). As the density f Y of each observation Y i turns out to be the m-fold self-convolution of f X , it is convenient to apply a Fourier approach. With respect to the characteristic functions f f t X and f f t Y , we have

f f t Y (t) = f f t X (t) m .
We assume that f X (x) = f X (-x) for almost all x; and, in addition,

f f t X (t) = 0 , ∀t . (2)
The necessity of a condition as (2) to ensure identifiability of f X is shown in the following:

We define the density

f 0 (x) = 1 -cos(x) /(πx 2 ) , (3)
having the triangle-shaped Fourier transform f f t 0 (t) = 1 -|t| on t ∈ [-1, 1]. Therefore, the densities f ± (x) = f 0 (x) ± (1/2)f 0 (x) cos(2x) possess the Fourier transforms

f f t ± (t) = f f t 0 (t) ± 1 4 f f t 0 (t + 2) + f f t 0 (t -2) , hence we have f f t + (t) 2 = f f t -(t)
2 for all t. This proves that f X = f + cannot be uniquely reconstructed from the observation density f Y = f + * f + in this example; where * denotes convolution.

The empirical Fourier transform is denoted by f f t Y (t) = 1 n n k=1 exp(itY k ). We define the estimator of f X by Fourier inversion,

fX (x) = 1 2π exp(-itx)K f t (th) f f t Y (t) 1/m dt , (4)

A c c e p t e d m a n u s c r i p t

where K denotes a square-integrable kernel function with K f t (0) = 1 and K f t is compactly supported; parameter h denotes the bandwidth.

In order to establish rates of convergence, we propose common smoothness conditions on f X by assuming a uniform upper bound on the Sobolev norm of f X . We introduce the class F βC of even densities which satisfy (2) and

|f f t X (t)| 2 (1 + t 2 ) β dt ≤ C ,
where β describes the smoothness degree. Further, we consider so-called supersmooth densities whose Fourier transforms satisfy

|f f t X (t)| 2 exp C 0 |t| γ dt ≤ C 1 .
Those densities are collected into the class G C 0 C 1 γ . We give the following theorem.

Theorem 1 Let • denote the L 2 (R)-norm. As the kernel function K, we choose the sinc kernel with

K f t (t) = χ [-1,1] (t), i.e. the indicator function of the interval [-1, 1].
(a) Take estimator f X as in (4). Then, we have 2β+1) where c n > 0 is bounded away from both ∞ and 0;

sup f X ∈F βC E fX -f X 2 = O n -2β/ m(2β+1) , sup f X ∈G C 0 C 1 γ E fX -f X 2 = O (ln n) 1/γ n -1/m , when selecting h = c n n -1/ m(
and h = d(ln n) -1/γ with a constant d ≤ C 1/γ 0 , respectively. (b) Assume an arbitrary estimator f of f X based on the data Y 1 , . . . , Y n . Then, for γ ∈ (0, 1) and β > 1/2, C, C 1 sufficiently large, there is a constant c > 0 so that sup f X ∈F βC E f -f X 2 ≥ c • n -2β/ m(2β+1) , sup f X ∈G C 0 C 1 γ E f -f X 2 ≥ c • (ln n) 1/γ n -1/m .
Hence, we have established rate optimality of our estimation procedure. We notice deterioration of the convergence rate compared to density estimation based on direct data

A c c e p t e d m a n u s c r i p t

where the well-known rates n -2β/(2β+1) and (ln n) 1/γ n -1 , resp., occur; they are included into our framework for m = 1. The rates become worse when m increases. We mention that Theorem 1(a) can be extended to more general kernels K as long as

K f t (t) = 1 on
an open interval around t = 0 and K f t is compactly supported.

Adaptive estimation

The choice of the bandwidth h as given in Theorem 1 leads to optimal rates; however it requires knowledge of the parameters β, γ, C 0 . Therefore our goal is to find a fully data-driven bandwidth selector. In classical density estimation, cross-validation (CV) is a famous procedure for adaptive bandwidth choice. As mentioned in [START_REF] Linton | Nonparametric estimation with aggregated data[END_REF] there is no straight-forward extension of the underlying theory to aggregated data problems as the estimators are non-linear.

Nevertheless, we can apply CV to estimate the observation density f Y based on the direct data Y 1 , . . . , Y n . The outcome bandwidth is denoted by ĥC , see [START_REF] Hall | Extent to which least-squares cross-validation minimises integrated square error in nonparametric density estimation[END_REF] for the methodology and theory for that problem. In this section we restrict our consideration to densities whose Fourier transforms satisfy

C 2 |t| -β-1/2 ≤ |f f t X (t)| ≤ C 3 |t| -β-1/2 , ∀|t| ≥ T (5)
for some T; and |f f t X (t)| ≥ |f f t X (T )| for all |t| ≤ T . This smoothness assumption is closely related to f X ∈ F βC with appropriate constants; indeed, the optimal convergence rates are the same under the corresponding constraints. Therefore, the mean integrated squared error for the estimation of f Y is minimised by h = h 0 ∼ n -1/(2βm+m) . Surprisingly, this selection rule also minimises the MISE in our aggregated data problem when estimating f X according to Theorem 1. That inspires us to employ ĥ = ĥC as the bandwidth selector.

The resulting estimator is denoted by fX, ĥ. With respect to the convergence rates, we give a weak individual version. We write const. for a generic positive constant.

Proposition 1 Assume f X satisfies (5) where β > 7/2; and f X , f X are integrable. We apply the sinc kernel in fX, ĥ. Then, for all c > 0, we have

lim sup n→∞ P n 2β/ m(2β+1) fX, ĥ -f X 2 > c ≤ const. • c -1

A c c e p t e d m a n u s c r i p t

Therefore, the adaptive estimator fX, ĥ keeps the optimal rates from Theorem 1 under certain circumstances. The case of supersmooth f X is more difficult to address.

Skew densities

When f X is no longer assumed to be symmetric around zero, its Fourier transform f f t X is not real-valued. Therefore the inversion procedure becomes more difficult as we have

m different complex roots of f f t Y (t)
and its empirical version f f t Y (t). Let R(t), ϕ(t) denote the absolute value and the angle of f f t Y (t) in the polar representation of complex numbers. We face the problem that the angle is not uniquely defined. Our intention is to specify ϕ(t) so that those functions are continuous for all t and ϕ(0) = 0.

We introduce the intervals I j = ((j/2 -1)π, (j/2 + 1)π] for integer j. Assuming that R(t) = 0, for all t, the angle ϕ(t) is uniquely determined by f f t Y (t) if its image is restricted to I j for any j. Therefore, we denote the angle within I j by ϕ j (t). Setting ϕ(t) = ϕ j(t) (t), we start with t 0 = 0, j(t 0 ) = 0. Then, given t k , we denote by t k+1 the smallest t > t k where ϕ j(t k ) (t) crosses either (j(t k )/2 -1/2)π or (j(t k )/2 + 1/2)π; in the first case, we put j(t k+1 ) = j(t k ) -1; in the latter case, we define j(t k+1 ) = j(t k ) + 1; otherwise, put t k+1 = ∞. It follows from there that ϕ(t k ) = j(t k )π/2 holds for any k. The sequence (t k ) k>0 tends to infinity as, otherwise, the continuity of f f t Y (t) is violated at the limit of (t k ) k>0 . Then we define ϕ(t) = ϕ j(τ (t)) (t) with τ (t) = max{t k : t k ≤ t} for t ≥ 0; for t < 0 we set ϕ(t) = -ϕ(-t). Then ϕ is a continuous function on the whole real line with ϕ(0) = 0.

When determining an empirical version φ(t) for ϕ(t), we must consider that f f t Y (t) may have some isolated zeros. Therefore we introduce a parameter ρ n ∈ (0, 1). We realise that the set N = {t : | f f t Y (t)| ≤ ρ n } may be written as the disjoint union of countably many intervals, [τ j , τ j+1 ] say. We introduce a function f f t Y (t) which is equal to f f t Y (t) outside the set N ; while on any interval [τ j , τ j+1 ] we put f f t Y equal to the shortest connection between

τ j , f f t Y (τ j ) and τ j+1 , f f t Y (τ j+1 ) under the constraint | f f t Y (t)| = ρ n for all t ∈ [τ j , τ j+1 ]. Then, define R(t) = | f f t Y (t

A c c e p t e d m a n u s c r i p t

Further, we define the empirical version of f f t X (t) by

f f t X (t) = R1/m (t) exp it φ(t)/m , (6)
and, accordingly,

fX (x) = 1 2π exp(-itx)K f t (th) f f t X (t) dt , (7)
while stipulating that

K f t is supported on [-1, 1].
In order to give convergence rates we need more restrictive conditions compared to symmetric densities, namely |x| 2 f X (x)dx ≤ C 7 and (5). As an analogue for (5) for supersmooth densities we use

C 5 |t| (γ-1)/2 exp(-C 4 |t| γ ) ≤ |f f t X (t)| ≤ C 6 |t| (γ-1)/2 exp(-C 4 |t| γ ) , ∀|t| ≥ T . (8)
Densities satisfying those conditions are collected into the class F βC 7 C 2 C 3 , which corresponds to F βC . When assuming (8) instead of ( 5), we call the density class

G γC 7 C 4 C 5 C 6 , as the analogue of G C 0 C 1 γ .
Proposition 2 Take estimator fX as defined in ( 7) and K as in Theorem 1. Choose

ρ n = exp(-n). Under the constraint f X ∈ F βC 7 C 2 C 3 , select h = h n = [C h (ln n)/n] 1/[m(2β+1)] with a constant C h > (576m)/C 2m 2 . Then, for β > 1/2, we have sup f X ∈F βC 7 C 2 C 3 E fX -f X 2 = O (ln n/n) 2β/[m(2β+1)] . If f X ∈ G γC 7 C 4 C 5 C 6 and max{1, γ} < m, select h = h n = 2C 4 m/ 1-ν(ln ln n)/ ln n 1/γ • (ln n) -1/γ with ν ∈ (1 -m(γ -1)/γ, m/γ] to obtain sup f X ∈G γC 7 C 4 C 5 C 6 E fX -f X 2 = O (ln n) 1/γ n -1/m .
Therefore, the rates are kept in the case of f X ∈ F βC 7 C 2 C 3 from the symmetric constraint f X ∈ F βC (see Theorem 1) up to a logarithmic factor; while, for f X ∈ G γC 7 C 4 C 5 C 6 , they are exactly the same as for symmetric f X ∈ G C 0 C 1 γ in Theorem 1.

A c c e p t e d m a n u s c r i p t

Proofs

Proof of Theorem 1: (a) By Parseval's identity and Fubini's theorem, we obtain, for

f X ∈ F βC , E fX -f X 2 = |K f t (t/h)| 2 E | f f t Y (t)| 1/m -|f f t X (t)| 2 dt + O(h 2β ) . (9)
In the case of

f X ∈ G C 0 C 1 γ , the bias term in (9) changes from O(h 2β ) to O exp(-C 0 h -γ ) .
For the variance term we use the inequalities |x 1/my 1/m | m ≤ |x -y| for x, y > 0 and even m as well as |x 1/my 1/m | m ≤ 2 m |x -y| for all x, y in the case of odd m; combined with Jensen's inequality, we obtain

E | f f t Y (t)| 1/m -|f f t X (t)| 2 ≤ const. • E f f t Y (t) -f f t Y (t) 2 1/m = O n -1/m ,
independently of t. Therefore, the mean integrated squared error is bounded above by

O h -1 n -1/m , h 2β for f X ∈ F βC and O h -1 n -1/m , exp(-C 0 h -γ ) for f X ∈ G C 0 C 1 γ . Choos-
ing h as stated in the theorem leads to the given rates.

(b) First we consider f X ∈ F βC . Take f 0 as in (3). Furthermore, we introduce the supersmooth Cauchy density

f 1 (x) = π -1 (1 + x 2 ) -1 with f f t 1 (t) = exp(-|t|).
From there, we construct the following subclass of densities

f n,θ (x) = 1 2 f 1 (x) + f 0 (x) + const. • 2kn≥j≥kn θ j j -β-(1/2) cos(2jx)f 0 (x) , (10)
where k n denotes a positive integer still to be determined and const. is sufficiently small; and all θ j ∈ {0, 1} are i.i.d. random variables with P (θ j = 0) = 1/2. The corresponding Fourier transforms are given by

f f t n,θ (t) = 1 2 f f t 1 (t) + f f t 0 (t) + const. • 2kn≥j≥kn θ j j -β-(1/2) f f t 0 (t -2j) + f f t 0 (t + 2j) .
By that we may verify that all f f t n,θ are non-vanishing and f n,θ ∈ F βC . Now assume an arbitrary estimator f (x) = f (x; Y 1 , . . . , Y n ). In the sequel, we write θ j,b = θ kn , . . . , θ j-1 , b, θ j+1 , . . . , θ 2kn where b ∈ {0, 1} and f Y,θ j,b = f * ,n n,θ j,b , i.e. the n-fold self-convolution of f n,θ j,b , which, therefore, denotes the density of each observation Y j . We

A c c e p t e d m a n u s c r i p t

consider its mean integrated squared error, using Parseval's identity and Fubini's theorem, sup

fx∈F βC E f -f X 2 ≥ const. 2kn≥j≥kn E θ E f n,θ j 2j+1 2j-1 f f t (t) -f f t n,θ (t) 2 dt ≥ const. 2kn≥j≥kn E θ 2j+1 2j-1 f f t n,θ j,0 (t) -f f t n,θ j,1 (t) 2 dt • • • • min{f Y,θ j,0 (y 1 ) • • • f Y,θ j,0 (y n ) , f Y,θ j,1 (y 1 ) • • • f Y,θ j,1 (y n ) dy 1 • • • dy n ≥ const. 2kn≥j≥kn j -2β-1 ≥ const.k -2β n , (11)
if the χ 2 -distance between f Y,θ j,0 and f Y,θ j,1 satisfies

χ 2 (f Y,θ j,0 , f Y,θ j,1 ) = f Y,θ j,0 (x) -f Y,θ j,1 (x) 2 f Y,θ j,0 (x) -1 dx = O(1/n) (12)
holds for all k n ≤ j ≤ 2k n . There we have used a result in [START_REF] Fan | Adaptively local one-dimensional subproblems with application to a deconvolution problem[END_REF] which has been developed for the classical deconvolution problem. That method was influenced by Donoho & Liu (1991a,b). From the definition of the density subclass, we derive f Y,θ j,0 (x) ≥

2 -m f * ,m 1 (x). That implies f Y,θ j,0 (x) ≥ const.(1 + x 2 ) -1 ,
and, from there, we see that (12) may be represented equivalently in the Fourier domain by

2j+1 2j-1 g j (t) 2 + g j (t) 2 dt = O(1/n) , (13)
where

|g j (t)| = f Y,θ j,1 -f Y,θ j,0 f t (t) = 1 2 f f t 1 (t) + j -β-1/2 f f t 0 (t + 2j) m - 1 2 f f t 1 (t) m ≤ const. • j m(-β-1/2) ,
while j ≥ k n and |t| ≥ 2j -1. We can derive the same bound for |g j (t)|. Then, (13) and, hence, (11) follow when selecting

k n = n 1/[m(2β+1)] .
In the case f X ∈ G C 0 C 1 γ , we replace j -β-(1/2) by n -1/(2m) and set k n = (1/4) • ln n/(mC ) 1/γ with some constant C > C 0 in (10). Then the proof follows analogously.

A c c e p t e d m a n u s c r i p t

Proof of Proposition 1: We consider

P n 2β/ m(2β+1) • fX, ĥ -f X 2 > c ≤ P n 2β/ m(2β+1) |t|≤(dh 0 ) -1 f f t X, ĥ(t) -f f t X (t) 2 dt > c/2 + P ( ĥ < dh 0 ) + P n 2β/ m(2β+1) ĥ2β > c/2 (14)
for some d ∈ (0, 1). The first addend is seen to be bounded above by const.•c -1 when using

Markov's inequality. The second and the third term are bounded above by P | ĥ -

h 0 | ≥ const. • h 0 .
Both n 3/10 ( ĥ0h 0 ) and n 3/10 ( ĥ -ĥ0 ) are asymptotically normal distributed, due to Theorem 2.1 in [START_REF] Hall | Extent to which least-squares cross-validation minimises integrated square error in nonparametric density estimation[END_REF]; where ĥ0 is the bandwidth which minimises the integrated squared error; the conditions are satisfied for β > 7/2, f X , f X integrable and an appropriate kernel; also note that the kernel used for deriving ĥC need not coincide with K in estimator fX, ĥ. It follows from there that the second and the third term in ( 14) converge to zero as n → ∞ so that the proposition follows.

Proof of Proposition 2: Note that, for |t| ≤ 1/h, we have

| f f t Y (t) -f f t Y (t)| ≥ | f f t Y (t) -f f t Y (t)| -2ρ n and hence | f f t Y (t) -f f t Y (t)| ≥ (1/2)| f f t Y (t) -f f t Y (t)|. With f f t X as in (6), we derive E f f t X (t) -f f t X (t) 2 = E R1/m (t) exp i φ(t)/m -R 1/m (t) exp iϕ(t)/m 2 ≤ E R1/m (t) -R 1/m (t) 2 + E R1/m (t)R 1/m (t) 1 -cos ( φ(t) -ϕ(t))/m = O n -1/m + E R1/m (t)R 1/m (t) 1 -cos ( φ(t) -ϕ(t))/m • χ {| φ(t)-ϕ(t)|>π} + E R1/m (t)R 1/m (t) 1 -cos ( φ(t) -ϕ(t))/m • χ {| φ(t)-ϕ(t)|≤π} .
Considering the fact that the function 1cos(•/m) / 1cos(•) is bounded on the interval [-π, π], some standard techniques lead to The term above as well as the third addend in ( 16) are negligible when selecting M = C M h -m(β+1/2)-1 with an appropriate constant C M > 0.

E f f t X (t) -f f t X (t) 2 = O n -1/m , R 2
In the case f X ∈ G γC 7 C 4 C 5 C 6 , the proof follows by replacing C 2 by C 5 , h 2β by exp(-2C 4 h -γ ) in ( 15) and h m(β+1/2) by exp(-C 4 mh -γ ) • h -m(γ-1)/2 in the first two lines of ( 16) and, accordingly, in the sequel.

Then the specific choice of h as stated in the theorem leads to the desired rates in the view of (15).

  studies a procedure for estimating the density of a component of a moving average process (MA(1)) where the scaled with some factor ϑ with |ϑ| < 1. The problem of estimating the density of independent components of a Poisson sum is considered in van Es et al.

  leads to E fXf X 2 = O h -1 n -1/m , h 2β , sup |t|≤1/h P (| φ(t)ϕ(t)| > π) f X ∈ F βC 7 C 2 C 3 ; the latter term needs more careful consideration for any t ∈ [-1/h, 1/h]. Consider that | φ(t)ϕ(t)| > π implies the existence of at least one s ∈ [0, t] so that | φ(s)ϕ(s)| = π/2 and hence f f t Y (s)f f t Y (s) ≥ (1/2) • |f f t Y (s)|.Applying the sequence (s j ) j , j = 1, . . . , M -1 of equidistant points s j = jt/M where M = M n , we derive thatP (| φ(t)ϕ(t)| > π) ≤ P ∃j = 1, . . . , M -1 : f f t Y (s j )f f t Y (s j ) ≥ (C m 2 /6) • h m(β+1/2) ≥ const. • h m(β+1/2) + χ {h -1 M -1 ≥const.•h m(β+1/2) } .(16)By Hoeffding's inequality, the first addend in (16) has the upper boundO(M ) • exp -1 288 C 2m 2 h m(2β+1) n .By Markov's inequality, the second addend is bounded above by χ {h -1 M -1 ≥const.•h m(β+1/2) } + O(1/n)


Petroula M Mavrikiou

Kolmogorov inequalities for the partial sum of independent bernoulli random variables

Keywords:

In this paper two Kolmogorov inequalities are presented for the sample average of independent (but not necessarily identically distributed) Bernoulli random variables.

Introduction

For a sequence of independent and identically distributed Bernoulli random variables X 1 , X 2 , ..., with E(X 1 ) = p , [START_REF] Kolmogorov | On tables of random numbers[END_REF] provided the following inequality

P (sup k≥n | Xk -p| > ) ≤ 2e -2n 2 (1-) where Xk = 1 k k i=1 X i , > 0.
Improvements, extensions and related results can be found in [START_REF] Hoeffding | Probability inequalities for sums of bounded random variables[END_REF], [START_REF] Young | A note on a Kolmogorov inequality[END_REF], [START_REF] Turner | Improved Kolmogorov inequalities for the binomial distribution[END_REF], [START_REF] Young | A note on a Kolmogorov inequality for the Binomial distribution[END_REF], [START_REF] Christofides | Probability inequalities with exponential bounds for U-statistics[END_REF], [START_REF] Christofides | A Kolmogorov inequality for U-statistics based on Bernoulli kernels[END_REF], [START_REF] Kambo | On exponential bounds for binomial probabilities[END_REF] and [START_REF] Banjevic | On a Kolmogorov inequality[END_REF]. In this paper, two Kolmogorov inequalities are provided for the case of independent but not necessarily identically distributed Bernoulli random variables.

Preliminaries

The following results will be used:

Lemma 2.1

Let p 1 , ..., p n be positive real numbers and p = 1

n n i=1 p i . Then for t > 0 n i=1 (p i e t + 1 -p i ) ≤ (pe t + 1 -p) n .
Proof: From the arithmetic-geometric mean inequality, we have

n i=1 (p i e t + 1 -p i ) ≤ [ n i=1 1 n (p i e t + 1 -p i )] n = (pe t + 1 -p) n .
1

A c c e p t e d m a n u s c r i p t

The following result is due to [START_REF] Christofides | A Kolmogorov inequality for U-statistics based on Bernoulli kernels[END_REF].

Lemma 2.2

Let < 1 2 and

g(p, ) = (1 -p -)ln( 1 -p - 1 - p ) + (p + )ln( p + p ). Then for p + < 1 2 or 1 2 + 1 3 ≤ p ≤ 1 g(p, ) ≥ - 1 2 ln(1 -4 2 ). Lemma 2.3 Let x = 2(p + ) -1 and y = 1 -2p. Then, ∞ r=1 1 2r(2r -1) [x 2r + (2r -1)y 2r + 2rxy 2r-1 ] = g(p, )
where g(p, ) is the function defined in Lemma 2.2.

Proof: We have p + = x+1 2 and p = 1-y 2 . Then, g(p, ) from Lemma 2.2 is

g(p, ) = ( 1 -x 2 )ln( 1 -x 1 + y ) + ( 1 + x 2 )ln( 1 + x 1 -y ).
Using the Taylor series expansions

ln(1 + x) = ∞ r=1 (-1) r-1 x r r and ln(1 -x) = - ∞ r=1
x r r we have that

g(p, ) = (1 -x) 2 [- ∞ r=1 x r r - ∞ r=1 (-1) r-1 y r r ]+( 1 + x 2 )[ ∞ r=1 (-1) r-1 x r r + ∞ r=1 y r r ]

A c c e p t e d m a n u s c r i p t

and after algebraic manipulations we arrive at

g(p, ) = ∞ r=1 1 2r(2r -1) [x 2r + (2r -1)y 2r + 2rxy 2r-1 ]. Lemma 2.4
Let ν be a positive integer and x > 1. Then

2 2ν-2 (x 2ν-1 + 1) -(x + 1) 2ν-1 > 0. Proof: Let F (x) = 2 2ν-2 (x 2ν-1 + 1) -(x + 1) 2ν-1 . Then F (x) = 2 2ν-2 (2ν -1)x 2ν-2 -(2ν -1)(x + 1) 2ν-2 = (2ν -1)[2 2ν-2 x 2ν-2 -(x + 1) 2ν-2 ] = (2ν -1)[(2x) 2ν-2 -(x + 1) 2ν-2 ] > 0 since x > 1.
Thus, F (x) is an increasing function and F (x) > F (1) = 0.

Lemma 2.5

Let y ≥ 1 and ν = 1, 2, ... . Then

H(y) = y 2ν + 2ν -1 + 2νy - 4 2 2ν (y + 1) 2ν ≥ 0. Proof: H (y) = 2ν.y 2ν-1 + 2ν - 4 2 2ν 2ν(y + 1) 2ν-1 = 2ν 4 2 2ν [ 2 2ν 4 y 2ν-1 + 2 2ν 4 -(y + 1) 2ν-1 ]
A c c e p t e d m a n u s c r i p t

= 2 3-2ν .ν[2 2ν-2 y 2ν-1 + 2 2ν-2 -(y + 1) 2ν-1 ] = 2 3-2ν .ν[2 2ν-2 (y 2ν-1 + 1) -(y + 1) 2ν-1 ].
By Lemma 2.4 H (y) > 0 implying that H is increasing and therefore H(x) ≥ H(1) = 0.

Lemma 2.6

Let x, y be as in Lemma 2.3 and r = 1, 2, ... . Then

x 2r + (2r -1)y 2r + 2rxy 2r-1 ≥ 4( x + y 2 ) 2r .
Proof: By lemma 2.5 for c ≥ 1

c 2r + 2r -1 + 2rc - 4 2 2r (c + 1) 2r ≥ 0.
Take c = x y . Then x 2r y 2r + 2r -1 + 2r

x y ≥ 4 2 2r ( x y + 1) 2r

x 2r + (2r -1)y 2r + 2rxy 2r-1 ≥ 4( x + y 2 ) 2r .

Main Results

Theorem 3.1

Let Y 1 , Y 2 , ..., Y n be a sequence of independent Bernoulli random variables with E(Y i ) = p i , i = 1, ..., n and < 1 2 . Then, for p+ < 1 2

or p ≥ 1 2 + 1 3 , P ( Ȳ -p > ) ≤ (1 -4 2 ) n 2 where Ȳ = 1 n n i=1 Y i and p = 1 n n i=1 p i .

A c c e p t e d m a n u s c r i p t

Proof: Let s > 0. Then,

P ( Ȳ -p > ) = P [s( Ȳ -p -) > 0] ≤ E[e s( Ȳ -p-) ] = e -s(p+ ) E[e (s Ȳ ) ] = e -s(p+ ) E(e s 1 n n i=1 Y i ) = e -s(p+ ) n i=1 E(e s n Y i ) = e -s(p+ ) n i=1 (p i e s n + 1 -p i ) ≤ e -s(p+ ) (pe s n + 1 -p) n (1) = e -f (s)
where f (s) = s(p + ) -nln(pe s n + 1 -p) and the last inequality follows from Lemma 2.1.

The function f is maximized at s max = nln[ (1-p)(p+ ) p(1-p-) ] and

f (s max ) = n(p + )ln (p + )(1 -p) p(1 -p -) -nln 1 - p 1 -p - = n[(1 -p -)ln( 1-p- 1-p ) + (p + )ln( p+ p )] = n.g(p, ).
Thus, by Lemma 2.2

f (s max ) = ng(p, ) ≥ - n 2 ln(1 -4 2 )
A c c e p t e d m a n u s c r i p t and therefore

P ( Ȳ -p > ) ≤ (1 -4 2 ) n 2 .
Notice that the restrictions on p and arise from Lemma 2.2 which is due to [START_REF] Christofides | A Kolmogorov inequality for U-statistics based on Bernoulli kernels[END_REF].

The following theorem gives an exponential bound under different conditions on p and .

Theorem 3.2

Let Y 1 , Y 2 , ..., Y n be a sequence of independent Bernoulli random variables, with E(Y i ) = p i , i = 1, ..., n. Then for p + > 1 2 or p < 1 2 and ∀ < 1,

P ( Ȳ -p > ) ≤ e -n(2 2 + 1 3 4 e 2 4 )
where Ȳ = 1 n n i=1 Y i and p = 1 n n i=1 p i .

Proof: From the proof of Theorem 3.1 (p, ) where

P ( Ȳ -p > ) ≤ e -ng
g(p, ) = (1 -p -)ln( 1 -p - 1 - p ) + (p + )ln( p + p ). By Lemma 2.3 g(p, ) = ∞ r=1 1 2r(2r -1) [x 2r + (2r -1)y 2r + 2rxy 2r-1 ].
x and y are positive values, imposing the restrictions on p and . Now, since x + y = 2 and using Lemma 2.6 we have Proof: Using Lemma 1 of Turner, [START_REF] Turner | A Kolmogorov inequality for the sum of independent Bernoulli random variables with unequal means[END_REF] we can arrive at inequality (1) of Theorem 3.1 having the quantity P [sup k≥n ( Ȳk -pk ) > ] as our left-hand side. Then, one can follow the same steps of theorems 3.1 and 3.2 to reach the right-hand side, and the proof is complete.

Remarks

(1) Theorem 3.3 provides sharper bounds than that of the main result of [START_REF] Turner | A Kolmogorov inequality for the sum of independent Bernoulli random variables with unequal means[END_REF], under of course restrictions on p and .

(2) Theorem 3.1 is an extension of Corollary 3.2 of [START_REF] Christofides | A Kolmogorov inequality for U-statistics based on Bernoulli kernels[END_REF].

Corresponding author: Tel: +357

345159 Ext 132; Fax: +357 22 438234 Email address: bus.mp@fit.ac.cy (P. M. Mavrikiou)

.4( 2 2 )

The last inequality is clearly true for r = 2, and can be established (e.g., by induction) for r > 2. Then,

Thus

and the proof is complete.

Theorem 3.3

Let Y 1 , Y 2 , ..., be a sequence of independent Bernoulli random variables with E(Y i ) = p i . Then,

(1) ∀ < 1 2 and for pk + < 1 2 or pk ≥

n 2 n = 0, 1, 2, ...


Georgi N Boshnakov
email: georgi.boshnakov@manchester.ac.uk

Some measures for asymmetry of distributions

Keywords: asymmetry, asymmetry curve, asymmetry index, confidence characteristic

We propose several measures, functional and scalar, for asymmetry of distributions by comparing the behavior of probability densities to the right and left of the mode(s) and show how to generate classes of equivalent distributions from a given distribution, allowing for varying asymmetry but retaining some information theoretic properties of the original distribution, such as the entropy.

Introduction

The numerical characteristic normally employed to characterize the lack of symmetry of a distribution is the coefficient of skewness, a standardized version of the third central moment. A variant of this takes expectation with respect to the median rather than the mean. Various measures of asymmetry are discussed by [START_REF] Macgillivray | Skewness and asymmetry: Measures and orderings[END_REF]. The monograph Dharmadhikari and Joag-dev (1988) does not treat asymmetry directly but is still a very useful reference for its comprehensive study of closely related topics.

We propose some alternatives which measure asymmetry with respect to modes rather than means or medians. The proposed measures always exist and seem quite intuitive. Our approach provides a systematic way to create, from a reference distribution, an entire family of distributions having different asymmetries but sharing some fundamental properties of the original distribution, such as the differential entropy. This may be of some interest for kernel estimation as a source of asymmetric kernels.

Notation

We consider absolutely continuous distributions and measure asymmetry by comparing how long does it take for the density to fall to a given value on the two sides of the modes. A very useful tool for this is the confidence transformation [START_REF] Boshnakov | Confidence characteristics of distributions[END_REF] which produces, from a given source distribution, a new distribution (called the confidence characteristic) whose density is, effectively, a rearrangement of the values of the source density in decreasing order. This transformation preserves some important properties of the source distribution and so it makes sense to say that distributions having the same confidence density belong to a family, see [START_REF] Boshnakov | Confidence characteristics of distributions[END_REF] for details and related notions. The distribution function and the probability density of the Given a distribution, we denote its distribution function, probability density, confidence distribution function, and confidence density by F , f , G, and g, respectively.

Asymmetry curves and coefficients

Suppose that f is unimodal with mode m. Let, for l > 0, x = x(l) and y = y(l) be such that f (m + x) = f (m -y) and l = x + y. Then x = x(l) and y = y(l) are monotonically increasing as functions of l > 0. If f is strictly monotonically decreasing on each side of the mode then l represents the length of the region where f is greater than f (m+x), while x and y represent the length of the part of this region where f is monotonically descreasing and increasing, respectively. To accomodate multi-modal distributions we use this property to define x(l) and y(l).

For any c > 0 let S c = {z : f (z) ≥ c}. Let also S I,c (respectively, S D,c ) be the subregion of S c where the density f is strictly increasing (decresing). Denote the lengths of these regions by l, l incr , and l decr , respectively. We will assume that l = l incr + l decr .

If f is symmetric then l decr = l incr for all l. So, for a symmetric distribution the parametric plot of l decr = l decr (l) and l incr = l incr (l) as a function of l is a straight line with slope one. For asymmetric distributions such a plot provides comprehensive information about asymmetry. So, we introduce the following definition.

Definition 1. The curve (l decr (l), l incr (l)), l ≥ 0, is said to be the asymmetry curve of f and (l decr (l), l incr (l) -l 2 ) its detrended asymmetry curve. The detrended asymmetry curves of symmetric distributions coincide with the positive x-axis. Asymmetry may be defined in terms of functions of (l decr (l), l incr (l)) as follows.

Definition 2. The functions l decr (l)/l, l incr (l)/l, and l decr (l)/l incr (l), where l > 0, are called right-asymmetry, left-asymmetry, and odds-asymmetry, respectively.

Symmetric distributions may be thought of as having a constant asymmetry equal to zero. More generally: Definition 3. A distribution is said to have constant asymmetry if its asymmetry curve is a straight line.

A c c e p t e d m a n u s c r i p t

For many common distributions the asymmetry curve for small l is close to a straight line with slope one, indicating "approximate" symmetry in a neghbourhood of the mode. Figure 1 shows the asymmetry curves of several Gamma distributions. It makes sense to say that the distribution whose curve is "closest" to the line with slope 1 is the most symmetric one, Gamma(1,120) in this case. This is expected here since increasing the second parameter of the Gamma distribution results in a distribution closer to normal, see also the example in Section 4.3.

Various summary characteristics of the above functions may be considered candidates for the title "coefficient (or index) of asymmetry". Let

r pos = E g l decr l , r neg = E g l incr l , r assym = r pos -r neg , (1)
where E g denotes expectation with respect to the confidence density g(l). We will call r assym , r pos , and r neg the mean asymmetry, the mean positive asymmetry and the mean negative asymmetry, respectively. The expectations above exist since l decr l and l incr l are positive and less than one. Moreover, Theorem 1. The coefficients r pos , r neg , and r assym are always finite, r pos and

r neg are in the interval [0, 1], r assym is in [-1, 1].
The mean odds-asymmetry is defined by

r odds = E g l decr (l) l incr (l) , ( 2
)
and infinity is a possible value for it.

From the definitions above it is easy to see that symmetric distributions have constant asymmetry.

Theorem 2. Let f be symmetric with median m, i.e., f (m + x) = f (m -x). Then f has constant asymmetry; the right-asymmetry, left-asymmetry, r pos , and r neg are equal to l 2 ; r odds = 1, r assym = 0. Non-symmetric distributions may also have constant asymmetry. A typical example may happen when the density has several modes and decreases symmetrically around each one of them.

It is clear that the introduced measures of asymmetry are invariant with respect to a shift. Some are invariant with respect to change of scale as well, others are not. Indeed, let Y = cX, where c > 0 is a positive constant, X is a random variable. We will use the above notation with an additional index

A c c e p t e d m a n u s c r i p t

x or y for the probability characteristics of X and Y . The two densities are related by f y (y) = 1 c f x ( y c ). Hence f y (a) = f y (a + l) if and only if f x (a/c) = f x ((a + l)/c). In the particular case when f x is unimodal this shows that l x is transformed to cl y . Similarly, it can be seen that l y,incr (l) = cl x,incr (l/c) and l y,decr (l) = cl x,decr (l/c). Thus a change of scale leads, in general, to a scale change in the discussed asymmetry measures. It is easy to see however that the change of the asymmetry curve, for example, corresponds to changing the units of its plot. The odds-asymmetry is invariant under a scale transformation. The asumptions of unimodality can be removed with the help of Corollary 1 from [START_REF] Boshnakov | Confidence characteristics of distributions[END_REF].

Asymmetry of some distributions 4.1 Symmetric unimodal distributions

Let F be a symmetric unimodal distribution function with mode M . Then its confidence density is g

(l) = f (M + l 2 ) = f (M -l
2 ), i.e., l decr (l) = l incr (l) = l 2 in this case. So, the asymmetry curve is ( l 2 , l 2 ), r pos = r neg = 1 2 , r assym = 0, r odds = 1.

Triangular distribution

Let

f (x) = 2x/H 0 ≤ x ≤ H, 2(1 -x)/(1 -H) H ≤ x ≤ 1,
where H ∈ (0, 1). The system f (H -l incr ) = f (H +l decr ), l = l decr +l incr , gives l decr (l) = l(1 -H) and l incr (l) = lH. So, l decr (l) = 1-H H l incr (l), l decr (l)/l = 1 -H, l incr (l)/l = H, l decr (l)/l -l incr (l)/l = 1 -2H, where l ∈ [0, 1]. Hence, the asymmetry of the triangular distribution is constant. The distribution is symmetric if H = 1 2 , skewed to the right if H < 1 2 and skewed to the left otherwise. The odds-asymmetry is equal to 1-H H .

Γ-distribution

Let f be a Γ-density [START_REF] Johnson | Continuous univariate distributions[END_REF],

f (x) = λ α Γ(α) x α-1 e -λx , x ≥ 0.
(3)

A c c e p t e d m a n u s c r i p t

We assume here that α > 1. In this case the distribution is unimodal with mode M = (α -1)/λ such that

f (M ) > f (x), for every x = M . ( 4
)
The equation f (M + l decr ) = f (M -l incr ), where l decr ≥ 0 and l incr ≥ 0, will be satisfied if

(M + l decr ) α-1 e -λ(M +l decr ) = (M -l incr ) α-1 e -λ(M -l incr ) ,
which can be written as

M + l decr M -l incr α-1 = e λ(M +l decr -M +l incr ) = e λl .
(5)

Also, M + l decr M -l incr = (M -l incr ) + (l incr + l decr ) M -l incr = 1 + l M -l incr . (6)
From equations ( 5) and ( 6) we get

1 + l M -l = e λl/(α-1) .
Hence, l M -l incr = e λl/(α-1) -1.

So,

M -l incr = l e λl/(α-1) -1 .

Finally, using the identity l decr + l incr = l, we get

l incr = M - l e λl/(α-1) -1 l decr = l -M + l e λl/(α-1) -1
Hence, when α > 1 the confidence density of the Γ-distribution is

g(l) = f (M + l decr ), l ≥ 0 = λ α Γ(α) (l + l e λl/(α-1) -1 ) α-1 e -λ(l+ l e λl/(α-1) -1 )

A c c e p t e d m a n u s c r i p t

From the way we defined l decr and l incr , they should have the following limiting behavior

l decr ---→ l→∞ ∞ l decr --→ l→0 0 l incr ---→ l→∞ M l incr --→ l→0 0
This is indeed so, since it is easy to verify that for any c > 0

l e cl -1 ---→ l→∞ 0 and l e cl -1 --→ l→0 1 c .
We also have g(0) = f (M ), as expected. Thus, the odds-asymmetry tends to infinity as l → ∞, the right-asymmetry and the asymmetry tend to 1. We see that the asymmetry curve of the Γ-distribution depends only on M = (α -1)/λ. In other words, for Γ-distributions having the same mode the asymmetry curves and all measures of asymmetry derived from it are the same. For comparison, the usual coefficient of skewness is equal to 2/ √ α and so depends on α but not on λ.

Distributions with given asymmetry

The confidence transformation has a number of desirable properties. In particular, it preserves the entropy and other information theoretic properties of the original distribution, see [START_REF] Boshnakov | Confidence characteristics of distributions[END_REF] for details. It is therefore justifiable to classify distributions by their confidence characteristics. By reverting the above process distributions with specified asymmetry properties may be generated.

Constant asymmetry

Suppose that g is the confidence density of some unimodal distribution and we wish to create a distribution with the same confidence characteristic but with odds-asymmetry c > 0. Using the established notation, the following relations should be satisfied:

l decr l incr = c, l decr + l incr = l, f (M + l decr ) = f (M -l incr ) = g(l). (7)
Let x ≥ 0, y ≥ 0. The required density is

f (M + x) = g( 1+c c x) f (M -y) = g((1 + c)y).

An asymmetric normal family

The confidence distribution function and confidence density of the standard normal distribution are G(l) = 2Φ(l/2) -1 and g(l/2) = ϕ(l/2) = exp(-l 2 /8)/ √ 2π. The above formulae then give

f (M + x) = 1 √ 2π e -( 1+c c ) 2 x 2 /8 , f (M -y) = 1 √ 2π e -(1+c) 2 y 2 /8 .
The distributions obtained by varying c and M form a family of distributions having constant asymmetry and the same confidence characteristic as the standard normal distribution.

Another family

It may be more convenient in some circumstances to express l and l incr in terms of l decr . So, let

l incr = u(l decr ), l = l decr + u(l decr ), (8)
where u(.) is an appropriate function. Then we may define a new density f by f (M + l decr ) = f (M -l incr ) = g(l) = g(l decr + u(l decr )). ( 9)

For example if we take g(z) = λe -λz to be the the exponential density and u(x) = x 2 , then we get

f (M + l decr ) = f (M -l 2 decr ) = g(l decr + u(l decr )) = λe -λ(l decr +l 2 decr ) .
This can be written as

f (M + l decr ) = g(l decr + l 2 decr ) = λe -λ(l decr +l 2 decr ) , if l decr ≥ 0 g(|l decr | + |l decr |) = λe -λ(|l decr |+ √ |l decr |) , if l decr < 0.
Now the asymmetry is not constant. For example, the odds-asymmetry is l decr /l 2 decr = 1/l decr .

Conclusion

We defined some measures of asymmetry which provide useful information about this type of property. The asymmetry curve and its variants provide comprehensive information, while their averaged counterparts summarize it to single numbers. These measures provide a systematic way to generate asymmetric distributions. Figure 1: Asymmetry curves, (l decr (l), l incr (l)), l > 0, of Gamma(1,4), Gamma(1,20) and Gamma(1,120). When the second parameter is increased, the curves are closer to the straight line with slope 1 (almost the same in the case of Gamma(1,120) over the plotted range), the interpretation being that the corresponding distributions become more symmetric. All pictured curves however have horizontal asymptotes because the left limit of the support of the Gamma distribution is finite.

  called the confidence distribution function and confidence density, respectively.


Eric Beutner

Udo Kamps

Random convex combinations of order statistics

Keywords: Order statistics, Random convex combinations, Beta distribution, Record values

For order statistics X i,n , X j,n , X k,n with 1 ≤ i < j < k ≤ n, and

is established for uniform distributions as well as seen to be characteristic.

Introduction

Let X 1,n ≤ . . . ≤ X n,n be the order statistics based on a random sample of size n from a continuous distribution function F . These are denoted by U 1,n ≤ . . . ≤ U n,n , if the underlying distribution is uniform. * Email addresses: beutner@stochastik.rwth-aachen.de, kamps@stochastik.rwthaachen.de

1

A c c e p t e d m a n u s c r i p t

The deterministic convex combination

(1 -c) • X i,n + c • X i+1,n , c ∈ (0, 1), (1)
of order statistics appears in distribution-free confidence intervals for quantils (cf. [START_REF] David | Order Statistics[END_REF], Chapter 7.1).

Motivated by (1), [START_REF] Jones | On fractional uniform order statistics[END_REF] shows that some fractional order statistic Ũi+c,n , c ∈ (0, 1), (see [START_REF] Stigler | Fractional order statistics, with applications[END_REF]) in the uniform case may be represented by a random convex combination

(1 -C) • U i,n + C • U i+1,n , ( 2
)
where C is beta distributed with parameters c and 1 -c, Beta(c, 1 -c) for short, and independent of the order statistics U i,n and

U i+1,n , 1 ≤ i ≤ n -1,
from the uniform distribution on (0, 1), denoted by U (0, 1).

By considering equation (2) the question arises, whether an (ordinary) order statistic U j,n from a uniform distribution can be represented by a random convex combination of its neighbours.

Random convex combinations

It is shown that the distributional equation

U j,n d = (1 -C) • U i,n + C • U k,n (3)
holds true for uniform order statistics with 1 ≤ i < j < k ≤ n, where C is an appropriately chosen beta distributed random variable, independent [START_REF] Aggarwala | Progressive Censoring. Theory, Methods, and Applications[END_REF], p. 32). Hence,

It is well known that U r,n d = 1 - r l=1 B l , where B l ∼ Beta(n -l + 1, 1), 1 ≤ l ≤ n, (cf.
(1 -C) • U i,n + C • U k,n d = (1 -C) • 1 - i l=1 B l + C • 1 - k l=1 B l = 1 - i l=1 B l 1 -C • (1 - k l=i+1 B l ) . From k l=i+1 B l ∼ Beta(n-k +1, k -i) which implies (1- k l=i+1 B l ) ∼ Beta(k - i, n -k + 1), we have C • (1 - k l=i+1 B l ) ∼ Beta(j -i, n -j + 1). Therefore, (1 -C • (1 - k l=i+1 B l )) ∼ Beta(n -j + 1, j -i). Since i l=1 B l ∼ Beta(n -i + 1, i) we obtain that i l=1 B l • 1 -C • (1 - k l=i+1 B l ) ∼ Beta(n -j + 1, j)

A c c e p t e d m a n u s c r i p t

which implies the assertion because U j,n is beta distributed with parameters j and n -j + 1.

Due to relation (3), a random convex combination of order statistics may serve as an unbiased prediction of a missing observation of some order statistic U j,n .

By analogy with Theorem 1, a related result may also be stated for record values from an exponential distribution.

Remark 2 Let R i , R j and R k , 0 ≤ i < j < k, be records from the exponential distribution with parameter λ. Then, for n = 0, 1, 2, . . . , R n d

= 1 λ • n l=0 X l ,
where {X l , l ≥ 0} is a sequence of i.i.d. standard exponential distributed random variables (cf. [START_REF] Arnold | Records[END_REF], p. 20). Hence,

(1 -C) • R i + C • R k d = (1 -C) • λ -1 • i l=0 X l + C • λ -1 • k l=0 X l = λ -1 • i l=0 X l + C • λ -1 • k l=i+1 X l . ( 4
)
The first summand in ( 4) is gamma distributed with parameters i + 1 and λ.

If C is independent of the records and beta distributed with parameters j -i and k-j , then the second term in (4) is gamma distributed with parameters j -i and λ. Therefore, (4) is gamma distributed with parameters j + 1 and

λ which implies R j d = (1 -C) • R i + C • R k .

A c c e p t e d m a n u s c r i p t

The following Theorems provide characterizations of uniform distributions based on the distributional equation (3).

Theorem 3 Let X j-1,n , X j,n , X j+1,n be order statistics from an absolutely continuous and strictly increasing distribution function F , and let C be independent of X j-1,n and X j+1,n . Furthermore, let the distribution of C be concentrated on (0, 1) with an absolutely continuous distribution function

F C . Then (X j-1,n , Z, X j+1,n ) d = (X j-1,n , X j,n , X j+1,n ) where Z = (1 -C) • X j-1,n + C • X j+1,n iff C ∼ U (0, 1) and F ≡ U (a, b), for some a < b.
Proof. Due to the assumed independence of C and X j-1,n , X j+1,n the joint distribution function f X j-1,n ,Z,X j+1,n of Z and X j-1,n , X j+1,n is

f X j-1,n ,Z,X j+1,n (x j-1 , z, x j+1 ) = 1 x j+1 -x j-1 f X j-1,n ,X j+1,n (x j-1 , x j+1 ) × f C z -x j-1 x j+1 -x j-1 I {x j-1 ≤z≤x j+1 } = 1 x j+1 -x j-1 • n! (j -2)!(n -j -1)! F j-2 (x j-1 )f (x j-1 )(F (x j+1 ) -F (x j-1 )) ×f (x j+1 )(1 -F (x j+1 ) n-j-1 f C z -x j-1 x j+1 -x j-1 I {x j-1 ≤z≤x j+1 } .
The joint density of X j-1,n , X j,n , X j+1,n is given by

f X j-1,n ,X j,n ,X j+1,n (x j-1 , z, x j+1 ) = n! (j -2)!(n -j -1)! F j-2 (x j-1 )f (x j-1 ) × f (z)(1 -F (x j+1 )) n-j-1 f (x j+1 ) • I {x j-1 ≤z≤x j+1 } .

A c c e p t e d m a n u s c r i p t

Equating the densities leads to

F (x j+1 ) -F (x j-1 ) x j+1 -x j-1 • f C z -x j-1 x j+1 -x j-1 = f (z), ∀ x j-1 ≤ z ≤ x j+1 . (5)
It is directly seen that ( 5) holds true for F ≡ U (a, b) and C ∼ U (0, 1), which proves the if-part.

To prove the only if-part we divide equation ( 5) by F (x j+1 )-F (x j-1 ), which is greater than zero due our assumption, and integrate both sides from x j-1 to u to obtain

F C (v) -F C (0) = F (ṽ) -F (x j-1 ) F (x j+1 ) -F (x j-1 ) , ∀ v ∈ (0, 1), ∀ x j-1 ≤ x j+1 , (6) if we put v = u-x j-1
x j+1 -x j-1 and ṽ = v(x j+1 -x j-1 ) + x j-1 .

By differentiating (6) with respect to x j+1 we find vf (ṽ)(F (x j+1 ) -F (x j-1 )) = (F (ṽ) -F (x j-1 ))f (x j+1 ), ∀ x j-1 ≤ x j+1 , which implies

F C (v) -F C (0) = v • f (v(x j+1 -x j-1 ) + x j-1 ) f (x j+1 ) , ∀ x j-1 ≤ x j+1 . ( 7
)
Since the left hand side in ( 7) is independent of x j-1 , we obtain by differentiating (7) with respect to

x j-1 v f (x j+1 ) f (ṽ)(1 -v) = 0, ∀ x j-1 ≤ x j+1 , ∀ v ∈ (0, 1).
Hence, f must be constant, i.e., there exists a < b such that F ≡ U (a, b).

From ( 7) we now obtain that C ∼ U (0, 1).

A c c e p t e d m a n u s c r i p t

Theorem 4 Let X j-1,n , X j,n , X j+1,n be order statistics from an absolutely continuous distribution function F , which is assumed to be strictly increasing and convex or concave. Moreover, let C ∼ U (0, 1) and independent of X j-1,n and X j+1,n . Then

(1 -C) • X j-1,n + C • X j+1,n d = X j,n (8)
iff F ≡ U (a, b), for some a < b.

Proof. Notice that

P ((1 -C) •X j-1,n + C • X j+1,n ≤ x) = P ((1 -C) • x j-1 + C • x j+1 ≤ x)dP X j-1,n ,X j+1,n (x j-1 , x j+1 ) = x -x j-1 x j+1 -x j-1 I {x j-1 ≤x≤x j+1 } dP X j-1,n ,X j+1,n (x j-1 , x j+1 ) and P (X j,n ≤ x) = P (X j,n ≤ x|X j-1,n = x j-1 , X j+1,n = x j+1 )dP X j-1
,n ,X j+1,n (x j-1 , x j+1 ) = F (x) -F (x j-1 ) F (x j+1 ) -F (x j-1 ) I {x j-1 ≤x≤x j+1 } dP X j-1,n ,X j+1,n (x j-1 , x j+1 ).

The distributional equation ( 8) now implies 0 = F (x) -F (x j-1 ) F (x j+1 ) -F (x j-1 ) -x -x j-1 x j+1 -x j-1 ×I {x j-1 ≤x≤x j+1 } dP X j-1,n ,X j+1,n (x j-1 , x j+1 ).

A c c e p t e d m a n u s c r i p t

Due to the fact that F is assumed to be strictly increasing as well as convex or concave (cf. [START_REF] Hewitt | Real and Abstract Analysis[END_REF], p. 272) we conclude F (x) -F (x j-1 ) F (x j+1 ) -F (x j-1 ) -x -x j-1 x j+1 -x j-1 = 0, ∀ x j-1 ≤ x ≤ x j+1 .

Hence, F is linear, i.e., there exists a < b such that F ≡ U (a, b).

  statistics. Moreover, it is shown that uniform distributions are characterized by relation (3). Theorem 1 Let U i,n , U j,n and U k,n be order statistics from the uniform distribution on (a, b) with 1 ≤ i < j < k ≤ n, and let the random variable C be independent of the order statistics with C ∼ Beta(j -i, k -j). Then relation (3) is valid. Proof. Without loss of generality, a standard uniform distribution may be assumed, since multiplying relation (3) by b -a, b > a, and adding a then leads to the required distributional identity for order statistics from the uniform distribution on (a, b).


Fabio Bellini
email: fabio.bellini@unimib.it

Leonardo Bottolo
email: l.bottolo@imperial.ac.uk

Stationarity domains for -power garch process with heavy tails

Keywords: Power Garch models, Strict and weak stationarity, α-stable distributions, Heavy tails distributions

 can be extended to general innovations, regardless of the existence of their δ-moments. We prove some general properties of these domains and analyze some cases particularly relevant in applications.

Power Garch models

The δ-power Garch(1, 1) model is a natural extension of the usual Garch(1, 1) model, described by the equations

X t = σ t ε t σ δ t = w + α |X t-1 | δ + βσ δ t-1 (1)
with ω > 0, α ≥ 0, β ≥ 0, δ > 0 and ε t independent and identically distributed (i.i.d.). The case δ = 2 with ε t standard normal corresponds to the usual Garch(1, 1) model [START_REF] Bollerslev | Generalized autoregressive conditional heteroschedasticity[END_REF] while the case δ < 2 with ε t α-stable has been recently discussed in [START_REF] Mittnik | Stationarity of stable power-GARCH processes[END_REF], MP&R hereafter. The model [START_REF] Bellini | Domain and misspecification issues in fitting Garch(1, 1) models[END_REF] in its generality has been considered

A c c e p t e d m a n u s c r i p t amongst others in [START_REF] Ding | A long memory property of stock market returns and a new model[END_REF] and [START_REF] Liu | Maximum likelihood estimation of a Garch-stable model[END_REF]. The specification of a varying exponent δ has the merit of allowing "long memory" in the shocks of conditional variance (see [START_REF] Ding | A long memory property of stock market returns and a new model[END_REF], for a broad discussion).

As in the case of traditional Garch models, one general question is for which parameters' values we have a stationary solution (in the weak or strong sense) of the equations [START_REF] Bellini | Domain and misspecification issues in fitting Garch(1, 1) models[END_REF]. More precisely, we define the δSS domain in the (α, β) space as the region where a strictly stationary solution of (1) exists, while the δWS domain is the region where we have a weakly stationary (covariance-stationary) solution. For the usual case δ = 2 we will simply refer to 'SS domain' and 'WS domain'. For these models typically the δSS domain is bigger than the δWS domain, since for some parameters' values there are stationary solutions with infinite unconditional variance. This is in contrast with the situation of traditional (light-tailed) time series models, where strict stationarity usually implies weak stationarity.

Below we briefly review some classic and more recent results on stationarity domains. [START_REF] Nelson | Stationarity and persistence in the Garch(1,1) model[END_REF] in the case δ = 2 proved the following theorem (where ε is the common

distribution of ε t ). Theorem If ω > 0, if E[log(αε 2 + β)
] exists (possibly infinite) and ε 2 is non degenerate, then a necessary and sufficient condition for the strict stationarity of [START_REF] Bellini | Domain and misspecification issues in fitting Garch(1, 1) models[END_REF] is that

E[log(αε 2 + β)] < 0 (2)
Remark 1 For this characterization of the SS domain it is not necessary to assume that E[ε 2 ] < +∞; see for example the second part of Theorem 6 in [START_REF] Nelson | Stationarity and persistence in the Garch(1,1) model[END_REF] where ε has a Cauchy distribution.

Remark 2 If E[ε 2 ] = 1, a sufficient condition for (2) is α + β < 1 (3) since E[log(αε 2 + β)] ≤ E[(αε 2 + β -1)] = α + β -1 < 0.
As it is very well known since Bollerslev's seminal paper, [START_REF] Bollerslev | Generalized autoregressive conditional heteroschedasticity[END_REF], the condition (3) also identifies exactly the WS domain. Hence for the case δ = 2 with standard normal innovations the WS domain is strictly included in the SS domain.

A c c e p t e d m a n u s c r i p t

Remark 3 In [START_REF] Nelson | Stationarity and persistence in the Garch(1,1) model[END_REF] the SS domain has been explicitly computed for ε a standard normal or Cauchy distribution. See [START_REF] Bellini | Domain and misspecification issues in fitting Garch(1, 1) models[END_REF] and the discussion in Section 4 on the impact of the shape of the SS domain in the MLE.

In MP&R the following result has been proved (albeit in a more general form for Garch(p, q) models). We restate here Proposition 1 of MP&R for the Garch(1, 1) case.

Theorem If ε is α-stable with 1 < α S < 2 and δ < α S , then a necessary and sufficient condition for strict stationarity is

E[log(α |ε| δ + β)] < 0 ( 4
)
Remark 4 In this case the sufficient condition for the δWS domain becomes (Proposition 2, MP&R)

αλ α S ,β S ,δ S + β < 1 (5)
where λ α S ,β S ,δ S = E[|ε| δ ] < +∞, since δ < α S with α S , β S and δ S the characteristic exponent or tail index, the skewness and the location parameter of a stable distribution.

This condition follows immediately using the same argument as before:

E[log(α |ε| δ +β)] ≤ E[α |ε| δ + β -1] = αλ α S ,β S ,δ S + β -1 < 0.
We can think of this region as the natural counterpart of the WS region for the δ = 2 case.

Since the hypothesis of the Bougerol and Picard (1992a) theorem do not require integrability of the innovations but only of their logarithm, it seems natural to extend this theorem to a general δ-power Garch, by removing the two assumptions of MP&R that ε is α-stable and δ < α S . In the footpaths of MP&R we will find that in general the condition

E[log(α |ε| δ + β)] < 0
is necessary and sufficient for strict stationarity provided that the integral exists, while if the innovations have finite δ-moment E[|ε| δ ] < +∞ a sufficient condition is given by

αλ α S ,β S ,δ S + β -1 < 0
The rest of the paper is organized as follows: in Section 2 we will extend MP&R's results while in Section 3 we will show some properties of the general δSS together with some illustrative examples. Finally Section 4 will contain some concluding remarks and comments.

A c c e p t e d m a n u s c r i p t 2 The δSS domain for general δ-power Garch models

We follow the same strategy of Bougerol and Picard (1992b) and MP&R. The general δ-power Garch (p, q) model can be written as

X t = σ t ε t σ δ t = ω + p i=1 α i |X t-j | δ + q j=1 β j σ δ t-i (6)
with ω > 0, α i ≥ 0, β j ≥ 0, and ε t i.i.d. with E[ε t ] = 0 while no additional assumptions on the existence of other moments of ε t are made. It is of course possible to add a timevarying mean µ t but we avoid it for notational simplicity. The equations ( 6) can be written in the following form

Y t+1 = A t Y t + B ( 7
)
where

A t =            α 1 |ε| δ + β 1 1 0 ... 0 α 2 |ε| δ + β 2 0 1 ... 0 ... ... ... ... ... α m-1 |ε| δ + β m-1 0 0 ... 1 α m |ε| δ + β m 0 0 ... 0            , B =            ω 0 ... 0 0            (8)
and m = min(p, q).

The advantages of this representation with respect to the original work of Bougerol and Picard (1992b) are described in MP&R. In this way the problem becomes the existence of the solution of a stochastic recurrence equation of the form

Y n+1 = A n Y n + B n where
A n and B n are i.i.d. random matrices. The basic result is the following (see Bougerol and Picard, 1992b).

Theorem If E[log + A n ] < +∞, E[log + B n ] < +∞,
λ = inf n>0 1 n E[log A 1 A 2 • • • A n ] (9)
A c c e p t e d m a n u s c r i p t is negative.

Remark 5 We define log + (x) = max(0, log x) and • any matrix norm on R d . The most important point for us is that only a logarithmic moment condition is required, that can be fulfilled also without other (power) moment conditions. In the δ-power Garch(1, 1) case, condition (9) simply becomes

λ = E[log(α |ε| δ + β)] < 0
Using this theorem it is easy to prove the following generalized form of Mittnik's result. 3 Some general properties of the δSS domain

In the previous Section we saw how it is possible to generalize Nelson's and Mittnik's results in order to cover cases of very heavy tails. In this Section we will analyze more deeply the δSS domains by means of some numerical simulations. In the following remarks we point out some general properties of the δSS domain for δ-power Garch(1, 1).

Remark 7 Scaling the innovations with a factor η > 0 has the effect of inflating the δSS domain along the α-axis by a factor η -δ . Hence in order to compare different domains it is necessary to fix a "scale parameter" for the innovations. When δ = 2 it is natural to set E[ε 2 ] = 1 but for other innovations it is not clear which is the best "normalization value". Theorem 2 If ω > 0, the δSS domain is a convex subset of the (α, β) space whose boundary is given by the axes and a decreasing function β(α) starting from (0, 1). If

E[|ε| δ ] < +∞, the boundary is tangent in (0, 1) to the δWS domain while if E[|ε| δ ] = +∞
it has vertical tangent in (0, 1). Finally, we have

β (α) = V ar[Y ] α 2 (E[Y ]) 3
where

Y = 1 α |X| δ + β
Proof The convexity follows immediately from the concavity of the logarithm. Since the boundary of the δSS domain is given by

G(α, β) = E[log(α |ε| δ + β)] = 0 we have β (α) = - G α G β = - E[ |ε| δ (α|ε| δ +β)] ] E[ 1 (α|ε| δ +β)] ] < 0 and β (0) = -E[|ε| δ ]
since β(0) = 1. From Dini's theorem we have

β (α) = - [1, β (α)]HG[1, β (α)] T G β
where

HG = -    E[ |ε| 2δ (α|ε| δ +β) 2 ] E[ |ε| δ (α|ε| δ +β) 2 ] E[ |ε| δ (α|ε| δ +β) 2 ] E[ 1 (α|ε| δ +β) 2 ]    By posing X = |ε| δ (α|ε| δ +β)
and Y = 1 (α|ε| δ +β) 2 and noting that αX + βY = 1, we get

β (α) = [1, -E[X ] E[Y ] ]   E[X 2 ] E[X Y ] E[X Y ] E[Y 2 ]   [1, -E[X ] E[Y ] ] T E[Y ]
with straightforward computations we get the result.

A c c e p t e d m a n u s c r i p t

Apart from a qualitative characterization of the δSS for a δ-power Garch(1, 1) process given by Theorem 2, it is neither straightforward to derive the analytic expression of the δSS nor the δWS. However using Monte Carlo integration it is easy to represent the shape of both domains for a wide range of distributions commonly used in practice.

The normal case

In the normal case it is possible to calculate analytically the integral in (2) (see [START_REF] Nelson | Stationarity and persistence in the Garch(1,1) model[END_REF]) although the resulting expression is quite complex. If ε ∼ N (0, 1), β > 0, then

E[log(αε 2 + β)] = log(2α) + ϕ( 1 2 ) + ( 2πβ α ) 1 2 Φ( 1 2 , 3 2 ; β 2α ) -( β α ) 2 F 2 (1, 1; 2, 3 2 ; β 2α )
where

Φ(a, b; z) = ∞ k=0 (a) k (b) k z k k! with (a) k = Γ(a + k) Γ(a)
is the increasing factorial. The

The t-Student case

When ε has a normalized t-Student distribution with ν degrees of freedom it is in general not possible to compute analytically the integral in (2). Instead we computed SS domains for different values of ν using Monte Carlo integration with 100000 simulations from a t 3 , t 6 and t 9 . The results are showed in Fig. 2; the SS domain becomes bigger while the tail of the innovations become heavier, i.e. for smaller values of ν.

A c c e p t e d m a n u s c r i p t 3.3 The GED case

Another "classic" parametric family for Garch innovations is the Generalized Error Distribution (GED) with probability density function

f (x) = ν exp -1 2 x ν ν βΓ 1 ν 2 (1+ 1 ν ) with β = 2 -2 ν Γ 1 ν Γ 3 ν 1/2
where ν is a shape parameter; with ν = 2 we have a normal distributions, for ν < 2 the tails are heavier than the normal while for ν > 2 the tails are lighter. Again we used Monte Carlo integration with 100000 simulations from a GED distribution ν = 1.5, 3 and 5; the results are showed in Fig. 3. Again, the heavier the tails, the bigger the SS domain is. The following remark provides an heuristic explanation of this phenomenon.

Remark 9

We chose the t-Student family as a prototypical supergaussian family and the GED family as a prototypical subgaussian family. By visual comparing the various SS domains, we saw that in all cases the heavier the tails of the innovations, the bigger the SS domain is. A (quite crude) heuristic argument to motivate this behavior is as follows:

using the approximation log(1 + x) ∼ x -x 2 /2 in (2), it follows that

E[log(αZ 2 + β)] ∼ = α + β -1 - 1 2 (α + β -1) 2 - 1 2 α 2 E[Z 4 ] -1
showing that as E[Z 4 ] increases, the stationarity domain becomes bigger. However it is worth notice that this line of reasoning is valid only if the variance of the innovations is finite.

Remark 10 For all the distribution considered henceforth we have that E[ε 2 ] < +∞; it follows from Theorem 2 that the WS domain is tangent to the SS domain and strictly smaller.

The α-stable case

In this Subsection we consider the most important case for which E[ε 2 ] = +∞, the α-stable innovations. We use symmetric α-stable innovations with normalized scaling factor. More

A c c e p t e d m a n u s c r i p t

precisely, in the notations of [START_REF] Nolan | Stable Distributions -Models for Heavy Tailed Data[END_REF] we use the "0 parametrization" with different values of α S , β S = 0, γ S = 1 and δ S = 0, the characteristic exponent, the skewness, the scale and location parameter respectively. Also in this case the integral in (2) has been computed through Monte Carlo integration. The generation of stable variates has been performed with the STABLE package downloadable from the J. P. Nolan's homepage:

http://academic2.american.edu/∼jpnolan/stable/stable.html. As a special case, α S = 1, we have the Cauchy distribution; here it is easy to calculate analytically the SS domain (see [START_REF] Nelson | Stationarity and persistence in the Garch(1,1) model[END_REF] since if ε has a Cauchy distribution, then

E[log(αε 2 + β)] = 2 log( √ α + β)
Hence the boundary of the SS domain is given by the function

β(α) = 1 + α -2 √ α
It is easy to check that β (0) = -∞ as prescribed by Theorem 2.

For the sake of comparison the SS domains of stable innovations are plotted in Fig. 4 for α = 1 (Cauchy), α = 1.2, α S = 1.5, α S = 1.9 and α S = 2 (normal). It is interesting to note that in this case the SS domain becomes smaller when the tails of the innovations become heavier, i.e. when α S → 1 + . Hence the qualitative behavior is the opposite than in the case of finite variance innovations. This was argued in MP&R on the basis of the comparison of the corresponding δWS domains.

The case δ = 2

In the last Subsection we analyzed the shapes and the magnitudes of SS domains for usual Garch models, i.e. with δ = 2 in the "variance" equation, for general α-stable innovation with different values of the characteristic exponent α S . We found that there is a sharp distinction between the finite variance and the infinite variance case.

In this Subsection we want to investigate the situation for a generic δ. More precisely, for two fixed innovations, normal and α-stable with α S = 1.5, we compare the δSS domains for different values of δ. The results are shown in Fig. 5 for the normal case and in Fig. 6 for what concerns the 1.5-stable. In the normal case we have that as δ ↓ 1 the δSS region becomes significantly smaller, while in the stable case we have the opposite behavior. We we see that there is no qualitative change in the δSS domain, in agreement with the fact that the existence of the δ moment of the innovations is not important for the shape of the δSS domain.

Discussion

In this paper we showed that the strict stationarity condition for δ-power Garch models introduced by MP&R can be generalized to more general innovations, even with infinite δ moment. The two simplest examples of this behaviur are Cauchy innovations and α-stable innovations with α S < δ. In general, we saw that the δSS domain has a typical convex shape as provided by Theorem 2, with two possibilities according to the finiteness or not of the δ moment of the innovations.

In a previous paper, [START_REF] Bellini | Domain and misspecification issues in fitting Garch(1, 1) models[END_REF], we investigated the importance of the correct specification of the domain for the MLE. The usual practice is to consider only the δWS domain, whose shape doesn't depend on the density of the innovations but only on their δ moment. In the common case of δ = 2 this corresponds to the 'α 1 + β 1 < 1' constraint.

Our analysis shows two things:

• in the case of finite δ moment, the δSS domain can be much bigger than the δWS domain. Hence the choice of this boundary can affect the estimation process (for example giving rise to a spurious 'IGARCH' effect, see [START_REF] Bellini | Domain and misspecification issues in fitting Garch(1, 1) models[END_REF];

• this procedure is not applicable when the δ moment of the innovations is infinite. In this case there is no counterpart to the δWS region so it is necessary to determine numerically the δSS region, before performing MLE.

The same analysis should be extended to higher dimensions where to our knowledge the shape of SS domains is still not very well known even for normal innovations. for δ = 1, 1.5, 1.9 (blue) and 2 (red) from right to left.

A c c e p t e d m a n u s c r i p t Remark 8

 8 If |ε 2 | ≥ |ε 1 | a.s., then the δ SS domain of ε 1 is strictly bigger than that of ε 2 since E[log(α |ε 1 | δ + β)] < E[log(α |ε 2 | δ + β)] < 0 The general shape of the δSS domain in the two cases E[|ε| δ ] < +∞ and E[|ε| δ ] = +∞ is characterized in the following theorem.

function 2 F 2 F 2

 222 2 (a, b; c, d; z) is a generalized hypergeometric defined by (a, b; c, d; z) = ∞ k=0 (a) k (b) k (c) k (d) k z k k! while ϕ(z) is Euler 'psi' function with ϕ(1/2) = -1.96351. The SS and WS domains are shown in Fig. 1.

A c c e p t e d m a n u s c r i p t considered δ = 1

 1 , 1.5, 1.9, hence we have both situations E[|ε| δ ] < +∞ and E[|ε| δ ] = +∞;

Figure 1 :

 1 Figure 1: The SS domain (blue) and the WS domain (red) for standard normal innovations.

Figure 2 :

 2 Figure 2: The SS domains for a t-Student distribution with ν = 3, 6 and 9 degrees of freedom, from right to left, (black) compared with the normal SS domain (blue) and the WS domain (red).

Figure 3 :

 3 Figure 3: The SS domains for a GED distribution with ν = 1.5, 3 and 4, from left to right, (black) compared with the normal SS domain (blue) and the WS domain (red).

Figure 4 : 1 (

 41 Figure 4: The SS domains for α-stable innovations (black) with tail indexes α S = 1 (Cauchy), 1.2, 1.5, 1.9 and 2 (normal) (blue) from left to right.

Figure 5 :

 5 Figure 5: A comparison of SS domains for standard normal innovations for δ = 1, 1.5, 1.9 and 2, from left to right, (blue) and the WS domain (red).

Figure 6 :

 6 Figure 6: A comparison of SS domains for α-Stable innovations with tail index α S = 1.5

A c c e p t e d m a n u s c r i p t


Philip J Boland

Nóra Ní Chuív

Optimal Times for Software Release when Repair is Imperfect

Keywords: primary 90B25, secondary 68M20; 68N30 Software Testing, Software Reliability, Optimal Release Time, Birth Processes, Cost Models, Imperfect Repair

The determination of the optimal release time for a new piece of software is of primary importance in the process of software development. We study a model where initially there are N faults in the software each with failure detection rate λ(t), but where the probability of a perfect repair of a fault when found is p (in general repair is not perfect). We investigate various cost models for this situation and give some insight into how the optimal release times and costs for the software vary with p and λ(t).

Introduction

Software reliability is a crucial concern in our rapidly developing information technology world.

It is often defined as the probability of failure free operation of a computer programme in a specified environment for a specified period of time. For many years statisticians and software engineers have been developing models for failures and breakdowns in computer software (Farr, 1999;[START_REF] Musa | Software Reliability: Measurement, Prediction, Application[END_REF]Okumoto, 1987, Singpurwalla and[START_REF] Singpurwalla | Statistical Methods in Software Engineering: Reliability and Risk[END_REF][START_REF] Xie | Software Reliability Models -Past, Present and Future[END_REF][START_REF] Boland | Determining the Optimal Release Time for Software in the Geometric Poisson Reliability Model[END_REF].

A c c e p t e d m a n u s c r i p t

One may (broadly speaking) break the software life-cycle into the four stages of: requirements and specification, design, coding, and testing. Although the time spent in the different stages clearly varies with products and customers, it is not unreasonable to expect testing to take in the region of 30% of the development time. Failures in software are not physical in nature and usually result as a consequence of errors in logic. They can lead to results inconsistent with performance specifications and expectations, and hence knowledge about them often influences important economic decisions concerning software release. Software is tested in order to detect faults and consequently reduce the chances of failures occurring during its post release operation.

If software is tested for too little time it will involve high costs of customer dissatisfaction, while if it is tested for too long it will involve large costs of testing and opportunity lost as a result of releasing the software at a later stage.

Software reliability models may be broken into two broad categories ( In this paper we will restrict our attention to some generalizations of the Jelinski Moranda model. In this classic model one assumes there are a finite number of faults N in the software and that each fault has an equal and independent failure (discovery) rate of λ. Here one assumes that once a fault is discovered it is perfectly repaired and the new failure rate for the software is reduced by λ, and also that the time between the i th and (i + 1) st observed failures are exponentially distributed with parameter (N -i)λ. The total number of faults M (t) detected by time t is a birth process and letting

P n (t) = P (M (t) = n), it follows that P n (t) = N n e -λ(N -n)t (1-e -λt ) n and E(M (t)) = N (1 -e -λt ).
Imperfect Repair may occur in attempting to remove a fault from software. The Irish software testing company Lionbridge Technologies suggests that typically in the region of 15% of faults are not correctly repaired. We consider a modification of the Jelinski Moranda model to take account of imperfect repair and let p represent the probability that a fault which is discovered or detected is perfectly repaired. In extending earlier work of Boland and Ní Chuív (2001), we also assume that the detection (failure) rate of a fault is a function of time λ(t) which is not

Numbers of Faults Detected and Removed in Software Testing

We shall assume that initially there are N faults in the software each with failure rate λ(t),

resulting in an initial failure rate for the system of N λ(0). We let Λ(t) = t 0 λ(s) ds for t > 0, and in particular we will consider the fault detection rate functions λ(t) = be -bt and λ(t) = b/(1+bt), in addition to the standard λ(t) = b. When a fault is detected in the software at time t, we assume it is perfectly repaired with probability p, in which case the number of faults in the software is reduced by 1 (and consequently the failure rate is reduced by λ(t)). Otherwise the number of faults in the software (and the failure rate for the software) remains the same.

The probability density function for the time x k+1 that the (k + 1) th fault is successfully removed, given the time x k of removal of the k th fault, is given by

f (x k+1 |x k ) = p(N -k) λ(x k+1 ) e -(N -k)p [Λ(x k+1 )-Λ(x k )] (1)
Moreover the probability density function for the time x k of the k th fault to be removed is given by

f (x k ) = N ! (N -k)!(k -1)! p λ(x k ) e -pΛ(x k )(N -k+1) 1 -e -p Λ(x k ) k-1 = N ! (N -k)!(k -1)! e -p Λ(x k )(N -k) 1 -e -p Λ(x k ) k-1 d -e -p Λ(x k )
Note that this is the density of the k th order statistic, where the failures of each of the faults are independently identically distributed with failure rate λ(t).

For any interval of time (0, t], we let M (t) be the total number of faults detected or met by time t, and R(t) be the number of faults which have been removed (that is detected and perfectly repaired) by time t. The number of unsuccessful repairs in the interval (0, t] will be denoted by

U (t) = M (t) -R(t).
In our model the times between removals of faults are still independent (in the classic Jelinski Moranda model where λ(t) = b they are exponentially distributed with different parameters). For any integers u and r (where 0 ≤ r ≤ N ) we let P r (t) = P (R(t) = r)

and P u,r (t) = P (U (t) = u, R(t) = r). We derive useful expressions for these quantities in terms of t, Λ(t), and the parameters N and p. Using a standard differential equations approach, it may be established that for any

1 ≤ r ≤ N P r (t) = -(N -r) λ(t) p P r (t) + (N -r + 1) λ(t) p P r-1 (t). (2)

A c c e p t e d m a n u s c r i p t

With the initial conditions that P r (0) = 0 for 1 ≤ r ≤ N and P 0 (t) = e -N pΛ(t) for t ≥ 0, one may establish (by induction or otherwise) that

P r (t) = N r e -p(N -r) Λ(t) (1 -e -pΛ(t) ) r . ( 3
)
Note that R(t) is a binomial birth process (R(t) is Bin(N, 1 -e -ptΛ(t) ) for any t > 0), and when p = 1 (perfect repair) and λ(t) = b we obtain the standard Jelinski Moranda model.

We now derive an expression for P u,r (t) using an integral approach. We let

x 1 ≤ x 2 ≤ • • • ≤ x r
be the occurrence times of the r faults which are detected and successfully removed in the time interval (0, t], and u 0 , u 1 , . . . , u r be the number of faults which are encountered but not removed in the time intervals (0, x 1 ], (x 1 , x 2 ], . . . , (x r , t]. For convenience we will furthermore set x 0 = 0 and x r+1 = t. Letting x = (x 0 , x 1 , x 2 , . . . , x r ), u = (u 0 , u 1 , . . . , u r ), and u = Σ r i=0 u i , we have that the joint density function for x and u is given by

f (x, u) = p r r i=0 [(N -i)(Λ(x i+1 ) -Λ(x i ))] u i u i ! e -(N -i)(Λ(x i+1 )-Λ(x i )) (1 -p) u i r-1 i=0 (N -i)λ(x i ) = p r (1 -p) r i=0 u i N ! (N -r)! (e -r i=0 (N -i)(Λ(x i+1 )-Λ(x i )) ) r i=1 λ(x i ) • r i=0 [(N -i)(Λ(x i+1 ) -Λ(x i )] u i u i ! = p r (1 -p) u N r r! (e [-(N -r)Λ(t)-Λ(x r )-•••-Λ(x 1 )] ) • r i=1 λ(x i ) r i=0 [(N -i)(Λ(x i+1 ) -Λ(x i ))] u i u i ! .
On integration it follows that

P u,r (t) = u 0 +•••+u r =u t 0 xr 0 • • • x 2 0 f (x, u) dx 1 . . . dx r = N r r!p r (1 -p) u t 0 x r 0 • • • x 2 0 e -[(N -r)Λ(t)+Λ(x r )+•••+Λ(x 1 )] r i=1 λ(x i ) • 1 u! u 0 +•••+u r =u r i=0 u![(N -i)(Λ(x i+1 -Λ(x i ))] u i u i ! dx 1 . . . dx r .
Using the multinomial expansion

u 0 +•••+ur=u r i=0 u![(N -i)(Λ(x i+1 ) -Λ(x i ))] u i u i ! = r i=0 (N -i)(Λ(x i+1 ) -Λ(x i ))
u the expression for P u,r (t) simplifies to

P u,r (t) = N r p r (1 -p) u r! t 0 x r 0 • • • x 2 0 [(N -r)Λ(t) + Λ(x r ) + • • • + Λ(x 1 )] u u! •e -[(N -r)Λ(t)+Λ(xr)+•••+Λ(x 1 )] r i=1 λ(x i ) dx 1 • • • dx r .

A c c e p t e d m a n u s c r i p t

Hence

P r (t) = ∞ u=0 P u,r (t) = ∞ u=0 N r p r (1 -p) u r! • t 0 xr 0 • • • x 2 0 [(N -r)Λ(t) + Λ(x r ) + • • • + Λ(x 1 )] u u! e -[(N -r)Λ(t)+Λ(x r )+•••+Λ(x 1 )] r i=1 λ(x i ) dx 1 • • • dx r . = N r p r r! t 0 x r 0 • • • x 2 0 ∞ u=0 [(1 -p)[(N -r)Λ(t) + Λ(x r ) + • • • + Λ(x 1 )] u u! e -[(N -r)Λ(t)+Λ(xr)+•••+Λ(x 1 )] r i=1 λ(x i ) dx 1 • • • dx r . = N r p r r! t 0 xr 0 • • • x 2 0 e (1-p)[(N -r)Λ(t)+Λ(x r )+•••+Λ(x 1 )] e -[(N -r)Λ(t)+Λ(x r )+•••+Λ(x 1 )] r i=1 λ(x i ) dx 1 • • • dx r . = N r p r r! t 0 x r 0 • • • x 2 0 e -p[(N -r)Λ(t)+Λ(x r )+•••+Λ(x 1 )] r i=1 λ(x i ) dx 1 • • • dx r . = N r p r r!e -p(N -r)Λ(t) t 0 xr 0 • • • x 2 0 e -pΛ(xr) λ(x r ) • • • e -pΛ(x 1 ) λ(x 1 ) dx 1 • • • dx r = N r p r r! e -p(N -r)Λ(t) (1 -e -pΛ(t) ) r 1 p r r! = N r e -pΛ(t)(N -r) (1 -e -pΛ(t) ) r

This is the same expression as (3). It follows immediately that E[R(t)]

= N (1 -e -pΛ(t) ).

Similarly

E[U (t)] = N r=0 ∞ u=0 uP u,r (t) = N r=0 N r p r r! t 0 x r 0 • • • x 2 0 ∞ u=0 u (1 -p) u [(N -r)Λ(t) + Λ(x r ) + • • • + Λ(x 1 )] u u! •e [-(N -r)Λ(t)-Λ(x r )-•••-Λ(x 1 )] r i=1 λ(x i ) dx 1 • • • dx r = N r=0 N r p r r! t 0 x r 0 • • • x 2 0 (1 -p)[(N -r)Λ(t) + Λ(x r ) + • • • + Λ(x 1 )] e -p[(N -r)Λ(t)+Λ(xr)•••+Λ(x 1 )] r i=1 λ(x i ) dx 1 • • • dx r = N r=0 N r p r r! t 0 x r 0 • • • x 2 0 (1 -p) -∂ ∂p e -p[(N -r)Λ(t)+Λ(x r )+•••+Λ(x 1 )] r i=1 λ(x i ) dx 1 • • • dx r

A c c e p t e d m a n u s c r i p t

= - N r=0 N r p r r! (1 -p) ∂ ∂p t 0 x r 0 • • • x 2 0 e -p[(N -r)Λ(t)+Λ(xr)+•••+Λ(x 1 )] r i=1 λ(x i ) dx 1 • • • dx r = -(1 -p) N r=0 N r p r r! ∂ ∂p {e -p(N -r)Λ(t) (1 -e -pΛ(t) ) r 1 p r r! } = -(1 -p) N r=0 N r p r r! 1 p r r! ∂ ∂p {e -p(N -r)Λ(t) (1 -e -pΛ(t) ) r } - r p r+1 r! {e -p(N -r)Λ(t) (1 -e -pΛ(t) ) r } = -(1 -p) (4) ∂ ∂p { N r=0 N r e -p(N -r)Λ(t) (1 -e -pΛ(t) ) r } - N r=0 N r r p {e -p(N -r)Λ(t) (1 -e -pΛ(t) ) r } = N (1 -p) p (1 -e -pΛ(t) ) N -1 r-1=0 N -1 r -1 e -pΛ(t)(N -1-(r-1)) (1 -e -pΛ(t) ) r-1 = N (1 -p) p (1 -e -pΛ(t) ). ( 5
)
Given that

M (t) = R(t) + U (t), it is clear therefore that E[M (t)] = (N/p)(1 -e -pΛ(t) ).
3 Optimal Release Time for Testing Software.

In our cost models we will consider two distinct cases. In the first case we consider a fixed software life-cycle time t 0 where it is desired that if the software is released at time T , it functions well in the period (T, t 0 ]. In the second case we will want the software to function for some fixed mission time τ after release (that is to say, in the period (T, T + τ ]). Let us use c 1 to be the cost associated with encountering a fault (whether successfully repaired or not) during the testing period (0, T ], c 2 to be the cost of dealing with a fault encountered after the release

C(T ) = c 1 M (T ) + c 2 [M (g(T )) -M (T )] + c 3 T (6)
where g(T ) = g 1 (T ) = t 0 in case 1 (where we are considering a fixed life cycle time t 0 ), and g(T ) = g 2 (T ) = T + τ in case 2 (when we consider a fixed mission time τ ). The total expected cost for releasing the software at time T therefore takes the form

E[C(T )] = c 1 E[M (T )] + c 2 ( E[M (g(T ))] -E[M (T )] ) + c 3 T

A c c e p t e d m a n u s c r i p t

= c 2 N p (1 -e -pΛ(g(T )) ) -(c 2 -c 1 ) N p (1 -e -pΛ(T ) ) + c 3 T, (7)
and the objective is to find a T * which minimizes expected costs. Differentiating [START_REF] Goel | Software Error Detection model with applications[END_REF] with respect to T one finds critical points, and the T * minimizing expected costs is a solution of

λ(T ) e -pΛ(T ) = c 3 N (c 2 -c 1 ) , (8)
for case 1, while in case 2 it is a solution of

(c 2 -c 1 ) λ(T ) e -pΛ(T ) -c 2 λ(T + τ ) e -pΛ(T +τ ) = c 3 N . ( 9
)
If λ(t) (failure rate or discovery rate of faults) is non increasing and c 2 > c 1 , then any solution T 1 to equation ( 8) is unique. This gives a minimum cost for release if

λ(0)N [c 2 -c 1 ] > c 3 , and
otherwise the optimal value is T * = 0. Remember that in case 1 we must have T * ≤ g(T ) = t 0 , thus, denoting the optimal choice of T by T * ,

T * = min (t 0 , max(0, T 1 )) , ( 10
)
where T 1 is the solution of (8).

In the case of a fixed mission time τ, the critical point is given by the solution of

c 2 [λ(T ) e -pΛ(T ) -λ(T + τ ) e -pΛ(T +τ ) ] -c 1 λ(T ) e -pΛ(T ) = c 3 N (11)
In the case of a non increasing failure rate λ(t) the solution gives a unique minimum if

1 - c 1 c 2 > p[λ(T + τ )] 2 -λ (T + τ ) p[λ(T )] 2 -λ (T ) e -p[Λ(T +τ )-Λ(T )] (12)

Numerical Examples

We now apply our general results to specific failure rate models. As our first example we take λ(t) = be -bt . Then Λ(t) = 1 -e -bt . We illustrate the example with values of N and b suggested by the failure data set DS1 discussed by [START_REF] Goel | Software Error Detection model with applications[END_REF]. This data gives rise to estimates of Table 1 gives the optimal testing time T * for different values of p (the probability of perfect fault repair), for the cases where the fixed life cycle time is t 0 = 100 and for the situation when we are interested in the software functioning well during a mission time of τ = 2, 5, 10, 20, 50, 100. For our second example we take λ(t) = b/(1 + bt). Then Λ(t) = log(1 + b t). For the case of a fixed life cycle time t 0 , the optimal testing time T * is given by min(t 0 , max(0, T 1 )), where

T 1 = 1 b   N b(c 1 -c 2 ) c 3 1 1+p -1  
For the case of a fixed mission time τ, T * = max(0, T 2 ) where T 2 is given as the solution of

(c 2 -c 1 )b (1 + bT ) 1+p - c 2 b (1 + b(T + τ )) 1+p = c 3 N . ( 13
)
The optimal testing times T * and the corresponding minimal expected costs are given in Tables 3 and4.

p t 0 = 100 τ = 2 τ = 5 τ = 10 τ = 20 τ = 50 τ =

A c c e p t e d m a n u s c r i p t

  necessarily constant. A study is made of the resulting optimal release times for software under a cost structure for development.

time T , and c 3

 3 to be the cost of testing per unit time during the testing period. A common cost function model used for releasing the software at time T (see Okumoto and Goel (1980), McDaid and Wilson (2001), and Boland and Singh (2002)) takes the form

N

  = 1348 and b = 0.124. We also use the same cost parameters of c 1 = 1, c 2 = 5, c 3 = 100 suggested by Okumoto and Goel (1980) (and used by Boland and Singh (2002) and (2003) ) in order to make comparisons with other studies.

Figure 1

 1 Figure 1 illustrates the behavior of T * as a function of p plotted for a number of values of τ .

A c c e p t e d m a n u s c r i p t p t 0

 0 = 100 τ = 2 τ = 5 τ = 10 τ = 20 τ = 50 τ = 100 0.

Figure 1 :Figure 2 :Figure 3 :A c c e p t e d m a n u s c r i p tA c c e p t e d m a n u s c r i p t

 123 Figure 1: Optimal Testing Times for 3 Models with t 0 = 100 and b = 0.124.

Table 1 :

 1 Optimal Testing Time T * when λ(t) = 0.124 e -0.124 tThe corresponding minimum expected costs are given in the following table:

	2	13.99	0	6.13	10.69	13.17	13.98	13.99
	0.4	12.76	0	5.94	9.80	12.01	12.74	12.76
	0.6	11.63	0	5.73	9.02	10.96	11.61	11.63
	0.8	10.60	0.61	5.53	8.34	10.02	10.59	10.60
	1.0	9.68	1.10	5.33	7.74	9.18	9.67	9.68
	p t 0 = 100 τ = 2 τ = 5 τ = 10 τ = 20 τ = 50 τ = 100
	0.2	3414	1448 2575	3068	3328	3412	3414
	0.4	3161	1417 2387	2835	3080	3159	3161
	0.6	2928	1387 2218	2624	2852	2927	2928
	0.8	2714	1351 2067	2433	2643	2712	2714
	1.0	2517	1306 1932	2261	2452	2516	2517

Table 2 :

 2 Minimum

Expected Cost when λ(t) = 0.124 e -0.124 t

Table 3 :

 3 Optimal Testing Time T

	100

* when λ(t) = 0.124/(1 + 0.124 t)

We compare these two examples with the results obtained in [START_REF] Boland | Cost Implications of Imperfect Repair in Software Reliability[END_REF]. In that case λ(t) = b. The results are tabulated in Tables 5 and6.


Thorsten Schmidt

Winfried Stute
email: winfried.stute@math.uni-giessen.de

Shot-noise processes and the minimal martingale measure

Keywords: Minimal martingale measure, shot-noise process, jump diffusion Germany

This article proposes a model for stock prices which incorporates shot-noise effects. This means, that sudden jumps in the stock price are allowed, but their effect may decline as time passes by. Our model is general enough to capture arbitrary effects of this type. Generalizing previous approaches to shot-noise we in particular allow the decay to be stochastic. This model describes an incomplete market, so that the martingale measure is not unique. We derive the minimal martingale measure via continuous time methods.

Introduction

In financial markets information often comes as a surprise. This usually leads to a jump in stock prices, may it be upward or downward. Quite often this is an over-reaction and the effect partially vanishes as time passes by.

A by now classical approach to incorporate jumps in a stock price model is via so-called jump-diffusion, which originated from the work of Merton (1973). In a jump-diffusion model the stock prices may jump to a new level and then follow a geometric Brownian motion. This means that the jump effect persists in the model. In reality this need not be the case since, e.g., an upward jump may lead to profit taking. As a result, the jump effect may partially fade away.

In [START_REF] Arai | Minimal martingale measures for jump diffusion processes[END_REF] a multidimensional jump-diffusion model has been considered and the minimal martingale measure was computed. This paper generalizes the one-dimensional case to shot-noise models in several aspects.

The setup of the paper is as follows. In section two we extend the model of [START_REF] Altmann | A shot noise model for financial assets[END_REF] by considering a general model for shot-noise processes. We analyze the model in continuous time and derive the minimal martingale measure. Martingale measures are important for pricing derivatives of financial assets. Actually, these prices may be obtained as expectations of the payoffs at maturity w.r.t. a martingale measure. Finally we give conditions which ensure that the minimal martingale measure is indeed a probability measure. A small simulation study will exhibit typical features of the model.

Setup

Consider a sequence of iid random variables (rv) η i , i = 1, 2, . . . with values in R d and a sequence of independent Brownian motions B i , i = 1, 2, . . . . Let (λ t ) t≥0 be a positive hazard process, such that P( T 0 λ t dt < ∞) = 1 for a fixed horizon T . Finally, let Ñ be a standard Poisson process independent of all η i , B i and λ. Define Λ t := t 0 λ s ds and set

N t := ÑΛ t .
Then N is a Cox process with intensity λ. Denote its jump times by τ i .

Let a : R

+ × R × R d → R, b : R + × R × R d → R + be
J i (t) = J i (0) + t 0 a(s, J i (s), η i ) ds + t 0 b(s, J i (s), η i ) dB i (s), (1)
where J i (0) are i.i.d. real valued rv. Then we call

Y t := N t i=1 J i (t -τ i ) (2)
a general shot-noise process. To shorten notation, we write a i t := a(t-τ i , J i (t-τ i ), η i ) and

b i t := b(t-τ i , J i (t -τ i ), η i
). Note that we do not require that J i (0) and η i are independent.

It will turn out to be useful at some time to have η equal to J i (0).

We first mention some important special cases of (1):

1. The processes J i are of the form U i •h(t), with iid U i and a deterministic differentiable function h. We assume h(0) = 1. Then

J i (t) = U i + t 0 U i h (s) ds.
Indeed, this is a special case of (1) with b ≡ 0, η i = U i = J i (0) and a(s, J, η i ) = η i h (s). An important special example for h is h(t) = e -at , in which case the process is Markovian.

2. In the previous example the decay parameter a was fixed. As an extension, we may want a to be random. For this, consider a differentiable function h :

R + × R → R, with h(0) = 1 and J i (t) = U i • h(t, ηi ), with iid ηi . Assume J i equals J i (t) = U i + t 0 U i h (s, ηi ) ds.
This time we take η i = (U i , ηi ) and the choice of a, b is obvious. We would typically assume that the investor observes η i at the ith jump τ i of N . The decline is random, but if the jump occurs, the future decay is known. Some kind of regime-switching could be incorporated by letting h(t, x) = 1 {t∈[0,x]} . In this case the above setting

A c c e p t e d m a n u s c r i p t

would imply that at the upward-jump it is also immediately known, when the process jumps down again.

3. Up to now, b was always 0. In the next example we consider a stochastic noise effect. For example, if some unexpected event occurs investors may overreact due to various reasons. As time passes by, further information my arise so that the market comes back to more reasonable levels. To incorporate this effect we assume a mean reversion of the noise effect. Let J i be equal to

J i (t) = J i (0) + t 0 κ ηi -J i (s) ds + t 0 σdB i (s).
Intuitively, this means that the process Y jumps up by J i (0) at τ i and then diffuses to the level Y (τ i -) + η i . Also the volatility after the jump will decline exponentially such that the effect of coming back is captured. Other examples are immediately found: if J i is required to be nonnegative, one may set b(t, J, η) = σ √ J so that

dJ i (t) = κ(θ -J i (t))dt + σ J i (t)dB i (t). ( 3
)
The process specified by ( 3) is a so-called [START_REF] Cox | A theory of the term structure of interest rates[END_REF] process.

For more details we refer to [START_REF] Björk | Arbitrage Theory in Continuous Time[END_REF].

The main tool for obtaining the minimal martingale measure is the semimartingale

representation of Y . The increment of Y equals Y t+∆ -Y t = Nt i=1 J i (t + ∆ -τ i ) -J i (t -τ i ) + N t+∆ i=Nt+1 J i (t + ∆ -τ i ) -J i (t -τ i ) .
Letting ∆ go to zero one obtains that

dY t = N t- i=1 dJ i (t -τ i ) + d N t i=1 J i (0) .

A c c e p t e d m a n u s c r i p t

Denoting m 1 := E(J 1 (0)) we obtain the semimartingale decomposition:

dY t = N t- i=1 a i t dt + N t- i=1 b i t dB i t + d N t i=1 J i (0) (4) = N t- i=1 a i t + λ t m 1 dt + N t- i=1 b i t dB i t + d N t i=1 J i (0) - t 0 λ s m 1 ds ,
where the last two terms are local martingales.

The Model

In this section we consider a model for stock prices which utilizes the previously introduced general shot-noise process. First we show how to add shot-noise effects to a standard model for stock prices. Assume that S is a stochastic process for stock prices. For example, S could be a geometric Brownian motion. Set

S t = St • exp Y t . (5)
-insert figure 1 here -

We give some examples in Figure 1. Both simulations start from the model given by (5). The left graph uses a shot-noise component of the form J i (t) = U i exp(-at). After the jumps (at t = 0.5 and t = 0.9) the model tends towards the pre-jump level. The right plot incorporates a shot-noise specification with a stochastic mean-reversion effect according to (3), where κ = 2, θ = 0.1 and σ = 0.4. After the jumps the process shows an increasing volatility which decreases after some time. This is due to the non-zero volatility of J i after the jumps.

Remark 3.1. The model presented in [START_REF] Altmann | A shot noise model for financial assets[END_REF] assumes

S t = S 0 exp µ - σ 2 2 t -σB t N t i=1 (1 + U i h(t -τ i )) = S 0 exp µ - σ 2 2 t -σB t exp Nt i=1 ln(1 + U i h(t -τ i )) .

A c c e p t e d m a n u s c r i p t

This is a special case of ( 5), with J i chosen as

J i (t) = ln(1 + U i ) + t 0 U i h (s) 1 + U i h(s)
ds

and ãt ≡ µ + σ 2 /2, bt ≡ σ.
A special case: exponential decay. The shot-noise process considered so far has a quite general behaviour which typically leads to non-Markovian processes. Considering

J i (t) = U i h(t)
there are important special cases where the process is Markovian. For example, the classical jump-diffusion is obtained setting h ≡ 1. More interestingly, the process is also Markovian if h(x + y) = h(x)h(y), i.e., h has the form e -ax . The main reason for this is that the shot-noise process then takes on the following form

S t = Nt i=1 U i h(t -τ i ) = h(t) Nt i=1 U i h(-τ i ).
This process behaves deterministically between the jump times, with velocity h and jumps at τ i . Moreover, as has been shown in [START_REF] Gaspar | Quadratic portfolio credit risk models with shot-noise effects[END_REF], the only h for which Markovianity holds, is h(t) = e -at .

The minimal martingale measure

We assume for simplicity that the market interest rate is zero. The minimal martingale measure Q as proposed in [START_REF] Föllmer | Hedging of contingent claims under incomplete information[END_REF] can be described by its density L T , where d Q := L T dP and P is the objective measure. The density is determined as follows. Assume that S has the semimartingale representation S t = A t + M t with a local martingale M and an increasing process A of bounded variation. Furthermore, assume there exists a process λ which satisfies

A t = t 0 λs-d M s . (6)

A c c e p t e d m a n u s c r i p t

Then the minimal martingale density is given by

L T = E - • 0 λs-dM s T , ( 7
)
where E(•) t denotes the Doléans-Dade exponential (see, for example, Protter ( 2004)).

Recall that S t = St exp(Y t ). Consider a Brownian motion B which is independent of Y and two positive processes (ã t ) t≥0 and ( bt ) t≥0 which are adapted to the natural filtration of B. We assume that St is the strong solution of

dS t = S t (ã t dt + bt dB t ).
Proposition 4.1. Set m e := E[exp(J 1 (0))-1] and m 2 := E[(exp(J i (0))-1) 2 ] and assume both are finite. The minimal martingale measure is given by the density L T in (7), where

λt = 1 S t- • ãt + N t- i=1 a i t + λ t m e + 1 2 N t- i=1 (b i t ) 2 ( bt ) 2 + N t- i=1 (b i t ) 2 + λ t m 2 .
Proof. First, we compute the semimartingale representation of S using Itô's formula. See Protter ( 2004) for a suitable version. From

d(exp(Y t )) = exp(Y t-) dY c t + 1 2 d Y c t + d N t i=1 exp(J i (0)) -1
we obtain by (4) and the independence of S and Y :

dS t = exp(Y t-) St ãt dt + bt dB t + St exp(Y t-) N t- i=1 a i t dt + N t- i=1 b i t dB i t + 1 2 N t- i=1 (b i t ) 2 dt + d N t i=1 exp(J i (0)) -1 = S t-ãt + N t- i=1 a i t + λ t m e + 1 2 N t- i=1 (b i t ) 2 dt (8) + bt dB t + N t- i=1 b i t dB i t + d Nt i=1 exp(J i (0)) -1 - t 0 λ s m e ds .
With the semimartingale representation S t = A t + M t at hand we now compute M .

Using the independence of B and B i we obtain that

d M t = S 2 t-( bt ) 2 + N t- i=1 (b i t ) 2 + λ t m 2 dt. (9)
Using ( 8) and ( 9) in ( 6) we get the stated representation of λ.

A c c e p t e d m a n u s c r i p t

Once λ is obtained the density follows immediately:

L T = N t i=1 1 -λ τ i ∆M τ i exp . . . = N t i=1 1 -λ τ i S τ i -e J i (0) -1 exp . . . = Nt i=1 1 -l τ i e J i (0) -1 exp . . . , (10)
where we set

l t := ãt + N t- i=1 a i t + 1 2 (b i t ) 2 + λ t m e ( bt ) 2 + N t- i=1 (b i t ) 2 + λ t m 2 . ( 11
)
Conditions for Q to be a probability measure. It is well known, that the minimal martingale measure need not be a probability measure. It will turn out that this is often the case for the considered model. On the other side, there are special cases which guarantee that the minimal martingale measure is a probability measure. Note that the denominator in ( 11) is always positive. For example, consider J i t = U i , i.e., the case of a jump diffusion.

In this case the numerator equals ãt + λ t m e and so Q is a probability measure if

ãt + λ t m e ( bt ) 2 + λ t m 2 e J i (0) -1 < 1. ( 12
)
This condition corresponds to Assumption 3.1 from [START_REF] Arai | Minimal martingale measures for jump diffusion processes[END_REF]. One possibility to guarantee positivity of L in a jump-diffusion model is therefore to assume that J i < J and choose the parameters such that the above inequality ( 12) is satisfied.

In the general case we have the following result.

Proposition 4.2. Assume that J i (0) ≤ J with J > 0 and set K := 1/(exp( J) -1). Then L T > 0 a.s., if the following inequalities are satisfied for all i ≥ 1 and t ∈ [0, T ]:

a i t + (b i t ) 2 2 ≥ 0, ãt + λ t m e ≥ 0 ( 13
)
a i t + (b i t ) 2 1 2 -K < 0 (14) ãt + λ t m e -K bt + λ t m 2 < 0. ( 15
)
A c c e p t e d m a n u s c r i p t

Proof. First, K > 0. It follows from ( 14) that

N t- i=1 a i t + (b i t ) 2 1 2 -K < 0.
From ( 15) we obtain that

K bt + λ t m 2 -ãt -λ t m e > 0 > N t- i=1 a i t + (b i t ) 2 1 2 -K and therefore ãt + λ t m e + N t- i=1 a i t + (b i t ) 2 2 < K bt + λ t m 2 + N t- i=1 (b i t ) 2 .
Comparing the above inequality with (11), we obtain l t < K. Hence, if J i (0) ≥ 0 we have that

l τ i exp(J i (0)) -1 < 1. ( 16
)
On the other side, consider the case J i (0) < 0. First, (13) ensures that l t ≥ 0. Together with exp(J i (0)) -1 < 0 we have that l τ i exp(J i (0)) -1 ≤ 0 < 1.

Inserting ( 16) in (10) guarantees that L T > 0, and the proof is completed.

Stochastic noise This paragraph describes an example with stochastic noise effect, i.e., b i t is not zero. Furthermore the model is constructed such that the minimal martingale measure is an equivalent probability measure. Assume that J i t = η i -I i t , where each I i follows a Cox-Ingersoll-Ross process with mean reversion level zero:

dI i (t) = -κI i (t)dt + σ I i (t)dB i t ( 17
)
where κ, σ > 0 and random positive I i (0). This means that the process J jumps up at τ i by η i + I i (0) and then comes down to the level η i in a random way. We then have the following result.

Proposition 4.3. Assume that (15) as well as ãt + λ t m e ≥ 0 holds and, furthermore,

κ + σ 2 1 2 -K < 0. ( 18
)
A c c e p t e d m a n u s c r i p t

Then L T > 1 with probability 1.

Proof. First, it is well-known that zero is a reflecting boundary for I i . Second, observe that a i τ i +t = κ(η i -J i t ) ≥ 0 and (b i τ i +t ) 2 = σ 2 (η i -J i t ). Hence (13) holds. From (18) it follows that ( 14) holds. Applying Proposition 4.2 shows the assertion.

Example 4.4. Consider a Black-Scholes model for St , i.e., ãt = µ ≥ 0, and assume that λ t ≡ λ > 0. Then K(σ + λm 2 ) > µ + λm e implies (15). If additionally µ + λm e ≥ 0 and (18) holds, then Q is a probability measure.

Example 4.5. For σ = 0 and η i ≡ 0 we obtain the classical shot-noise process with negative jumps as a special case, in which J i t = η i -I i 0 exp(-κt) = -I i 0 exp(-κt). U i exp(-2t) with jumps at 0.4 and 0.9; right: shot-noise model with J i (t) according to (3).

1

  Winfried Stute, Department of Mathematics, University of Giessen, Arndtstr. 2, D-35392 Giessen,

  smooth functions, and let the processes J i , i = 1, 2, . . . be given by the unique strong solutions of the stochastic

Figure 1 :

 1 Figure 1: Two simulations of shot-noise processes according to (5). S is a geometric Brownian motion starting in 100 with σ = 0.4. Left: shot-noise model with J i (t) =