## Indicators

• Similars in SciELO

## On-line version ISSN 1807-0302

### Comput. Appl. Math. vol.31 no.2 São Carlos  2012

#### http://dx.doi.org/10.1590/S1807-03022012000200002

Solutions to the recurrence relation un+1 = vn+1 + un vn in terms of Bell polynomials

Christopher S. WithersI; Saralees NadarajahII, *

IApplied Mathematics Group, Industrial Research Limited, Lower Hutt, New Zealand E-mails: c.withers@irl.cri.nz
IISchool of Mathematics, University of Manchester, Manchester M13 9PL, UK E-mails: Saralees.Nadarajah@manchester.ac.uk

ABSTRACT

Motivated by time series analysis, we consider the problem of solving the recurrence relation un+1 = vn+1 + un vn for n 0 and un, given the sequence vn. A solution is given as a Bell polynomial. When vn can be written as a weighted sum of nth powers, then the solution un also takes this form.

Mathematical subject classification: 33E99.

Key words: autoregressive processes, Bell polynomials, convolution, maximum.

1 Introduction and summary

We define the convolution of sequences {an, bn : n > m} as

for n > 0. We consider the recurrence equation for un,

for n > 0, where u0 = v0, and vn is a given sequence.

The need for solutions to (1.1) arises with respect to the distribution of the maximum of first order autoregressive processes and more generally to that of autoregressive processes of any order, see Withers and Nadarajah (2010). There are many papers studying the distribution of the maximum of autoregressive processes, see, for example, Chernick and Davis (1982), McCormick and Mathew (1989), McCormick and Park (1992), Borkovec (2000), Ol'shanskii (2004) and Elek and Zempléni (2008). However, the results either give the limiting extreme value distributions or assume that the errors come from a specific class (for example, uniform distributed errors, negative binomial errors, ARCH(1) errors, etc). We are aware of no work giving the exact distribution of the maximum of autoregressive processes. As explained in Withers and Nadarajah (2010), solutions to (1.1) will lead to the exact distribution of the maximum of autoregressive processes of any order.

The aim of this note to provide solutions for (1.1) accessible for all scientists, not just mathematicians. These solutions are given in terms of Bell polynomials. In-built routines for Bell polynomials are available in most computer algebra packages. For example, see BellB in Mathematica and IncompleteBellPoly in Matlab. So, the solutions given will be accessible to most practitioners.

In Section 2, a solution for (1.1) is given as a Bell polynomial. In order to investigate the behavior of un for large n, an assumption needs to made on the behavior of vn for large n. In Section 3, we show that when vn can be written as a weighted sum of nth powers, then the solution un also has this form. This assumption is extended in Section 4 to the case when the weights are not constants, but polynomials in n. Some conclusions and future work are noted in Section 5.

2 The solution as a Bell polynomial

Theorem 2.1 provides an explicit solution of the recurrence relation (1.1) in terms of the complete ordinary Bell polynomial, n(w), a function of (w1, ..., wn), defined for any sequence w = (w1, w2, ...) by the formal generating function

where

The solution given by Theorem 2.1 can be computed by a single call to the in-built routine, BellB, in Mathematica or some other equivalent computer package. We believe that this is the most direct and the most efficient way to calculate a solution for the recurrence relation (1.1). However, in the absence of computer packages, the complete ordinary Bell polynomial can be calculated using the recurrence relation derived by Theorem 2.2.

Theorem 2.1. The recurrence relation (1.1) has the solution:

for n > 0, where wn = vn-1.

Proof. Multiply (1.1) by tn and sum from n = 0 to obtain

(U(t) – V (t)) /t = U(t) V(t),

where

since W(t) = tV(t), and

Taking the coefficient of tn in the last line gives the explicit solution (2.1).

Theorem 2.2. A recurrence relation for bn = n(w) is given by

bn = wn bn

for n > 1, where w0 = 0 and b0 = 1. For example,

 b1 = w0b1 + w1b0 = w1, b2 = w0b2 + w1b1 + w2b0 = w1b1 + w2 = + w2, b3 = w0b3 + w1b2 + w2b1 + w3b0 = w1b2 + w2b1 + w3 = + 2w1w2 + w3,

and so on.

Proof. Follows by taking the coefficient of tn in (1 – x)-1 -1 = x(1 - x)-1, where x = W(t).

An alternative to the complete ordinary Bell polynomial is the partial ordinary Bell polynomial, nj(w), defined by

