Computation of the Frequency Function of a Quadratic Form in Random Normal Variables

One of the outstanding problems in the theory of time series analysis is the distribution problem in spectral analysis for small samples. When the number of observations is sufficiently large for the Central Limit Theorem to be applicable, the normal approximation can be used to advantage. In recent years, work has aimed at the discovery of more generally applicable approximate methods and of more rational criteria for the sample size at which the large sample theory becomes useful; the state of the art is summarized in two recent papers by Grenander, Pollak and Slepian [1] and Freiberger and Grenander [2]. Consider a sample <italic>x</italic> = (<italic>x</italic><subscrpt>1</subscrpt>, <italic>x</italic><subscrpt>2</subscrpt>, ··· , <italic>x<subscrpt>n</subscrpt></italic>) of successive values taken from a discrete-valued stationary time series ··· , <italic>y</italic><subscrpt>-1</subscrpt> , <italic>y</italic><subscrpt>0</subscrpt> , <italic>y</italic><subscrpt>1</subscrpt> , ··· , of normally distributed random variables with mean zero and covariance matrix <italic>R</italic> with elements <italic>r<subscrpt>&ngr;μ</subscrpt></italic> = <italic>E</italic>(<italic>x<subscrpt>&ngr;</subscrpt>x<subscrpt>μ</subscrpt></italic>). For stationarity we have: <italic>r<subscrpt>i,i</italic>+<italic>&ngr;</italic></subscrpt>≡ <italic>r<subscrpt>&ngr;</subscrpt></italic> = <italic>E</italic>(<italic>y<subscrpt>i</subscrpt>y<subscrpt>i</italic>+<italic>&ngr;</italic></subscrpt>), <italic>i, &ngr;</italic> = 0, ±1, ±2, ···. The Fourier-Stieltjes representation of the covariances is <italic>r<subscrpt>&ngr;</subscrpt></italic> = 1/2π ∫<supscrpt>π</supscrpt><subscrpt>π</subscrpt> <italic>e</italic><supscrpt><italic>i</italic>λ<italic>&ngr;</italic></supscrpt> <italic>dF</italic>(λ) (1) where the spectral distribution function <italic>F</italic>(<italic>λ</italic>) is bounded and nondecreasing and can, if it is absolutely continuous, be expressed in terms of a spectral density ƒ(λ): <italic>F</italic>(λ) = ∫<supscrpt>λ</supscrpt><subscrpt>π</subscrpt>ƒ(λ) (1) The practical determination of ƒ(λ) from observations of the process is effected by the introduction of a quadratic form <italic>Q</italic> = ∑<supscrpt><italic>n</italic></supscrpt><subscrpt><italic>&ngr;</italic>,λ=1</subscrpt><italic>w</italic><subscrpt><italic>&ngr;</italic>-<italic>μ</italic><italic></subscrpt>x<subscrpt>&ngr;</subscrpt>x<subscrpt>μ</subscrpt></italic> which is taken as an estimate for ƒ(λ); the coefficients <italic>w<subscrpt>&ngr;</subscrpt></italic> can be written in terms of a spectral weight function or “spectral window” <italic>w</italic>(λ): <italic>w<subscrpt>&ngr;</subscrpt></italic> = 1/2π ∫<supscrpt>π</supscrpt><subscrpt>-π</subscrpt> <italic>x</italic>(λ)<italic>e</italic><supscrpt><italic>i&ngr;</italic>λ</supscrpt> <italic>d</italic>λ. (4) For details, see the cited references. In order to obtain confidence limits for these estimates <italic>Q</italic> of the spectral density ƒ(λ), it is important to have methods for computing the distribution of <italic>Q</italic>. This paper deals with one such method which has proved very efficient for digital computer application. It is well known and shown, for instance, by H. Cramèr [3, p. 118], that the characteristic function of (3), where <italic>w<subscrpt>&ngr;μ</subscrpt></italic> ≡ <italic>w<subscrpt>&ngr;</italic>-<italic>μ</italic></subscrpt> are elements of a non-negative definite symmetric matrix <italic>W</italic>, is given by <italic>φ</italic>(<italic>z</italic>) = <italic>Ee<supscrpt>izQ</italic> = | <italic>I</italic> - 2<italic>izRW</italic> | <supscrpt>-1/2</supscrpt> = ∏<supscrpt>n</supscrpt><subscrpt><italic>j</italic>=1</subscrpt> (1 - 2<italic>iλ<subscrpt>j</subscrpt>z</italic>)<supscrpt>-1/2 (5) where the λ<subscrpt><italic>j</italic></subscrpt> are the <italic>n</italic> eigenvalues of the matrix product <italic>RW</italic>. Although <italic>RW</italic> is not necessarily a symmetric matrix, both <italic>R</italic> and <italic>W</italic> are symmetric and non-negative definite, so that all the λ<subscrpt><italic>j</italic></subscrpt> are real and non-negative. The frequency function <italic>g</italic>(<italic>x</italic>) of <italic>Q</italic> then follows from Fourier's inversion formula: <italic>g</italic>(<italic>x</italic>) = 1/2π ∫<supscrpt>∞</supscrpt><subscrpt>∞</subscrpt> <italic>e</italic></supscrpt>-<italic>ixz</italic></supscrpt> ∫<supscrpt><italic>n</italic></supscrpt><subscrpt><italic>j</italic>=1 (1 - 2<italic>iz</italic>λ<subscrpt><italic>j</italic></subscrpt>)<supscrpt>-1/2</supscrpt> <italic>dz</italic>. (6) Several ways have been suggested for evaluating <italic>g</italic> (<italic>x</italic>). Slepian [4] obtains a sum of finite integrals by deforming the contour of integration into a set of circles enclosing pairs of the branch points <italic>z<subscrpt>j</subscrpt></italic> = -<italic>i</italic>/(2λ<subscrpt><italic>j</italic></subscrpt>) of (6), and collapsing the circles. This method, which was also used in [2], works well when the eigenvalues cluster toward zero, but not otherwise. A method of repeated convolution of the frequency functions of its individual terms to obtain the frequency function of a quadratic form was developed in [5] and programmed for the IBM 650. It was found to be slowly convergent in most cases. Taking the logarithmic derivative of (5) and applying the inverse Fourier transform yields the following singular integral equation for <italic>g</italic>(<italic>x</italic>): <italic>xg</italic>(<italic>x</italic>) = ∫<italic><supscrpt>x</supscrpt><subscrpt>0</subscrpt> g</italic>(<italic>x</italic> - <italic>y</italic>)<italic>h</italic>(<italic>y</italic>) <italic>dy</italic> (7) where <italic>h</italic>(<italic>x</italic>) = 1/2 ∑<supscrpt><italic>n</italic></supscrpt><italic>&ngr;</italic>=1 <italic>e</italic><supscrpt>-<italic>i</italic>/2λ<subscrpt><italic>&ngr;</italic></subscrpt>. This observation forms the basis of an ingeneous method [1] for a computational scheme to derive the distribution of <italic>Q</italic>. The difficulties of obtaining initial values, associated with the high order zero of the frequency function <italic>g</italic>(<italic>x</italic>) at <italic>x</italic> = 0, are solved and the solution of (7) is discussed in detail in [1]. This method is suitable only for large-scale computers of the order of the IBM 704, and becomes, as does the method in [4], slowly convergent when the eigenvalues are densely spaced. Gurland [6] expanded the frequency function in terms of a Laguerre series and ([7]) presented a convergent series for the distribution fucntion using Laguerre polynomials. This method, which was found to converge slowly, formed the basis for the following procedure which has proved fast, accurate and reliable for a variety of problems. Essentially, the Laguerre expansion is now taken around Rice's approximation [8, p. 99], which is a type III distribution with appropriately chosen parameters. The following definition of the Laguerre polynomials <italic>L<subscrpt>n</subscrpt></italic><supscrpt>(α)</supscrpt> (<italic>x</italic>) = ∑<supscrpt><italic>n</italic></supscrpt><subscrpt><italic>&ngr;</italic>=0</subscrpt> (<italic>n</italic>+<italic>α</italic><italic>n</italic>-<italic>&ngr;</italic>)(-<italic>x)<supscrpt>&ngr;</supscrpt></italic>/<italic>&ngr;</italic>! (8) satisfies the orthogonality relations ∫<supscrpt>∞</supscrpt><subscrpt>0</subscrpt><italic>e</italic><supscrpt>-<italic>x</italic></supscrpt> <italic>x</italic><supscrpt><italic>α</italic></supscrpt><italic>L<subscrpt>n</subscrpt></italic><supscrpt>(<italic>α</italic>)</supscrpt> (<italic>x</italic>) <italic>L<subscrpt>m</subscrpt></italic><supscrpt>(<italic>α</italic>)</supscrpt> (<italic>x</italic>) <italic>dx</italic> = {Γ(<italic>α</italic> =1)(<italic>n</italic>+<italic>αn</italic>); <italic>m</italic> = <italic>n</italic><italic>m</italic> ≠ <italic>n</italic> (See Szegö [9]). Replacing <italic>x</italic> by λ<italic>x</italic> gives ∫<supscrpt>∞</supscrpt><subscrpt>0</subscrpt> <italic>e</italic><supscrpt>-λ<italic>x</italic></supscrpt><italic>x</italic><supscrpt><italic>α</italic></supscrpt><italic>L</italic><supscrpt>(<italic>α</italic></supscrpt>)<subscrpt><italic>n</italic></subscrpt>(Γ<italic>α</italic>) <italic>dx</italic>={Γ(<italic>α</italic>+1/λ<supscrpt><italic>α</italic>+1)</supscrpt> (<italic>n</italic>+<italic>αn</italic>);<italic>m</italic> = <italic>n</italic><italic>m</italic> ≠ <italic>n</italic> (9) which has the desired weight. The frequency function may now be expanded in a modified Laguerre series <italic>g</italic>(<italic>x</italic>) = <italic>Kx</italic><supscrpt><italic>α</italic></supscrpt><italic>e</italic><supscrpt>-λ<italic>x</italic></supscrpt>[<italic>c</italic><subscrpt>0</subscrpt> + <italic>c</italic><subscrpt>1</subscrpt><italic>L</italic><subscrpt>1</subscrpt><supscrpt>(<italic>α</italic>)</supscrpt> (λ<italic>x</italic>) + <italic>c</italic><subscrpt>2</subscrpt><italic>L</italic><subscrpt>2</subscrpt><supscrpt>(α)</supscrpt> (λ<italic>x</italic>) + …] (10) where <italic>K</italic> is chosen so that the weight integrates to one: ∫<supscrpt>∞</supscrpt><subscrpt>0</subscrpt><italic>Kx</italic><supscrpt>α</supscrpt><italic>e</italic><supscrpt>-λ<italic>x</italic></supscrpt> <italic>dx</italic> = 1, ∴ <italic>K</italic> = λ<supscrpt>α+1</supscrpt>/Γ(α + 1). (11) Multiplying both sides of (10) by <italic>L<subscrpt>n</subscrpt></italic><supscrpt>(α)</supscrpt> (λ<italic>x</italic>) and integrating from 0 to ∞, gives <italic>c<subscrpt>n</subscrpt></italic> = 1/(<italic>n</italic> + <italic>αn</italic>) ∫<supscrpt>∞</supscrpt><subscrpt>0</subscrpt><italic>L</italic><supscrpt>(<italic>α</italic>)<subscrpt><italic>n</italic></subscrpt>(λ<italic>x</italic>)<italic>g</italic>(<italic>x</italic>) <italic>dx</italic>. (12) Using (8), <italic>c<subscrpt>n</subscrpt></italic> = ∑<supscrpt><itali