The function uses the method of Durbin (see Durbin (1960) and Golub and Van Loan (1996)). Optionally the mean square prediction errors and/or the partial correlation coefficients for each step can be returned.

References

Bunch J R (1985) Stability of methods for solving Toeplitz systems of equations SIAM J. Sci. Statist. Comput. 6 349–364

Bunch J R (1987) The weak and strong stability of algorithms in numerical linear algebra Linear Algebra Appl. 88/89 49–66

Cybenko G (1980) The numerical stability of the Levinson–Durbin algorithm for Toeplitz systems of equations SIAM J. Sci. Statist. Comput. 1 303–319

Durbin J (1960) The fitting of time series models Rev. Inst. Internat. Stat. 28 233

Golub G H and Van Loan C F (1996) Matrix Computations (3rd Edition) Johns Hopkins University Press, Baltimore

Parameters

Compulsory Input Parameters

1: $n$ – int64int32nag_int scalar: The order of the Toeplitz matrix $T$ .

Constraint: $n \geq 0$ . When $n = 0$ , then an immediate return is effected.
2: $t (0 : n)$ – double array: $t (0)$ must contain the value $τ_{0}$ of the diagonal elements of $T$ , and the remaining n elements of t must contain the elements of the vector $t$ .

Constraint: $t (0) > 0.0$ . Note that if this is not true, then the Toeplitz matrix cannot be positive definite.
3: $wantp$ – logical scalar: Must be set to true if the partial (auto)correlation coefficients are required, and must be set to false otherwise.
4: $wantv$ – logical scalar: Must be set to true if the mean square prediction errors are required, and must be set to false otherwise.

Optional Input Parameters

None.

Output Parameters

1: $x (n)$ – double array: The solution vector $x$ .
2: $p (:)$ – double array: The dimension of the array p will be $\max (1, n)$ if $wantp = true$ and $1$ otherwise

With wantp as true, the $i$ th element of p contains the partial (auto)correlation coefficient, or reflection coefficient, $p_{i}$ for the $i$ th step. (See Further Comments and Chapter G13.) If wantp is false, then p is not referenced. Note that in any case, $x_{n} = p_{n}$ .
3: $v (:)$ – double array: The dimension of the array v will be $\max (1, n)$ if $wantv = true$ and $1$ otherwise

With wantv as true, the $i$ th element of v contains the mean square prediction error, or predictor error variance ratio, $v_{i}$ , for the $i$ th step. (See Further Comments and Chapter G13.) If wantv is false, then v is not referenced.
4: $vlast$ – double scalar: The value of $v_{n}$ , the mean square prediction error for the final step.
5: $ifail$ – int64int32nag_int scalar: $ifail = 0$ unless the function detects an error (see Error Indicators and Warnings).

Error Indicators and Warnings

Note: nag_linsys_real_toeplitz_yule (f04fe) may return useful information for one or more of the following detected errors or warnings.

Errors or warnings detected by the function:

Cases prefixed with W are classified as warnings and do not generate an error of type NAG:error_n. See nag_issue_warnings.

$ifail = - 1$

On entry,	$n < 0$ ,
or	$t (0) \leq 0.0$ .

W $ifail > 0$

The principal minor of order (

ifail + 1

) of the Toeplitz matrix is not positive definite to working accuracy. If, on exit,

x_{ifail}

is close to unity, then the principal minor was close to being singular, and the sequence

τ_{0}, τ_{1}, \dots, τ_{ifail}

may be a valid sequence nevertheless. The first ifail elements of x return the solution of the equations

T_{ifail} x = - {(τ_{1}, τ_{2}, \dots, τ_{ifail})}^{T},

where

T_{ifail}

is the ifailth principal minor of

T

. Similarly, if wantp and/or wantv are true, then p and/or v return the first ifail elements of p and v respectively and vlast returns

v_{ifail}

. In particular if

ifail = n

, then the solution of the equations

T x = - t

is returned in x, but

τ_{n}

is such that

T_{n + 1}

would not be positive definite to working accuracy.

$ifail = - 99$: An unexpected error has been triggered by this routine. Please contact NAG.

$ifail = - 399$: Your licence key may have expired or may not have been installed correctly.

$ifail = - 999$: Dynamic memory allocation failed.

Accuracy

The computed solution of the equations certainly satisfies

r = T x + t,

where

{‖r‖}_{1}

is approximately bounded by

{‖r‖}_{1} \leq c ε (\prod_{i = 1}^{n} (1 + |p_{i}|) - 1),

c

being a modest function of

n

and

ε

being the machine precision. This bound is almost certainly pessimistic, but it has not yet been established whether or not the method of Durbin is backward stable. If

|p_{n}|

is close to one, then the Toeplitz matrix is probably ill-conditioned and hence only just positive definite. For further information on stability issues see Bunch (1985), Bunch (1987), Cybenko (1980) and Golub and Van Loan (1996). The following bounds on

{‖t^{- 1}‖}_{1}

hold:

\max (\frac{1}{v_{n - 1}}, \frac{1}{\prod_{i = 1}^{n - 1} (1 - p_{i})}) \leq {‖T^{- 1}‖}_{1} \leq \prod_{i = 1}^{n - 1} (\frac{1 + |p_{i}|}{1 - |p_{i}|}) .

Note:

v_{n} < v_{n - 1}

. The norm of

T^{- 1}

may also be estimated using function nag_linsys_real_gen_norm_rcomm (f04yd).

Further Comments

The number of floating-point operations used by nag_linsys_real_toeplitz_yule (f04fe) is approximately

2 n^{2}

, independent of the values of wantp and wantv.

The mean square prediction error,

v_{i}

, is defined as

v_{i} = (τ_{0} + (τ_{1} τ_{2} \dots τ_{i - 1}) y_{i - 1}) / τ_{0},

where

y_{i}

is the solution of the equations

T_{i} y_{i} = - {(τ_{1} τ_{2} \dots τ_{i})}^{T}

and the partial correlation coefficient,

p_{i}

, is defined as the

i

th element of

y_{i}

. Note that

v_{i} = (1 - p_{i}^{2}) v_{i - 1}

Example

This example finds the solution of the Yule–Walker equations

T x = - t

, where

T = (\begin{array}{l} 4 & 3 & 2 & 1 \\ 3 & 4 & 3 & 2 \\ 2 & 3 & 4 & 3 \\ 1 & 2 & 3 & 4 \end{array}) and t = (\begin{array}{l} 3 \\ 2 \\ 1 \\ 0 \end{array}) .

Open in the MATLAB editor: f04fe_example

function f04fe_example


fprintf('f04fe example results\n\n');

% Solve Yule-Walker equations for symmetric definite Toeplitz system
n = int64(4);
t = [4     3     2     1     0];

wantp = true;
wantv = true;
[x, p, v, vlast, ifail] = ...
  f04fe(n, t, wantp, wantv);

disp('Solution vector');
disp(x');
disp('Reflection coefficients');
disp(p');
disp('Mean square prediction errors');
disp(v');

f04fe example results

Solution vector
   -0.8000    0.0000   -0.0000    0.2000

Reflection coefficients
   -0.7500    0.1429    0.1667    0.2000

Mean square prediction errors
    0.4375    0.4286    0.4167    0.4000

PDF version (NAG web site, 64-bit version, 64-bit version)

Chapter Contents

Chapter Introduction

NAG Toolbox