NAG Library Function Document

nag_rand_varma (g05pjc)

+− 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

+− 10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results

1 Purpose

nag_rand_varma (g05pjc) generates a realization of a multivariate time series from a vector autoregressive moving average (VARMA) model. The realization may be continued or a new realization generated at subsequent calls to nag_rand_varma (g05pjc).

2 Specification

#include <nag.h>

#include <nagg05.h>

void	nag_rand_varma (Nag_OrderType order, Nag_ModeRNG mode, Integer n, Integer k, const double xmean[], Integer p, const double phi[], Integer q, const double theta[], const double var[], Integer pdv, double r[], Integer lr, Integer state[], double x[], Integer pdx, NagError *fail)

3 Description

Let the vector

X_{t} = {(x_{1 t}, x_{2 t}, \dots, x_{k t})}^{T}

, denote a

k

-dimensional time series which is assumed to follow a vector autoregressive moving average (VARMA) model of the form:

\begin{matrix} X_{t} - μ = & ϕ_{1} (X_{t - 1} - μ) + ϕ_{2} (X_{t - 2} - μ) + \dots + ϕ_{p} (X_{t - p} - μ) + \\ ε_{t} - θ_{1} ε_{t - 1} - θ_{2} ε_{t - 2} - \dots - θ_{q} ε_{t - q} \end{matrix}

(1)

where

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

, is a vector of

k

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

Σ

. The components of

ε_{t}

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

ϕ_{i}

's and

θ_{j}

's are

k

k

matrices of parameters.

\{ϕ_{i}\}

, for

i = 1, 2, \dots, p

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

\{θ_{j}\}

, for

j = 1, 2, \dots, q

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

p

k

k

ϕ

-matrices, the

q

k

k

θ

-matrices, the mean vector

μ

and the residual error covariance matrix

Σ

. Let

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

where

I

denotes the

k

k

identity matrix.

The model (1) must be both stationary and invertible. The model 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.

For

k \geq 6

the VARMA model (1) is recast into state space form and a realization of the state vector at time zero computed. For all other cases the function computes a realization of the pre-observed vectors

X_{0}, X_{- 1}, \dots, X_{1 - p}

ε_{0}, ε_{- 1}, \dots, ε_{1 - q}

, from (1), see Shea (1988). This realization is then used to generate a sequence of successive time series observations. Note that special action is taken for pure MA models, that is for

p = 0

At your request a new realization of the time series may be generated more efficiently using the information in a reference vector created during a previous call to nag_rand_varma (g05pjc). See the description of the argument mode in Section 5 for details.

The function returns a realization of

X_{1}, X_{2}, \dots, X_{n}

. On a successful exit, the recent history is updated and saved in the array r so that nag_rand_varma (g05pjc) may be called again to generate a realization of

X_{n + 1}, X_{n + 2}, \dots

, etc. See the description of the argument mode in Section 5 for details.

Further computational details are given in Shea (1988). Note, however, that nag_rand_varma (g05pjc) uses a spectral decomposition rather than a Cholesky factorization to generate the multivariate Normals. Although this method involves more multiplications than the Cholesky factorization method and is thus slightly slower it is more stable when faced with ill-conditioned covariance matrices. A method of assigning the AR and MA coefficient matrices so that the stationarity and invertibility conditions are satisfied is described in Barone (1987).

One of the initialization functions nag_rand_init_repeatable (g05kfc) (for a repeatable sequence if computed sequentially) or nag_rand_init_nonrepeatable (g05kgc) (for a non-repeatable sequence) must be called prior to the first call to nag_rand_varma (g05pjc).

4 References

Barone P (1987) A method for generating independent realisations of a multivariate normal stationary and invertible ARMA

(p, q)

process J. Time Ser. Anal. 8 125–130

Shea B L (1988) A note on the generation of independent realisations of a vector autoregressive moving average process J. Time Ser. Anal. 9 403–410

5 Arguments

1: order – Nag_OrderTypeInput

On entry: the order argument specifies the two-dimensional storage scheme being used, i.e., row-major ordering or column-major ordering. C language defined storage is specified by

order = Nag_RowMajor

. See Section 3.2.1.3 in the Essential Introduction for a more detailed explanation of the use of this argument.

Constraint:

order = Nag_RowMajor

Nag_ColMajor

2: mode – Nag_ModeRNGInput

On entry: a code for selecting the operation to be performed by the function.

$mode = Nag_InitializeReference$: Set up reference vector and compute a realization of the recent history.
$mode = Nag_GenerateFromReference$: Generate terms in the time series using reference vector set up in a prior call to nag_rand_varma (g05pjc).
$mode = Nag_InitializeAndGenerate$: Combine the operations of $mode = Nag_InitializeReference$ and $Nag_GenerateFromReference$ .
$mode = Nag_ReGenerateFromReference$: A new realization of the recent history is computed using information stored in the reference vector, and the following sequence of time series values are generated.

mode = Nag_GenerateFromReference

Nag_ReGenerateFromReference

