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 g05pjc.
The function may be called by the names: g05pjc, nag_rand_times_mv_varma or nag_rand_varma.
3Description
Let the vector , denote a -dimensional time series which is assumed to follow a vector autoregressive moving average (VARMA) model of the form:
(1)
where , is a vector of residual series assumed to be Normally distributed with zero mean and covariance matrix . The components of are assumed to be uncorrelated at non-simultaneous lags. The 's and 's are matrices of parameters. , for , are called the autoregressive (AR) parameter matrices, and , for , the moving average (MA) parameter matrices. The parameters in the model are thus the
-matrices, the
-matrices, the mean vector and the residual error covariance matrix . Let
where denotes the identity matrix.
The model (1) must be both stationary and invertible. The model is said to be stationary if the eigenvalues of lie inside the unit circle and invertible if the eigenvalues of lie inside the unit circle.
For 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 , , 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 .
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 g05pjc. See the description of the argument mode in Section 5 for details.
The function returns a realization of . On a successful exit, the recent history is updated and saved in the array r so that g05pjc may be called again to generate a realization of , etc. See the description of the argument mode in Section 5 for details.
Further computational details are given in Shea (1988). Note, however, that 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 g05kfc (for a repeatable sequence if computed sequentially) or g05kgc (for a non-repeatable sequence) must be called prior to the first call to g05pjc.
4References
Barone P (1987) A method for generating independent realisations of a multivariate normal stationary and invertible ARMA 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
5Arguments
1: – 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 . See Section 3.1.3 in the Introduction to the NAG Library CL Interface for a more detailed explanation of the use of this argument.
Constraint:
or .
2: – Nag_ModeRNGInput
On entry: a code for selecting the operation to be performed by the function.
Set up reference vector and compute a realization of the recent history.
Generate terms in the time series using reference vector set up in a prior call to g05pjc.
Combine the operations of and .
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.
If or , 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 g05pjc.
Constraint:
, , or .
3: – IntegerInput
On entry: , the number of observations to be generated.
Constraint:
.
4: – IntegerInput
On entry: , the dimension of the multivariate time series.
Constraint:
.
5: – const doubleInput
On entry: , the vector of means of the multivariate time series.
6: – IntegerInput
On entry: , the number of autoregressive parameter matrices.
Constraint:
.
7: – const doubleInput
On entry: must contain the elements of the autoregressive parameter matrices of the model, . The th element of is stored in , for , and .
Constraint:
the elements of phi must satisfy the stationarity condition.
8: – IntegerInput
On entry: , the number of moving average parameter matrices.
Constraint:
.
9: – const doubleInput
On entry: must contain the elements of the moving average parameter matrices of the model, . The th element of is stored in
, for , and .
Constraint:
the elements of theta must be within the invertibility region.
10: – const doubleInput
Note: the dimension, dim, of the array
var
must be at least
.
where appears in this document, it refers to the array element
when ;
when .
On entry: must contain the ()th element of , for and . Only the lower triangle is required.
Constraint:
the elements of var must be such that is positive semidefinite.
11: – IntegerInput
On entry: the stride separating row or column elements (depending on the value of order) in the array var.
Constraint:
.
12: – doubleCommunication Array
On entry: if or , the array r as output from the previous call to g05pjc must be input without any change.
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: – doubleOutput
Note: the dimension, dim, of the array
x
must be at least
when ;
when .
where appears in this document, it refers to the array element
when ;
when .
On exit: will contain a realization of the th component of , for and .
16: – IntegerInput
On entry: the stride separating row or column elements (depending on the value of order) in the array x.
Constraints:
if ,
;
if , .
17: – NagError *Input/Output
The NAG error argument (see Section 7 in the Introduction to the NAG Library CL Interface).
6Error Indicators and Warnings
NE_ALLOC_FAIL
Dynamic memory allocation failed.
See Section 3.1.2 in the Introduction to the NAG Library CL Interface for further information.
NE_BAD_PARAM
On entry, argument 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, .
Constraint: .
On entry, lr is not large enough, : minimum length required .
On entry, .
Constraint: .
On entry, .
Constraint: .
On entry, . Constraint: .
On entry, . Constraint: .
On entry, .
Constraint: .
NE_INT_2
On entry, and .
Constraint: .
On entry, and .
Constraint: .
On entry, and .
Constraint: .
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.
See Section 7.5 in the Introduction to the NAG Library CL Interface for further information.
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_NO_LICENCE
Your licence key may have expired or may not have been installed correctly.
See Section 8 in the Introduction to the NAG Library CL Interface for further information.
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 and .
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.
7Accuracy
The accuracy is limited by the matrix computations performed, and this is dependent on the condition of the parameter and covariance matrices.
8Parallelism and Performance
Background information to multithreading can be found in the Multithreading documentation.
g05pjc is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
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 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.
9Further Comments
Note that, in reference to NE_INVERTIBILITY, 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 or , and for various values of and .
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
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 g13dxc may be used to check whether a VARMA model is stationary and invertible.
The time taken depends on the values of , and especially and .
10Example
This program generates two realizations, each of length , from the bivariate AR(1) model
with
and
The pseudorandom number generator is initialized by a call to g05kfc. Then, in the first call to g05pjc, in order to set up the reference vector before generating the first realization. In the subsequent call and a new recent history is generated and used to generate the second realization.