for 0 < j < n. For example,

where δ00 = 1 and δn0 = 0 for n > 0. These polynomials are tabled on page 309 of Comtet (1974) for 1 < n < 10. Recurrence formulas for them are also given by Comtet (1974). They can be computed by a single call to the in-built routine, IncompleteBellPoly, in Matlab.

Theorem 2.3 states the relationship between the complete ordinary Bell polynomial and the partial ordinary Bell polynomial. It also provides a recurrence relation for the latter. Corollary 2.1 derives the relationship between un and vn.

Corollary 2.2 derives the reciprocal relationship between vn and un.

Theorem 2.3. We have

A recurrence relation for bnj = nj (w) is

for j1, j2 > 0. For example,

for j > 1.

Proof. Take the coefficient of tn in the expansion for (1 - W(t))-1 to obtain (2.3). The recurrence relation (2.4) follows by taking the coefficient of tn in . The proof is complete.

Corollary 2.1. We have

Proof. Applying (2.3) to (2.1), and reading the partial polynomials fromComtet's table, we obtain the results.

Corollary 2.2. We have vn = -n+1(x) for n > 0, where xn = –un-1.

Proof. This follows from W(t) = 1 - (1 - X(t))-1, where X(t) = –tU(t) = xntn. Alternatively, it follows by writing (1.1) as u'n+1 = v'n+1+u'n v'n, n > 0, u'0 = v'0, where u'n = -vn, v'n = -un and applying (2.1).

3 The weighted sum of powers solution

The solution (2.1) gives no indication of the behavior of un for large n. To obtain this we need to make an assumption on the behavior of vn for large n. Then if V(t) is tractable, we can obtain un-1, n > 1, as the coefficient of tn in (2.2). Theorems 3.1 and 3.2 illustrate this by two methods.

Theorem 3.1. Suppose vn is a weighted sum of powers, at least for large enough n, say

for n > n0, where 1 < r < , n0 > 0. Suppose also {δj} are the roots of

where pn+1(δ) = δn+1 -vn⊗ δn. Assume that {bj} are all non-zero and that {νj} are all non-zero and distinct. Assume also that the r roots of (3.2) are all distinct. Then un has the form

for n > m0, where J < I' = r + n0, m0 = 2n0 -1 if n0 > 1 and m0 = 0 if n0 = 0.

Proof. Having found {δj}, {γj} are the roots of

for ν = ν1, ..., νI, where Ajn(ν) = δjn/(δj-ν) and qn+1(ν) = νn+1 + un ⊗ δn. Note (3.4) can be written

where (An)kj = Ajn(νk), Qn = (Qn1, ..., QnI)' and Qnk = qn(νk). So, if J = r, a solution is

If r = , numerical solutions can be found by truncating the infinite matrix An and infinite vectors (Qn, γ) to N × N matrix and N-vectors, then increasing N until the desired precision is reached. The proof, which is by substitution, assumes that {δj,νj} are all distinct. The proof relies on the fact that if = 0 for 1 < n < J and r1,..., rJ are distinct, then a1 = ... = aJ = 0, since det( : 1 < n, j < J) 0.

If J < r, a solution is given by dropping r - J rows of (5). If J > r, there are not enough equations for a solution by this method. If n0 > 2, the values of un for n < 2n0 - 2 can be found from (1.1) or (2.1) or the extension of (2.5).

Now suppose that n0 > 1 and n = 2n0 - 1. Then s2 = 0 so that the result remains true.

Set f(δ) equal to the left and right hand sides of (3.2). Then a sufficient condition that its roots are distinct, is that f(δ) = 0 implies (δ) = 0.

The second and more general method is to compute un from its generating function via the generating function of vn and (2.2). The advantages of this method are: (i) it always applies, provided that the generating functions exist near t = 0; (ii) we shall see that it becomes clear how to extend the method when multiple roots exist.

Theorem 3.2. Suppose (3.1) holds. Then un has the form

for some R, Tj and pj(n), a polynomial of some degree nj.

Proof. Again we begin with (3.1). The generating function is V(t) = V0 + V1, where V0 = vntn, V1 = /(1 - sj), sj = νj t. Set D = (1 - sj), N = DV1. Then D(1 - tV(t)) = L, where L = D(1 - tV0) - tN is a polynomial of degree J = r + n0 > r, assuming that n0 > 0. So, we can write L = (1 - ttk) say, and

1 + tU(t) = (1 - tV(t))-1 = D/L.