, then you must ensure that the reference vector r and the values of k, p, q, xmean, phi, theta, var and pdv have not been changed between calls to nag_rand_varma (g05pjc).

Constraint:

mode = Nag_InitializeReference

Nag_GenerateFromReference

Nag_InitializeAndGenerate

Nag_ReGenerateFromReference

3: n – IntegerInput

On entry:

n

, the number of observations to be generated.

Constraint:

n \geq 0

4: k – IntegerInput

On entry:

k

, the dimension of the multivariate time series.

Constraint:

k \geq 1

5: xmean[k] – const doubleInput

On entry:

μ

, the vector of means of the multivariate time series.

6: p – IntegerInput

On entry:

p

, the number of autoregressive parameter matrices.

Constraint:

p \geq 0

7: phi[ $k \times k \times p$ ] – const doubleInput

On entry: must contain the elements of the

p \times k \times k

autoregressive parameter matrices of the model,

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

. The

(i, j)

th element of

ϕ_{l}

is stored in

phi [(l - 1) \times k \times k + (j - 1) \times k + i - 1]

, for

l = 1, 2, \dots, p

i = 1, 2, \dots, k

and

j = 1, 2, \dots, k

Constraint: the elements of phi must satisfy the stationarity condition.

8: q – IntegerInput

On entry:

q

, the number of moving average parameter matrices.

Constraint:

q \geq 0

9: theta[ $k \times k \times q$ ] – const doubleInput

On entry: must contain the elements of the

q \times k \times k

moving average parameter matrices of the model,

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

. The

(i, j)

th element of

θ_{l}

is stored in

theta [(l - 1) \times k \times k + (j - 1) \times k + i - 1]

, for

l = 1, 2, \dots, q

i = 1, 2, \dots, k

and

j = 1, 2, \dots, k

Constraint: the elements of theta must be within the invertibility region.

10: var[ $\dim$ ] – const doubleInput

Note: the dimension, dim, of the array var must be at least

pdv \times k

Where

VAR (i, j)

appears in this document, it refers to the array element

$var [(j - 1) \times pdv + i - 1]$ when $order = Nag_ColMajor$ ;
$var [(i - 1) \times pdv + j - 1]$ when $order = Nag_RowMajor$ .

On entry:

VAR (i, j)

must contain the (

i, j

)th element of

Σ

, for

i = 1, 2, \dots, k

and

j = 1, 2, \dots, k

. Only the lower triangle is required.

Constraint: the elements of var must be such that

Σ

is positive semidefinite.

11: pdv – IntegerInput

On entry: the stride separating row or column elements (depending on the value of order) in the array var.

Constraint:

pdv \geq k

12: r[lr] – doubleCommunication Array

On entry: if

mode = Nag_GenerateFromReference

Nag_ReGenerateFromReference

, the array r as output from the previous call to nag_rand_varma (g05pjc) must be input without any change.

mode = Nag_InitializeReference

Nag_InitializeAndGenerate

, the contents of r need not be set.

On exit: information required for any subsequent calls to the function with

mode = Nag_GenerateFromReference

Nag_ReGenerateFromReference

. See Section 9.

13: lr – IntegerInput

On entry: the dimension of the array r.

Constraints:

if $k \geq 6$ , $lr \geq (5 r^{2} + 1) \times k^{2} + (4 r + 3) \times k + 4$ ;
if $k < 6$ , $lr \geq ({(p + q)}^{2} + 1) \times k^{2} + (4 \times (p + q) + 3) \times k + \max (k r (k r + 2), k^{2} {(p + q)}^{2} + l (l + 3) + k^{2} (q + 1)) + 4$ .

Where

r = \max (p, q)

and if

p = 0

l = k (k + 1) / 2

, or if

p \geq 1

l = k (k + 1) / 2 + (p - 1) k^{2}

See Section 9 for some examples of the required size of the array r.

14: state[ $\dim$ ] – IntegerCommunication Array

Note: the dimension,

\dim

, of this array is dictated by the requirements of associated functions that must have been previously called. This array MUST be the same array passed as argument state in the previous call to nag_rand_init_repeatable (g05kfc) or nag_rand_init_nonrepeatable (g05kgc).

On entry: contains information on the selected base generator and its current state.

On exit: contains updated information on the state of the generator.

15: x[ $\dim$ ] – doubleOutput

Note: the dimension, dim, of the array x must be at least

$\max (1, pdx \times n)$ when $order = Nag_ColMajor$ ;
$\max (1, k \times pdx)$ when $order = Nag_RowMajor$ .

Where

X (i, t)

appears in this document, it refers to the array element

$x [(t - 1) \times pdx + i - 1]$ when $order = Nag_ColMajor$ ;
$x [(i - 1) \times pdx + t - 1]$ when $order = Nag_RowMajor$ .

On exit:

X (i, t)

will contain a realization of the

i

th component of

X_{t}

, for

i = 1, 2, \dots, k

and

t = 1, 2, \dots, n

16: pdx – IntegerInput

On entry: the stride separating row or column elements (depending on the value of order) in the array x.

Constraints:

