g13dsf is a diagnostic checking routine suitable for use after fitting a vector ARMA model to a multivariate time series using g13ddf. The residual cross-correlation matrices are returned along with an estimate of their asymptotic standard errors and correlations. Also, g13dsf calculates the modified Li–McLeod portmanteau statistic and its significance level for testing model adequacy.

2 Specification

Fortran Interface

Subroutine g13dsf (

k, n, v, kmax, ip, iq, m, par, parhld, qq, ishow, r0, r, rcm, ldrcm, chi, idf, siglev, iw, liw, work, lwork, ifail)

Integer, Intent (In)	::	k, n, kmax, ip, iq, m, ishow, ldrcm, liw, lwork
Integer, Intent (Inout)	::	ifail
Integer, Intent (Out)	::	idf, iw(liw)
Real (Kind=nag_wp), Intent (In)	::	v(kmax,n), par((ip+iq)kk)
Real (Kind=nag_wp), Intent (Inout)	::	qq(kmax,k), r0(kmax,k), r(kmax,kmax,m), rcm(ldrcm,mkk)
Real (Kind=nag_wp), Intent (Out)	::	chi, siglev, work(lwork)
Logical, Intent (In)	::	parhld((ip+iq)kk)

C Header Interface

#include <nag.h>

void

g13dsf_ (const Integer *k, const Integer *n, const double v[], const Integer *kmax, const Integer *ip, const Integer *iq, const Integer *m, const double par[], const logical parhld[], double qq[], const Integer *ishow, double r0[], double r[], double rcm[], const Integer *ldrcm, double *chi, Integer *idf, double *siglev, Integer iw[], const Integer *liw, double work[], const Integer *lwork, Integer *ifail)

The routine may be called by the names g13dsf or nagf_tsa_multi_varma_diag.

3 Description

Let

W_{t} = {(w_{1 t}, w_{2 t}, \dots, w_{k t})}^{T}

, for

t = 1, 2, \dots, n

, denote a vector of

k

time series which is assumed to follow a multivariate ARMA model of the form

\begin{array}{l} W_{t} - μ = & ϕ_{1} (W_{t - 1} - μ) + ϕ_{2} (W_{t - 2} - μ) + \dots + ϕ_{p} (W_{t - p} - μ) \\ + ε_{t} - θ_{1} ε_{t - 1} - θ_{2} ε_{t - 2} - \dots - θ_{q} ε_{t - q}, \end{array}

(1)

where

ε_{t} = {(ε_{1 t}, ε_{2 t}, \dots, ε_{k t})}^{T}

, for

t = 1, 2, \dots, n

, is a vector of

k

residual series assumed to be Normally distributed with zero mean and positive definite covariance matrix

Σ

. The components of

ε_{t}

are assumed to be uncorrelated at non-simultaneous lags. The

ϕ_{i}

and

θ_{j}

are

k \times k

matrices of parameters.

{ϕ_{i}}

, for

i = 1, 2, \dots, p

, are called the autoregressive (AR) parameter matrices, and

{θ_{i}}

, for

i = 1, 2, \dots, q

, the moving average (MA) parameter matrices. The parameters in the model are thus the

p

(

k \times k

)

ϕ

-matrices, the

q

(

k \times k

)

θ

-matrices, the mean vector

μ

and the residual error covariance matrix

Σ

. Let

A (ϕ) = {[\begin{array}{l} ϕ_{1} & I & 0 & . & . & . & 0 \\ ϕ_{2} & 0 & I & 0 & . & . & 0 \\ . & . \\ . & . \\ . & . \\ ϕ_{p - 1} & 0 & . & . & . & 0 & I \\ ϕ_{p} & 0 & . & . & . & 0 & 0 \end{array}]}_{p k \times p k} and B (θ) = {[\begin{array}{l} θ_{1} & I & 0 & . & . & . & 0 \\ θ_{2} & 0 & I & 0 & . & . & 0 \\ . & . \\ . & . \\ . & . \\ θ_{q - 1} & 0 & . & . & . & I \\ θ_{q} & 0 & . & . & . & . & 0 \end{array}]}_{q k \times q k}

where

I

denotes the

k \times k

identity matrix.

The ARMA model (1) is said to be stationary if the eigenvalues of

A (ϕ)

lie inside the unit circle, and invertible if the eigenvalues of

B (θ)

lie inside the unit circle. The ARMA model is assumed to be both stationary and invertible. Note that some of the elements of the

ϕ

- and/or

