Solutions to the recurrence relation u n+1 = v n+1 + u n ⊗ v n in terms of Bell polynomials

Solutions to the recurrence relation u_n+1 = v_n+1 + u_n ⊗ v_n in terms of Bell polynomials

Christopher S. Withers^I; Saralees Nadarajah^II

^IApplied Mathematics Group, Industrial Research Limited, Lower Hutt, New Zealand. E-mail: [e-mail]

^IISchool of Mathematics, University of Manchester, Manchester M13 9PL, UK. E-mail: [e-mail]

ABSTRACT

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

Mathematical subject classification: 33E99.

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

1 Introduction

We define the convolution of sequences {a_n, b_n: n > m} as

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

for n > 0, where u₀ = v₀, and v_n 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 autoregressiveprocesses, 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 thelimiting 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 distributionof the maximum of autoregressive processes. As explained in Withers andNadarajah (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 u_n for large n, an assumption needs to made on the behavior of v_n for large n. In Section 3, we show that when v_n can be written as a weighted sum of nth powers, then the solution u_n 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 notedin Section 5.

2 The solution as a Bell polynomial

Theorem 2.1 provides an explicit solution of the recurrence relation (1.1) interms of the complete ordinary Bell polynomial, (w), a function of (w₁, ..., w_n), defined for any sequence w = (w₁, w₂, ...) by the formal generatingfunction

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 w_n = v_n–1.

Proof. Multiply (1.1) by tⁿ 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 tⁿ in the last line gives the explicit solution (2.1).

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

b_n = w_n b_n

for n > 1, where w₀ = 0 and b₀ = 1. For example,

and so on.

Proof. Follows by taking the coefficient of tⁿ 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, _j(w), defined by

for 0 < j < n. For example, ₀(w) = δ_n0, ₁(w) = w_n, _n(w) = w₁ⁿ, where δ₀₀ = 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 u_n and v_n. Corollary 2.2 derives the reciprocal relationship between v_n and u_n.

Theorem 2.3. We have

A recurrence relation for b_nj = _j(w) is

for j₁, j₂> 0. For example,

for j > 1.

Proof. Take the coefficient of tⁿ in the expansion for (1 - W(t))^-1 to obtain (2.3). The recurrence relation (2.4) follows by taking the coefficient of tⁿ in W(t)^j1+j2 = W(t)^j1W(t)^j2. 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 v_n = –₊₁(x) for n > 0, where x_n = –u_n–1.

Proof. This follows from W(t) = 1-(1 - X(t))^-1, where

Alternatively, it follows by writing (1.1) as u'_n+1 = v'_n+1 + u'_n v'_n, n > 0, u'₀ = v'₀, where u'_n = -v_n, v'_n = -u_n and applying (2.1).

3 The weighted sum of powers solution

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

Theorem 3.1. Suppose v_n is a weighted sum of powers, at least for largeenough n, say

for n > n₀, where 1 < r < ∞, n₀> 0. Suppose also {δ_j} are the roots of

where p_n+1(δ) = δⁿ⁺¹-v_n δⁿ. Assume that {b_j} 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 u_n has the form

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

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

for ν = ν₁, ..., ν_I, where A_jn( ν) = δ_jⁿ/(δ_j – ν) and q_n+1( ν) = νⁿ⁺¹ + u_n δⁿ. Note (3.4) can be written

where (A_n)_kj = A_jn( ν_k), Q_n = (Q_n1, ..., Q_nI)' and Q_nk = q_n( ν_k). So, if J = r, a solution is

If r = ∞, numerical solutions can be found by truncating the infinite matrix A_n and infinite vectors (Q_n, γ) 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

and r₁,..., r_J are distinct, then a₁ = ... = a_J = 0, since det(r_jⁿ: 1 < n, j < J) ≠ 0.

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

Now suppose that n₀> 1 and n = 2n₀-1. Then s₂ = 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 u_n from its generating function via the generating function of v_n 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 u_n has the form

for some R, T_j and p_j(n), a polynomial of some degree n_j.

Proof. Again we begin with (3.1). The generating function is V(t) = V₀ + V₁, where V₀ = , V₁ = /(1 – s_j), s_j = ν_j t. Set D = (1 – s_j), N = DV₁. Then D(1 – tV(t)) = L, where L = D(1 – tV₀) – tN is a polynomial of degree J = r + n₀ > r, assuming that n₀ > 0. So, we can write L = (1 – tt_k) say, and

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

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

where q_k(t) = D/L_k and L_k = L/(1 – tt_k). So,

for n > 1. If n₀ = 0 then V₀ = 0, U(t) = N/L = m_k(t_k)/(1 – tt_k), where m_k(t) = M/L_k so that

for n > 0. Now suppose first that {t_k} are not distinct. Then we can write L = (1 – tT_j)^n_j, where n_j = J and {T_k} are distinct. By thegeneral partial fraction expansion, (compare Section 2.10 of Gradshteyn and Ryzhik (2007)),

where c_{j, nj–k+1} = /(k – 1)! and Q_j(t) = (1 – tT_j)^njD/L.So, (6) holds, where

for n > 1. So, p_j(n) is a polynomial of degree n_j. If n₀ = 0 just replace D/L by N/L as for the case of distinct roots; then u_n is equal to the right hand side of (6) for n > 0.

Corollaries 3.1 to 3.3 deduce the asymptotics of u_n and v_n 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 n₀ = 1. By (3.1), V(t) = b_j(1 – ν_jt)^–1 = N(t)/D when | ν_jt| < 1, where N(t) = b₁(1 – ν₂t) + b₂(1 – ν₁t), D = d_itⁱ for d₀ = 1, d₁ = – ν₁ – ν₂ and d₂ = ν₁ν₂. So, 1 – tV(t) = M(t)/D, where M(t) = D – tN(t) = 1 – c₁ + c₂t², c₁ = ν₁ν₂ + b₁ + b₂ and c₂ = ν₁ν₂ + b₁ν₂ + b₂ν₁.

Corollary 3.2. If (3.3) holds then

as n → ∞, where

Corollary 3.3. If (3.6) holds, set T_j = r_jexp(i ψ_j), r =

r_j and N = max{n_j : r_j= r}. Then

for

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

Corollary 3.5. Suppose c₂ = 0 ≠ c₁. 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 – c₁t)^–1 –1, we obtain u_n–1 = for n > 1.