if $order = Nag_ColMajor$ , $pdx \geq k$ ;
if $order = Nag_RowMajor$ , $pdx \geq n$ .

17: fail – NagError *Input/Output

The NAG error argument (see Section 3.6 in the Essential Introduction).

6 Error Indicators and Warnings

NE_ALLOC_FAIL: Dynamic memory allocation failed.
NE_BAD_PARAM: On entry, argument $⟨value⟩$ had an illegal value.
NE_CLOSE_TO_STATIONARITY: The reference vector cannot be computed because the AR parameters are too close to the boundary of the stationarity region.
NE_INT: On entry, $k = ⟨value⟩$ .
Constraint: $k \geq 1$ .
On entry, lr is not large enough, $lr = ⟨value⟩$ : minimum length required $= ⟨value⟩$ .
On entry, $n = ⟨value⟩$ .
Constraint: $n \geq 0$ .
On entry, $p = ⟨value⟩$ .
Constraint: $p \geq 0$ .
On entry, $pdv = ⟨value⟩$ .
Constraint: $pdv > 0$ .
On entry, $pdx = ⟨value⟩$ .
Constraint: $pdx > 0$ .
On entry, $q = ⟨value⟩$ .
Constraint: $q \geq 0$ .
NE_INT_2: On entry, $pdv = ⟨value⟩$ and $k = ⟨value⟩$ .
Constraint: $pdv \geq k$ .
On entry, $pdx = ⟨value⟩$ and $k = ⟨value⟩$ .
Constraint: $pdx \geq k$ .
On entry, $pdx = ⟨value⟩$ and $n = ⟨value⟩$ .
Constraint: $pdx \geq n$ .
NE_INTERNAL_ERROR: 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.
NE_INVALID_STATE: On entry, state vector has been corrupted or not initialized.
NE_INVERTIBILITY: On entry, the moving average parameter matrices are such that the model is non-invertible.
NE_POS_DEF: On entry, the covariance matrix var is not positive semidefinite to machine precision.
NE_PREV_CALL: k is not the same as when r was set up in a previous call.
Previous value of $k = ⟨value⟩$ and $k = ⟨value⟩$ .
NE_STATIONARY_AR: On entry, the AR parameters are outside the stationarity region.
NE_TOO_MANY_ITER: An excessive number of iterations were required by the NAG function used to evaluate the eigenvalues of the covariance matrix.
An excessive number of iterations were required by the NAG function used to evaluate the eigenvalues of the matrices used to test for stationarity or invertibility.
An excessive number of iterations were required by the NAG function used to evaluate the eigenvalues stored in the reference vector.

7 Accuracy

The accuracy is limited by the matrix computations performed, and this is dependent on the condition of the argument and covariance matrices.

8 Parallelism and Performance

nag_rand_varma (g05pjc) is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.

nag_rand_varma (g05pjc) 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 Users' Note for your implementation for any additional implementation-specific information.

9 Further Comments

Note that, in reference to

fail . code =

NE_INVERTIBILITY, nag_rand_varma (g05pjc) will permit moving average parameters on the boundary of the invertibility region.

The elements of r contain amongst other information details of the spectral decompositions which are used to generate future multivariate Normals. Note that these eigenvectors may not be unique on different machines. For example the eigenvectors corresponding to multiple eigenvalues may be permuted. Although an effort is made to ensure that the eigenvectors have the same sign on all machines, differences in the signs may theoretically still occur.

The following table gives some examples of the required size of the array r, specified by the argument lr, for

k = 1, 2

3

, and for various values of

p

and

q

				$q$

		0	1	2	3

		13	20	31	46
	0	36	56	92	144
		85	124	199	310

		19	30	45	64
	1	52	88	140	208
		115	190	301	448
p
		35	50	69	92
	2	136	188	256	340
		397	508	655	838

		57	76	99	126
	3	268	336	420	520
		877	1024	1207	1426

Note that nag_tsa_arma_roots (g13dxc) may be used to check whether a VARMA model is stationary and invertible.

The time taken depends on the values of

p

q

and especially

n

and

k

10 Example

This program generates two realizations, each of length

48

, from the bivariate AR(1) model

X_{t} - μ = ϕ_{1} (X_{t - 1} - μ) + ε_{t}

with

ϕ_{1} = [\begin{matrix} 0.80 & 0.07 \\ 0.00 & 0.58 \end{matrix}],

μ = [\begin{matrix} 5.00 \\ 9.00 \end{matrix}],

and

Σ = [\begin{matrix} 2.97 & 0 \\ 0.64 & 5.38 \end{matrix}] .

The pseudorandom number generator is initialized by a call to nag_rand_init_repeatable (g05kfc). Then, in the first call to nag_rand_varma (g05pjc),

mode = Nag_InitializeAndGenerate

in order to set up the reference vector before generating the first realization. In the subsequent call

mode = Nag_ReGenerateFromReference

and a new recent history is generated and used to generate the second realization.

NAG Library Function Documentnag_rand_varma (g05pjc)

+− 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

10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results

NAG Library Function Document

nag_rand_varma (g05pjc)