θ

-matrices may have been fixed at pre-specified values (for example by calling g13ddf).

The estimated residual cross-correlation matrix at lag

l

is defined to the

k \times k

matrix

{\hat{R}}_{l}

whose

(i, j)

th element is computed as

{\hat{r}}_{i j} (l) = \frac{\sum_{t = l + 1}^{n} ({\hat{ε}}_{i t - l} - {\bar{ε}}_{i}) ({\hat{ε}}_{j t} - {\bar{ε}}_{j})}{\sqrt{\sum_{t = 1}^{n} {({\hat{ε}}_{i t} - {\bar{ε}}_{i})}^{2} \sum_{t = 1}^{n} {({\hat{ε}}_{j t} - {\bar{ε}}_{j})}^{2}}}, l = 0, 1, \dots, i ​ and ​ j = 1, 2, \dots, k,

where

{\hat{ε}}_{i t}

denotes an estimate of the

t

th residual for the

i

th series

ε_{i t}

and

{\bar{ε}}_{i} = \sum_{t = 1}^{n} {\hat{ε}}_{i t} / n

. (Note that

{\hat{R}}_{l}

is an estimate of

E (ε_{t - l} ε_{t}^{T})

, where

E

is the expected value.)

A modified portmanteau statistic,

Q_{(m)}^{*}

, is calculated from the formula (see Li and McLeod (1981))

Q_{(m)}^{*} = \frac{k^{2} m (m + 1)}{2 n} + n \sum_{l = 1}^{m} \hat{r} {(l)}^{T} ({\hat{R}}_{0}^{−1} \otimes {\hat{R}}_{0}^{−1}) \hat{r} (l),

where

\otimes

denotes Kronecker product,

{\hat{R}}_{0}

is the estimated residual cross-correlation matrix at lag zero and

\hat{r} (l) = vec ({\hat{R}}_{l}^{T})

, where

vec

of a

k \times k

matrix is a vector with the

(i, j)

th element in position

(i - 1) k + j

m

denotes the number of residual cross-correlation matrices computed. (Advice on the choice of

m

is given in Section 9.2.) Let

l_{C}

denote the total number of ‘free’ parameters in the ARMA model excluding the mean,

μ

, and the residual error covariance matrix

Σ

. Then, under the hypothesis of model adequacy,

Q_{(m)}^{*}

, has an asymptotic

χ^{2}

-distribution on

m k^{2} - l_{C}

degrees of freedom.

Let

\underset{̲}{\hat{r}} = (vec (R_{1}^{T}), vec (R_{2}^{T}), \dots, vec (R_{m}^{T}))

then the covariance matrix of

\underset{̲}{\hat{r}}

is given by

Var (\hat{\underset{̲}{r}}) = [Y - X {(X^{T} G G^{T} X)}^{−1} X^{T}] / n,

where

Y = I_{m} \otimes (Δ \otimes Δ)

and

G = I_{m} (G G^{T})

Δ

is the dispersion matrix

Σ

in correlation form and

G

a nonsingular

k \times k

matrix such that

G G^{T} = Δ^{- 1}

and

G Δ G^{T} = I_{k}

. The construction of the matrix

X

is discussed in Li and McLeod (1981). (Note that the mean,

μ

, plays no part in calculating

Var (\hat{\underset{̲}{r}})

and, therefore, is not required as input to g13dsf.)

4 References

Li W K and McLeod A I (1981) Distribution of the residual autocorrelations in multivariate ARMA time series models J. Roy. Statist. Soc. Ser. B 43 231–239

5 Arguments

The output quantities k, n, v, kmax, ip, iq, par, parhld and qq from g13ddf are suitable for input to g13dsf.

1: $k$ – Integer Input

On entry:

k

, the number of residual time series.

Constraint:

k \geq 1

2: $n$ – Integer Input

On entry:

n

, the number of observations in each residual series.

3: $v (kmax, n)$ – Real (Kind=nag_wp) array Input

On entry:

v (i, t)

must contain an estimate of the

i

th component of

ε_{t}

, for

i = 1, 2, \dots, k

and

t = 1, 2, \dots, n

Constraints:

no two rows of $v$ may be identical;
in each row there must be at least two distinct elements.

4: $kmax$ – Integer Input

On entry: the first dimension of the arrays v, qq, r and r0 and the second dimension of the array r as declared in the (sub)program from which g13dsf is called.

Constraint:

kmax \geq k

5: $ip$ – Integer Input

On entry:

p

, the number of AR parameter matrices.

Constraint:

ip \geq 0

6: $iq$ – Integer Input

On entry:

q

, the number of MA parameter matrices.

Constraint:

iq \geq 0

Note:

ip = iq = 0

is not permitted.

7: $m$ – Integer Input

On entry: the value of

m

, the number of residual cross-correlation matrices to be computed. See Section 9.2 for advice on the choice of m.

Constraint:

ip + iq < m < n

8: $par ((ip + iq) \times k \times k)$ – Real (Kind=nag_wp) array Input

On entry: the parameter estimates read in row by row in the order

ϕ_{1}, ϕ_{2}, \dots, ϕ_{p}

θ_{1}, θ_{2}, \dots, θ_{q}

Thus,

if $ip > 0$ , $par ((l - 1) \times k \times k + (i - 1) \times k + j)$ must be set equal to an estimate of the $(i, j)$ th element of $ϕ_{l}$ , for $l = 1, 2, \dots, p$ and $i = 1, 2, \dots, k$ ;
if $iq \geq 0$ , $par (p \times k \times k + (l - 1) \times k \times k + (i - 1) \times k + j)$ must be set equal to an estimate of the $(i, j)$ th element of $θ_{l}$ , for $l = 1, 2, \dots, q$ and $i = 1, 2, \dots, k$ .

The first

p \times k \times k

elements of par must satisfy the stationarity condition and the next

q \times k \times k

elements of par must satisfy the invertibility condition.

9: $parhld ((ip + iq) \times k \times k)$ – Logical array Input

On entry:

parhld (i)

must be set to .TRUE. if

par (i)

has been held constant at a pre-specified value and .FALSE. if

par (i)

is a free parameter, for

i = 1, 2, \dots, (p + q) \times k \times k

10: $qq (kmax, k)$ – Real (Kind=nag_wp) array Input/Output

On entry:

qq (i, j)

is an efficient estimate of the

(i, j)

th element of

Σ

. The lower triangle only is needed.

Constraint:

qq

must be positive definite.

On exit: if

ifail \neq 1

, then the upper triangle is set equal to the lower triangle.

11: $ishow$ – Integer Input

On entry: must be nonzero if the residual cross-correlation matrices

{{\hat{r}}_{i j} (l)}

and their standard errors

{se ({\hat{r}}_{i j} (l))}

, the modified portmanteau statistic with its significance and a summary table are to be printed. The summary table indicates which elements of the residual correlation matrices are significant at the

5 %

level in either a positive or negative direction; i.e., if

{\hat{r}}_{i j} (l) > 1.96 \times se ({\hat{r}}_{i j} (l))

then a ‘

+

’ is printed, if

{\hat{r}}_{i j} (l) < - 1.96 \times se ({\hat{r}}_{i j} (l))

then a ‘

-

’ is printed, otherwise a fullstop (.) is printed. The summary table is only printed if

k \leq 6

on entry.

The residual cross-correlation matrices, their standard errors and the modified portmanteau statistic with its significance are available also as output variables in r, rcm, chi, idf and siglev.

12: $r0 (kmax, k)$ – Real (Kind=nag_wp) array Output

On exit: if

i \neq j

, then

r0 (i, j)

contains an estimate of the

(i, j)

th element of the residual cross-correlation matrix at lag zero,

{\hat{R}}_{0}

. When

i = j

r0 (i, j)

contains the standard deviation of the

i

th residual series. If

ifail = 3

on exit then the first k rows and columns of r0 are set to zero.

13: $r (kmax, kmax, m)$ – Real (Kind=nag_wp) array Output

On exit:

r (l, i, j)

is an estimate of the

(i, j)

th element of the residual cross-correlation matrix at lag

l

, for

i = 1, 2, \dots, k

j = 1, 2, \dots, k

and

l = 1, 2, \dots, m

. If

ifail = 3

on exit then all elements of r are set to zero.

14: $rcm (ldrcm, m \times k \times k)$ – Real (Kind=nag_wp) array Output

On exit: the estimated standard errors and correlations of the elements in the array r. The correlation between

r (l, i, j)

and

r (l_{2}, i_{2}, j_{2})

is returned as

rcm (s, t)

where

s = (l - 1) \times k \times k + (j - 1) \times k + i

and

t = (l_{2} - 1) \times k \times k + (j_{2} - 1) \times k + i_{2}

except that if

s = t

, then

rcm (s, t)

contains the standard error of

r (l, i, j)

