NAG Library Function Document

nag_tsa_varma_diagnostic (g13dsc)


    1  Purpose
    7  Accuracy


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


#include <nag.h>
#include <nagg13.h>
void  nag_tsa_varma_diagnostic (Integer k, Integer n, const double v[], Integer kmax, Integer ip, Integer iq, Integer m, const double par[], const Nag_Boolean parhld[], double qq[], Integer ishow, const char *outfile, double r0[], double r[], double rcm[], Integer pdrcm, double *chi, Integer *idf, double *siglev, NagError *fail)


Let Wt = w1t,w2t,,wktT , for t=1,2,,n, denote a vector of k time series which is assumed to follow a multivariate ARMA model of the form
Wt-μ= ϕ1Wt-1-μ+ϕ2Wt-2-μ++ϕpWt-p-μ +εt-θ1εt-1-θ2εt-2--θqεt-q, (1)
where εt = ε1t,ε2t,,εktT , for t=1,2,,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 by k matrices of parameters. ϕi, for i=1,2,,p, are called the autoregressive (AR) parameter matrices, and θi, for i=1,2,,q, the moving average (MA) parameter matrices. The parameters in the model are thus the p (k by k) ϕ-matrices, the q (k by k) θ-matrices, the mean vector μ and the residual error covariance matrix Σ. Let
Aϕ= ϕ1 I 0 . . . 0 ϕ2 0 I 0 . . 0 . . . . . . ϕp-1 0 . . . 0 I ϕp 0 . . . 0 0 pk×pk   and  Bθ= θ1 I 0 . . . 0 θ2 0 I 0 . . 0 . . . . . . θq-1 0 . . . I θq 0 . . . . 0 qk×qk  
where I denotes the k by 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 nag_tsa_varma_estimate (g13ddc)).
The estimated residual cross-correlation matrix at lag l is defined to the k by k matrix R^l whose i,jth element is computed as
r ^ i j l = t = l + 1 n ε ^ i t - l - ε - i ε ^ j t - ε - j t = 1 n ε ^ i t - ε - i 2 t = 1 n ε ^ j t - ε - j 2 ,   l =0,1,,i​ and ​j=1,2,,k ,  
where ε^it denotes an estimate of the tth residual for the ith series εit and ε-i=t=1nε^it/n. (Note that R^l is an estimate of Eεt-lεtT, where E is the expected value.)
A modified portmanteau statistic, Q m *, is calculated from the formula (see Li and McLeod (1981))
Qm * = k2 mm+1 2n + n l=1 m r^ lT R^ 0 -1 R^ 0 -1 r^ l ,  
where  denotes Kronecker product, R^0 is the estimated residual cross-correlation matrix at lag zero and r^l=vec R^lT , where vec of a k by k matrix is a vector with the i,jth element in position i-1k+j. m denotes the number of residual cross-correlation matrices computed. (Advice on the choice of m is given in Section 9.2.) Let lC 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 mk2-lC degrees of freedom.
Let r^̲=vec R1T ,vec R2T ,,vec RmT  then the covariance matrix of r^̲ is given by
where Y=ImΔΔ and G=ImGGT. Δ is the dispersion matrix Σ in correlation form and G a nonsingular k by k matrix such that GGT=Δ-1 and GΔGT=Ik. The construction of the matrix X is discussed in Li and McLeod (1981). (Note that the mean, μ, plays no part in calculating Varr̲^ and therefore is not required as input to nag_tsa_varma_diagnostic (g13dsc).)


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


The output quantities k, n, v, kmax, ip, iq, par, parhld and qq from nag_tsa_varma_estimate (g13ddc) are suitable for input to nag_tsa_varma_diagnostic (g13dsc).
1:     k IntegerInput
On entry: k, the number of residual time series.
Constraint: k1.
2:     n IntegerInput
On entry: n, the number of observations in each residual series.
3:     v[kmax×n] const doubleInput
On entry: v[kmax×t-1+i-1] must contain an estimate of the ith component of εt, for i=1,2,,k and t=1,2,,n.
  • no two rows of v may be identical;
  • in each row there must be at least two distinct elements.