Corollary 3.6. Suppose c₁ = c₂ = 0. Then M(t) = 1 and the right hand side of (2.2) is equal to D, giving u₀ = d₁, u₁ = d₂and u_n= 0 for n > 2.

Corollary 3.7. Suppose c₂≠ 0 and M(t) has two equal roots, say t = t₁. The root satisfies M(t) = (t) = 0, giving t₁ = c₁/(2c₂) = 2/c₁and 4c₂ = c₁². So, M(t) = (1 – t/t₁)² and the right hand side of (2.2) is equal to

giving u_n–1 = 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 v_n 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 v₁ = w_i0. For this more general case, the method of obtaining u_n from v_n via its generating function, (2.2), still holds, as shown by Theorem 4.1.

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

for n > 1 for some J, c_k and t_k. Note (4.2) is of the form (3.3) with n₀ = 1.

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

where

for |s| < 1. So, setting s_i = ν_i t for |s_i| < 1, (4.1) gives the generating function for {v_n}

V(t) = 1 + V₀(t) + V₁(t),

where

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

Then V₀(t) is a polynomial of degree n₀, N_j(t) is a polynomial of degree M_j, D(t) is a polynomial of degree M, and

where

say, Q_j(t) is a polynomial of degree M – M_j, and N(t) is a polynomial ofdegree 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 tⁿ gives u_n–1 as a weighted sum nth powers of {t_k} with the weights polynomials in n. If {t_k} are all distinct, then the partial fraction expansion has the form

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

Suppose that m₁ = 2, ν₁≠ 0, (m_i, ν_i) = (1, 0) for i > 1. Then

R = 1, θ₁ = ν₁, M = M₁ = m₁ = 2, J = 4, M₂ = m₂ = 1

and

So, setting s = ν₁t and c₀ = w_i0,

So, if {t_i} are distinct, then

say, giving (4.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 u_n+1 = v_n+1+u_n v_n, n ≠ 0 for u_n, given the sequence v_n. The solution for u_n 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 u_n for large n by assuming some form for the behavior of v_n 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 ofthese assumptions always holds. So, the assumptions are not at all restrictive.

The work presented can be extended in several ways: 1) consider solving u_n+1 = v_n+1 + ω_n+1 +u_n v_n + u_n ω_n, n ≠ 0 for u_n, given the sequences v_nand ω_n; 2) consider solving u_n+1 = v_n+1 + ω_n+1 + µ_n+1 +u_n v_n + u_n ω_n + u_n µ_n, n ≠ 0 for u_n, given the sequences v_n, ω_n and µ_n; and, 3) consider solving multivariate forms of (1.1) taking the form

u_n+1 = v_n+1 + u_n

v_n,

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.

Received: 21/VIII/10.

Accepted: 02/V/11.

#CAM-249/10.

• [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).
• [2] M. Borkovec, Extremal behavior of the autoregressive process with ARCH(1) errors. Stochastic Processes and Their Applications, 85 (2000), 189-207.
• [3] M.R. Chernick and R.A. Davis, Extremes in autoregressive processes with uniform marginal distributions. Statistics and Probability Letters, 1 (1982), 85-88.
• [4] L. Comtet, Advanced Combinatorics Reidel, Dordrecht (1974).
• [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.
• [6] I.S. Gradshteyn and I.M. Ryzhik, Tables of Integrals, Series and Products, seventh edition. Academic Press, New York (2007).
• [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.
• [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.
• [9] K.A. Ol'shanskii, On the extremal index of a thinned autoregression process.(Russian) Vestnik Moskovskogo Universiteta. Seriya I. Matematika, Mekhanika, 70, 17-23. Translation in "Moscow University Mathematics Bulletin", 59 (2004), 18-24.
• [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).

Publication Dates

• Publication in this collection
  06 Jan 2012
• Date of issue
  2011

History

• Accepted
  02 May 2011
• Received
  21 Aug 2010