. If on exit,

ifail \geq 5

, then all off-diagonal elements of rcm are set to zero and all diagonal elements are set to

1 / \sqrt{n}

15: $ldrcm$ – Integer Input

On entry: the first dimension of the array rcm as declared in the (sub)program from which g13dsf is called.

Constraint:

ldrcm \geq m \times k \times k

16: $chi$ – Real (Kind=nag_wp) Output

On exit: the value of the modified portmanteau statistic,

Q_{(m)}^{*}

. If

ifail = 3

on exit then chi is returned as zero.

17: $idf$ – Integer Output

On exit: the number of degrees of freedom of chi.

18: $siglev$ – Real (Kind=nag_wp) Output

On exit: the significance level of chi based on idf degrees of freedom. If

ifail = 3

on exit, siglev is returned as one.

19: $iw (liw)$ – Integer array Workspace

20: $liw$ – Integer Input

On entry: the dimension of the array iw as declared in the (sub)program from which g13dsf is called.

Constraint:

liw \geq k \times \max (ip, iq)

21: $work (lwork)$ – Real (Kind=nag_wp) array Workspace

22: $lwork$ – Integer Input

On entry: the dimension of the array work as declared in the (sub)program from which g13dsf is called.

Constraint: if

lwork \geq k (n + kmax + 2) + m k^{2} (NPAR + m k^{2} + 1) + 3 k^{2} + (NPAR + 1) NPAR

NPAR = (p + q) k^{2}

23: $ifail$ – Integer Input/Output

On entry: ifail must be set to

0

−1

1

to set behaviour on detection of an error; these values have no effect when no error is detected.

A value of

0

causes the printing of an error message and program execution will be halted; otherwise program execution continues. A value of

−1

means that an error message is printed while a value of

1

means that it is not.

If halting is not appropriate, the value

−1

1

is recommended. If message printing is undesirable, then the value

1

is recommended. Otherwise, the value

−1

is recommended since useful values can be provided in some output arguments even when

ifail \neq 0

on exit. When the value $- 1$ or $1$ is used it is essential to test the value of ifail on exit.

On exit:

ifail = 0

unless the routine detects an error or a warning has been flagged (see Section 6).

6 Error Indicators and Warnings

If on entry

ifail = 0

−1

, explanatory error messages are output on the current error message unit (as defined by x04aaf).

Errors or warnings detected by the routine:

Note: in some cases g13dsf may return useful information.

$ifail = 1$: On entry, $ip = ⟨ value ⟩$ .
Constraint: $ip \geq 0$ .

On entry, $ip = 0$ and $iq = 0$ .
Constraint: $ip = iq = 0$ must not hold.

On entry, $iq = ⟨ value ⟩$ .
Constraint: $iq \geq 0$ .

On entry, $k = ⟨ value ⟩$ .
Constraint: $k \geq 1$ .

On entry, $kmax = ⟨ value ⟩$ and $k = ⟨ value ⟩$ .
Constraint: $kmax \geq k$ .

On entry, $ldrcm = ⟨ value ⟩$ , $m = ⟨ value ⟩$ and $k = ⟨ value ⟩$ .
Constraint: $ldrcm \geq m \times k \times k$ .

On entry, $liw = ⟨ value ⟩$ and the minimum size $required = ⟨ value ⟩$ .
Constraint: $liw \geq \max (ip, iq)$ .

On entry, $lwork = ⟨ value ⟩$ and the minimum size $required = ⟨ value ⟩$ .
Constraint: lwork must be sufficiently large.

On entry, $m = ⟨ value ⟩$ , $ip = ⟨ value ⟩$ and $iq = ⟨ value ⟩$ .
Constraint: $m > ip + iq$ .

On entry, $m = ⟨ value ⟩$ and $n = ⟨ value ⟩$ .
Constraint: $m < n$ .

$ifail = 2$: On entry, the AR parameter matrices are outside the stationarity region. To proceed you must supply different parameter estimates in the arrays par and qq.

On entry, the covariance matrix qq is not positive definite. To proceed you must supply different parameter estimates in the arrays par and qq.

On entry, the MA parameter matrices are outside the invertibility region. To proceed you must supply different parameter estimates in the arrays par and qq.

$ifail = 3$: On entry, at least one of the residual series in the array v has near-zero variance. In this case chi is set to zero, siglev to one and all the elements of r0 and r are set to zero.

On entry, at least two of the residual series are identical. In this case chi is set to zero, siglev to one and all the elements of r0 and r are set to zero.