4:     kmax IntegerInput
On entry: the first dimension of the arrays v, qq, r and r0 and the second dimension of the matrix r.
Constraint: kmaxk.
5:     ip IntegerInput
On entry: p, the number of AR parameter matrices.
Constraint: ip0.
6:     iq IntegerInput
On entry: q, the number of MA parameter matrices.
Constraint: iq0.
Note: ip=iq=0 is not permitted.
7:     m IntegerInput
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×k×k] const doubleInput
On entry: the parameter estimates read in row by row in the order ϕ1,ϕ2,,ϕp, θ1,θ2,,θq.
  • if ip>0, par[l-1×k×k+i-1×k+j-1] must be set equal to an estimate of the i,jth element of ϕl, for l=1,2,,p and i=1,2,,k;
  • if iq0, par[p×k×k+l-1×k×k+i-1×k+j-1] must be set equal to an estimate of the i,jth element of θl, for l=1,2,,q and i=1,2,,k.
The first p×k×k elements of par must satisfy the stationarity condition and the next q×k×k elements of par must satisfy the invertibility condition.
9:     parhld[ip+iq×k×k] const Nag_BooleanInput
On entry: parhld[i-1] must be set to Nag_TRUE if par[i-1] has been held constant at a pre-specified value and Nag_FALSE if par[i-1] is a free parameter, for i=1,2,,p+q×k×k.
10:   qq[kmax×k] doubleInput/Output
On entry: qq[kmax×j-1+i-1] is an efficient estimate of the i,jth element of Σ. The lower triangle only is needed.
Constraint: qq must be positive definite.
On exit: if fail.code= NE_G13D_AR, NE_G13D_DIAG, NE_G13D_FACT, NE_G13D_ITER, NE_G13D_MA, NE_G13D_RES, NE_G13D_ZERO_VAR or NE_NOT_POS_DEF, then the upper triangle is set equal to the lower triangle.
11:   ishow IntegerInput
On entry: must be nonzero if the residual cross-correlation matrices r^ijl and their standard errors ser^ijl, 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 r^ijl>1.96×ser^ijl then a ‘+’ is printed, if r^ijl<-1.96×ser^ijl then a ‘-’ is printed, otherwise a fullstop (.) is printed. The summary table is only printed if k6 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:   outfile const char *Input
On entry: the name of a file to which diagnostic output will be directed. If outfile is NULL the diagnostic output will be directed to standard output.
13:   r0[kmax×k] doubleOutput
On exit: if ij, then r0[kmax×j-1+i-1] contains an estimate of the i,jth element of the residual cross-correlation matrix at lag zero, R^0. When i=j, r0[kmax×j-1+i-1] contains the standard deviation of the ith residual series. If fail.code= NE_G13D_RES or NE_G13D_ZERO_VAR on exit then the first k rows and columns of r0 are set to zero.
14:   r[dim] doubleOutput
Note: the dimension, dim, of the array r must be at least kmax×kmax×m.
Where Rl,i,j appears in this document, it refers to the array element r[j-1×kmax×kmax+i-1×kmax+l-1].
On exit: Rl,i,j is an estimate of the i,jth element of the residual cross-correlation matrix at lag l, for i=1,2,,k, j=1,2,,k and l=1,2,,m. If fail.code= NE_G13D_RES or NE_G13D_ZERO_VAR on exit then all elements of r are set to zero.
15:   rcm[pdrcm×m×k×k] doubleOutput
On exit: the estimated standard errors and correlations of the elements in the array r. The correlation between Rl,i,j and Rl2,i2,j2 is returned as rcm[pdrcm×t+s] where s=l-1×k×k+j-1×k+i and t=l2-1×k×k+j2-1×k+i2 except that if s=t, then rcm[pdrcm×t+s] contains the standard error of Rl,i,j. If on exit, fail.code= NE_G13D_DIAG or NE_G13D_FACT, then all off-diagonal elements of rcm are set to zero and all diagonal elements are set to 1/n.
16:   pdrcm IntegerInput
On entry: the first dimension of the array rcm.
Constraint: pdrcmm×k×k.
17:   chi double *Output
On exit: the value of the modified portmanteau statistic, Q m *. If fail.code= NE_G13D_RES or NE_G13D_ZERO_VAR on exit then chi is returned as zero.
18:   idf Integer *Output
On exit: the number of degrees of freedom of chi.
19:   siglev double *Output
On exit: the significance level of chi based on idf degrees of freedom. If fail.code= NE_G13D_RES or NE_G13D_ZERO_VAR on exit, siglev is returned as one.
20:   fail NagError *Input/Output
The NAG error argument (see Section 3.7 in How to Use the NAG Library and its Documentation).