Suppose first that {tk} are distinct. Then by the usual partial fractions expansion,

where qk(t) = D/Lk and Lk = L/(1 - ttk). So,

for n > 1. If n0 = 0 then V0 = 0, U(t) = N/L = mk(tk)/(1 - ttk), where mk(t) = M/Lk so that

for n > 0. Now suppose first that {tk} are not distinct. Then we can write L = , where nj = J and {Tk} are distinct. By the general partial fraction expansion, (compare Section 2.10 of Gradshteyn and Ryzhik (2007)),

where cj, nj-k+1 = /(k - 1)! and Qj(t) = (1 - tTj)njD/L. So, (3.6) holds, where

for n > 1. So, pj(n) is a polynomial of degree nj. If n0 = 0 just replace D/L by N/L as for the case of distinct roots; then un is equal to the right hand side of (3.6) for n > 0.

Corollaries 3.1 to 3.3 deduce the asymptotics of un and vn as n → ∞ when (3.1), (3.3) and (3.6) are satisfied. Further particular cases are considered by Corollaries 3.4 to 3.7.

Corollary 3.1. If (3.1) holds then

as n → ∞, where

Suppose in addition r = 2 and n0 = 1. By (3.1), V(t) = bj(1 - νjt)-1 = N(t)/D when |νjt| < 1, where N(t) = b1(1 - ν2t) + b2(1 - ν1t), D = diti for d0 = 1, d1 = -ν1 - ν2 and d2 = ν1ν2. So, 1 - tV(t) = M(t)/D, where M(t) = D - tN(t) = 1 - c1 + c2t2, c1 = ν1ν2 + b1 + b2 and c2 = ν1ν2 + b1ν2 + b2ν1.

Corollary 3.2. If (3.3) holds then

as n → ∞, where

Corollary 3.3. If (3.6) holds, set Tj = rj exp(iψj), r = rj and N = max{nj : rj = r}. Then

for

Corollary 3.4. If c2 0 and M(t) has two distinct roots then the solution (3.3) holds with J = 2.

Corollary 3.5. Suppose c2 = 0 c1. Then M(t) has one root and (3.3) holds with J = 1. Alternatively, since the right hand side of (2.2) is equal to D(1 - c1t)-1 - 1, we obtain un-1 = for n > 1.

Corollary 3.6. Suppose c1 = c2 = 0. Then M(t) = 1 and the right hand side of (2.2) is equal to D, giving u0 = d1, u1 = d2 and un = 0 for n > 2.

Corollary 3.7. Suppose c2 0 and M(t) has two equal roots, say t = t1. The root satisfies M(t) = (t) = 0, giving t1 = c1/(2c2) = 2/c1 and 4c2 = . So, M(t) = (1 - t/t1)2 and the right hand side of (2.2) is equal to

giving un-1 = {n + 1 + nd1t1 + (n - 1)d2} for n > 1. So, in this case, the solution is a weighted power with the weight linear in n.

4 An extension to polynomial weights

The weighted sum of powers assumption for vn arose naturally in Withers and Nadarajah (2010) assuming that a certain matrix had diagonal Jordan form, or at least if the eigenvalues of the non-diagonal Jordan blocks are zero. When this is not the case, we showed in Withers and Nadarajah (2010) that

for n > 1, so that v1 = wi0. For this more general case, the method of obtaining un from vn via its generating function, (2.2), still holds, as shown by Theorem 4.1.

Theorem 4.1. Suppose (4.1) holds. Then un has the form

for n > 1 for some J, ck and tk. Note (4.2) is of the form (3.3) with n0 = 1.

Proof. Setting s = νt and D = d/ds,

where

for |s| < 1. So, setting si = νit for |si| < 1, (4.1) gives the generating function for {vn}

V(t) = 1 + V0(t) + V1(t),

where

Let {θj, j = 1,..., R} be the distinct non-zero values of {νi, i = 1,..., r}. Set

n0 = max{mi : vi = 0}, J = M + 1 + n0, Mj = max {mi : vi = θj},

Then V0(t) is a polynomial of degree n0, Nj(t) is a polynomial of degree Mj, D(t) is a polynomial of degree M, and

where

say, Qj(t) is a polynomial of degree M - Mj, and N(t) is a polynomial of degree M. So,

say, is a polynomial of degree J, giving