$ifail = 4$: Excessive iterations needed to find zeros of determinental polynomials.

$ifail = 5$: On entry, the AR operator has a factor in common with the MA operator. To proceed you must either supply different parameter estimates in the array qq or delete this common factor from the model. In this case, the off-diagonal elements of rcm are returned as zero and the diagonal elements set to $1 / \sqrt{n}$ . All other output quantities will be correct.

$ifail = 6$: The matrix rcm could not be computed because one of its diagonal elements was found to be non-positive. In this case, the off-diagonal elements of rcm are returned as zero and the diagonal elements set to $1 / \sqrt{n}$ .

$ifail = - 99$: An unexpected error has been triggered by this routine. Please contact NAG.
See Section 7 in the Introduction to the NAG Library FL Interface for further information.

$ifail = - 399$: Your licence key may have expired or may not have been installed correctly.
See Section 8 in the Introduction to the NAG Library FL Interface for further information.

$ifail = - 999$: Dynamic memory allocation failed.
See Section 9 in the Introduction to the NAG Library FL Interface for further information.

7 Accuracy

The computations are believed to be stable.

8 Parallelism and Performance

g13dsf is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.

g13dsf makes calls to BLAS and/or LAPACK routines, which may be threaded within the vendor library used by this implementation. Consult the documentation for the vendor library for further information.

Please consult the X06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this routine. Please also consult the Users' Note for your implementation for any additional implementation-specific information.

9 Further Comments

9.1 Timing

The time taken by g13dsf depends upon the number of residual cross-correlation matrices to be computed,

m

, and the number of time series,

k

9.2 Choice of $m$

The number of residual cross-correlation matrices to be computed,

m

, should be chosen to ensure that when the ARMA model (1) is written as either an infinite order autoregressive process, i.e.,

W_{t} - μ = \sum_{j = 1}^{\infty} π_{j} (W_{t - j} - μ) + ε_{t}

or as an infinite order moving average process, i.e.,

W_{t} - μ = \sum_{j = 1}^{\infty} ψ_{j} ε_{t - j} + ε_{t}

then the two sequences of

k \times k

matrices

{π_{1}, π_{2}, \dots}

and

{ψ_{1}, ψ_{2}, \dots}

are such that

π_{j}

and

ψ_{j}

are approximately zero for

j > m

. An overestimate of

m

is, therefore, preferable to an under-estimate of

m

. In many instances the choice

m = 10

will suffice. In practice, to be on the safe side, you should try setting

m = 20

9.3 Checking a ‘White Noise’ Model

If you have fitted the ‘white noise’ model

W_{t} - μ = ε_{t}

then g13dsf should be entered with

p = 1

q = 0

, and the first

k^{2}

elements of par and parhld set to zero and .TRUE. respectively.

9.4 Approximate Standard Errors

When

ifail = 5

6

all the standard errors in rcm are set to

1 / \sqrt{n}

. This is the asymptotic standard error of

{\hat{r}}_{i j} (l)

when all the autoregressive and moving average parameters are assumed to be known rather than estimated.

9.5 Alternative Tests

{\hat{R}}_{0}

is useful in testing for instantaneous causality. If you wish to carry out a likelihood ratio test then the covariance matrix at lag zero

({\hat{C}}_{0})

can be used. It can be recovered from

{\hat{R}}_{0}

by setting

\begin{array}{l} {\hat{C}}_{0} (i, j) & = {\hat{R}}_{0} (i, j) \times {\hat{R}}_{0} (i, i) \times {\hat{R}}_{0} (j, j), & for ​ i \neq j \\ = {\hat{R}}_{0} (i, j) \times {\hat{R}}_{0} (i, j), & for ​ i = j \end{array}

10 Example

This example fits a bivariate AR(1) model to two series each of length

48

μ

has been estimated but

ϕ_{1} (2, 1)

has been constrained to be zero. Ten residual cross-correlation matrices are to be computed.

g13ds: FL CL CPP AD

NAG FL Interfaceg13dsf (multi_​varma_​diag)

▸▿ Contents

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

6 Error Indicators and Warnings

7 Accuracy

8 Parallelism and Performance

9 Further Comments

9.1 Timing

9.2 Choice of m

9.3 Checking a ‘White Noise’ Model

9.4 Approximate Standard Errors

9.5 Alternative Tests

10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results

NAG FL Interface
g13dsf (multi_varma_diag)

9.2 Choice of $m$