Error Indicators and Warnings

Dynamic memory allocation failed.
See Section in How to Use the NAG Library and its Documentation for further information.
On entry, argument value had an illegal value.
On entry, the AR parameter estimates are outside the stationarity region.
On entry, ip=0 and iq=0.
The matrix rcm could not be computed because one of its diagonal elements was found to be non-positive.
On entry, the AR operator has a factor in common with the MA operator.
Excessive iterations needed to find zeros of determinental polynomials.
On entry, the MA parameter matrices are outside the invertibility region.
On entry, at least two of the residual series are identical.
On entry, at least one of the residual series in the array v has near-zero variance.
On entry, ip=value.
Constraint: ip0.
On entry, iq=value.
Constraint: iq0.
On entry, k=value.
Constraint: k1.
On entry, kmax=value and k=value.
Constraint: kmaxk.
On entry, m=value and n=value.
Constraint: m<n.
On entry, m=value, ip=value and iq=value.
Constraint: m>ip+iq.
On entry, pdrcm=value, m=value and k=value.
Constraint: pdrcmm×k×k.
An internal error has occurred in this function. Check the function call and any array sizes. If the call is correct then please contact NAG for assistance.
See Section 2.7.6 in How to Use the NAG Library and its Documentation for further information.
Your licence key may have expired or may not have been installed correctly.
See Section 2.7.5 in How to Use the NAG Library and its Documentation for further information.
Cannot close file value.
On entry, the covariance matrix qq is not positive definite.
Cannot open file value for writing.


The computations are believed to be stable.

Parallelism and Performance

nag_tsa_varma_diagnostic (g13dsc) is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
nag_tsa_varma_diagnostic (g13dsc) 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 function. Please also consult the Users' Note for your implementation for any additional implementation-specific information.

Further Comments


The time taken by nag_tsa_varma_diagnostic (g13dsc) depends upon the number of residual cross-correlation matrices to be computed, m, and the number of time series, k.

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.,
or as an infinite order moving average process, i.e.,
Wt-μ=j= 1ψjεt-j+εt  
then the two sequences of k by k matrices π1,π2, and ψ1,ψ2, 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.

Checking a ‘White Noise’ Model

If you have fitted the ‘white noise’ model
then nag_tsa_varma_diagnostic (g13dsc) should be entered with p=1, q=0, and the first k2 elements of par and parhld set to zero and Nag_TRUE respectively.

Approximate Standard Errors

When fail.code= NE_G13D_DIAG or NE_G13D_FACT all the standard errors in rcm are set to 1/n. This is the asymptotic standard error of r^ijl when all the autoregressive and moving average parameters are assumed to be known rather than estimated.

Alternative Tests

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 C^0 can be used. It can be recovered from R^0 by setting
C^0i,j =R^0i,j×R^0i,i×R^0j,j, for ​ij =R^0i,j×R^0i,j, for ​i=j  


This example fits a bivariate AR(1) model to two series each of length 48. μ has been estimated but ϕ12,1 has been constrained to be zero. Ten residual cross-correlation matrices are to be computed.

Program Text

Program Text (g13dsce.c)

Program Data

Program Data (g13dsce.d)

Program Results

Program Results (g13dsce.r)

© The Numerical Algorithms Group Ltd, Oxford, UK. 2017