Expanding in partial fractions (see Section 2.10 of Gradshteyn and Ryzhik (2007)) and taking the coefficient of tn gives un-1 as a weighted sum nth powers of {tk} with the weights polynomials in n. If {tk} are all distinct, then the partial fraction expansion has the form

so that (4.2) follows. If {tk} are not distinct, then proceed as for the case given in Section 3.

Suppose that m1 = 2, ν1 0, (mi,νi) = (1,0) for i > 1. Then R = 1, θ1 = ν1, M = M1 = m1 = 2, J = 4, M2 = m2 = 1 and vn = w10 + (n - 1)w11 for n > 2. So, setting s = ν1t and c0 = wi0,

V0(t) = c0t, V1(t) = w10t(1 - s)-1 + w11t2(1 - s)-2,

V(t) = v0 + c0t + N1(t)(1 - s)-2,

N(t) = N1(t) = w10t(1 - s) + w11t2, Q1(t) = 1,

D(t)(1 - tV(t)) = (1 - s)2(1 - tv0 - t2c0) - tN(t) = L(t) = (1 - tit).

So, if {ti} are distinct, then

say, giving (2) with J = 4. For analytic expressions for the roots of a quartic, see Section 3.8.3, page 17 of Abramowitz and Stegun (1964).

5 Conclusions

We have solved the recurrence relation un+1 = vn+1+ unvn, n 0 for un, given the sequence vn. The solution for un is in terms of Bell polynomials. This form is convenient since in-built routines for computing Bell polynomials are widely available. So, the solutions can be directly applied to derive the distribution of the maximum for autoregressive processes.

We have also established the behavior of un for large n by assuming some form for the behavior of vn for large n. The assumed forms are (3.1), a weighted sum of powers, and (4.1). As shown in Withers and Nadarajah (2010), one of these assumptions always holds. So, the assumptions are not at all restrictive.

The work presented can be extended in several ways: 1) consider solving un+1 = vn+1 + ωn+1 + unvn + un⊗ωn, n 0 for un, given the sequences vnand ωn; 2) consider solving un+1 = vn+1 + ωn+1 +µn+1 + unvn + un⊗ωn + un⊗µn, n 0 for un, given the sequences vn, ωn and µn; and, 3) consider solving multivariate forms of (1.1) taking the form

un+1 = vn+1 + unvn,

where the equality and the convolution operator are interpreted element wise. We hope to address some of these problems in a future paper.

Acknowledgments. The authors would like to thank the Editor, the Associate Editor and the referee for careful reading and for their comments which greatly improved the paper.

REFERENCES

[1] M. Abramowitz and I.A. Stegun, Handbook of Mathematical Functions. U.S. Department of Commerce, National Bureau of Standards, Applied Mathematics Series, 55 (1964).         [ Links ]

[2] M. Borkovec, Extremal behavior of the autoregressive process with ARCH(1) errors. Stochastic Processes and Their Applications, 85 (2000), 189-207.         [ Links ]

[3] M.R. Chernick and R.A. Davis, Extremes in autoregressive processes with uniform marginal distributions. Statistics and Probability Letters, 1 (1982), 85-88.         [ Links ]

[5] P. Elek and A. Zempléni, Tail behaviour and extremes of two-state Markov-switching autoregressive models. Computers and Mathematics with Applications, 55 (2008), 2839-2855.         [ Links ]

[6] I.S. Gradshteyn and I.M. Ryzhik, Tables of Integrals, Series and Products, seventh edition. Academic Press, New York (2007).         [ Links ]

[7] W.P. McCormick and G. Mathew, Asymptotic results for an extreme value estimator of the autocorrelation coefficient for a first order autoregressive sequence. In: Extreme Value Theory (Oberwolfach, 1987), pp. 166-180, Lecture Notes in Statistics, 51 (1989), Springer, New York.         [ Links ]

[8] W.P. McCormick and Y.S. Park, Asymptotic analysis of extremes from autoregressive negative binomial processes. Journal of Applied Probability, 29 (1992), 904-920.         [ Links ]

[9] K.A. Ol'shanskii, On the extremal index of a thinned autoregression process. (Russian) Vestnik Moskovskogo Universiteta. Seriya I. Matematika, Mekhanika, 70 (2004), 17-23. Translation in Moscow University Mathematics Bulletin, 59 (2004), 18-24.         [ Links ]

[10] C.S. Withers and S. Nadarajah, The distribution of the maximum of a first order autoregressive process: the continuous case. Metrika, doi: 10.1007/s00184-010-0301-0. (2010).         [ Links ]