G05 Chapter Contents
G05 Chapter Introduction
NAG Library Manual

# NAG Library Routine DocumentG05PJF

Note:  before using this routine, please read the Users' Note for your implementation to check the interpretation of bold italicised terms and other implementation-dependent details.

## 1  Purpose

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

## 2  Specification

 SUBROUTINE G05PJF ( MODE, N, K, XMEAN, IP, PHI, IQ, THETA, VAR, LDVAR, R, LR, STATE, X, LDX, IFAIL)
 INTEGER MODE, N, K, IP, IQ, LDVAR, LR, STATE(*), LDX, IFAIL REAL (KIND=nag_wp) XMEAN(K), PHI(K*K*IP), THETA(K*K*IQ), VAR(LDVAR,K), R(LR), X(LDX,N)

## 3  Description

Let the vector ${X}_{t}={\left({x}_{1t},{x}_{2t},\dots ,{x}_{kt}\right)}^{\mathrm{T}}$, denote a $k$-dimensional time series which is assumed to follow a vector autoregressive moving average (VARMA) model of the form:
 $Xt-μ= ϕ1Xt-1-μ+ϕ2Xt-2-μ+⋯+ϕpXt-p-μ+ εt-θ1εt-1-θ2εt-2-⋯-θqεt-q$ (1)
where ${\epsilon }_{t}={\left({\epsilon }_{1t},{\epsilon }_{2t},\dots ,{\epsilon }_{kt}\right)}^{\mathrm{T}}$, is a vector of $k$ residual series assumed to be Normally distributed with zero mean and covariance matrix $\Sigma$. The components of ${\epsilon }_{t}$ are assumed to be uncorrelated at non-simultaneous lags. The ${\varphi }_{i}$'s and ${\theta }_{j}$'s are $k$ by $k$ matrices of parameters. $\left\{{\varphi }_{i}\right\}$, for $\mathit{i}=1,2,\dots ,p$, are called the autoregressive (AR) parameter matrices, and $\left\{{\theta }_{j}\right\}$, for $\mathit{j}=1,2,\dots ,q$, the moving average (MA) parameter matrices. The parameters in the model are thus the $p$ $k$ by $k$ $\varphi$-matrices, the $q$ $k$ by $k$ $\theta$-matrices, the mean vector $\mu$ and the residual error covariance matrix $\Sigma$. 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 . . . 0 I θq 0 . . . 0 0 qk×qk$
where $I$ denotes the $k$ by $k$ identity matrix.
The model (1) must be both stationary and invertible. The model is said to be stationary if the eigenvalues of $A\left(\varphi \right)$ lie inside the unit circle and invertible if the eigenvalues of $B\left(\theta \right)$ lie inside the unit circle.
For $k\ge 6$ the VARMA model (1) is recast into state space form and a realisation of the state vector at time zero computed. For all other cases the routine computes a realisation of the pre-observed vectors ${X}_{0},{X}_{-1},\dots ,{X}_{1-p}$, ${\epsilon }_{0},{\epsilon }_{-1},\dots ,{\epsilon }_{1-q}$, from (1), see Shea (1988). This realisation 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 realisation of the time series may be generated more efficiently using the information in a reference vector created during a previous call to G05PJF. See the description of the parameter MODE in Section 5 for details.
The routine returns a realisation 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 G05PJF may be called again to generate a realisation of ${X}_{n+1},{X}_{n+2},\dots$, etc. See the description of the parameter MODE in Section 5 for details.
Further computational details are given in Shea (1988). Note, however, that G05PJF 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 routines G05KFF (for a repeatable sequence if computed sequentially) or G05KGF (for a non-repeatable sequence) must be called prior to the first call to G05PJF.
Barone P (1987) A method for generating independent realisations of a multivariate normal stationary and invertible ARMA$\left(p,q\right)$ 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  Parameters

1:     MODE – INTEGERInput
On entry: a code for selecting the operation to be performed by the routine.
${\mathbf{MODE}}=0$
Set up reference vector and compute a realisation of the recent history.
${\mathbf{MODE}}=1$
Generate terms in the time series using reference vector set up in a prior call to G05PJF.
${\mathbf{MODE}}=2$
Combine the operations of ${\mathbf{MODE}}=0$ and $1$.
${\mathbf{MODE}}=3$
A new realisation of the recent history is computed using information stored in the reference vector, and the following sequence of time series values are generated.
If ${\mathbf{MODE}}=1$ or $3$, then you must ensure that the reference vector R and the values of K, IP, IQ, XMEAN, PHI, THETA, VAR and LDVAR have not been changed between calls to G05PJF.
Constraint: ${\mathbf{MODE}}=0$, $1$, $2$ or $3$.
2:     N – INTEGERInput
On entry: $n$, the number of observations to be generated.
Constraint: ${\mathbf{N}}\ge 0$.
3:     K – INTEGERInput
On entry: $k$, the dimension of the multivariate time series.
Constraint: ${\mathbf{K}}\ge 1$.
4:     XMEAN(K) – REAL (KIND=nag_wp) arrayInput
On entry: $\mu$, the vector of means of the multivariate time series.
5:     IP – INTEGERInput
On entry: $p$, the number of autoregressive parameter matrices.
Constraint: ${\mathbf{IP}}\ge 0$.
6:     PHI(${\mathbf{K}}×{\mathbf{K}}×{\mathbf{IP}}$) – REAL (KIND=nag_wp) arrayInput
On entry: must contain the elements of the ${\mathbf{IP}}×{\mathbf{K}}×{\mathbf{K}}$ autoregressive parameter matrices of the model, ${\varphi }_{1},{\varphi }_{2},\dots ,{\varphi }_{p}$. If PHI is considered as a three-dimensional array, dimensioned as ${\mathbf{PHI}}\left({\mathbf{K}},{\mathbf{K}},{\mathbf{IP}}\right)$, then the $\left(i,j\right)$th element of ${\varphi }_{\mathit{l}}$ would be stored in ${\mathbf{PHI}}\left(i,j,\mathit{l}\right)$; that is, ${\mathbf{PHI}}\left(\left(\mathit{l}-1\right)×k×k+\left(j-1\right)×k+i\right)$ must be set equal to the $\left(i,j\right)$th element of ${\varphi }_{\mathit{l}}$, for $\mathit{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.
7:     IQ – INTEGERInput
On entry: $q$, the number of moving average parameter matrices.
Constraint: ${\mathbf{IQ}}\ge 0$.
8:     THETA(${\mathbf{K}}×{\mathbf{K}}×{\mathbf{IQ}}$) – REAL (KIND=nag_wp) arrayInput
On entry: must contain the elements of the ${\mathbf{IQ}}×{\mathbf{K}}×{\mathbf{K}}$ moving average parameter matrices of the model, ${\theta }_{1},{\theta }_{2},\dots ,{\theta }_{q}$. If THETA is considered as a three-dimensional array, dimensioned as THETA(K,K,IQ), then the $\left(i,j\right)$th element of ${\theta }_{\mathit{l}}$ would be stored in ${\mathbf{THETA}}\left(\mathit{i},\mathit{j},\mathit{l}\right)$; that is, ${\mathbf{THETA}}\left(\left(\mathit{l}-1\right)×k×k+\left(\mathit{j}-1\right)×k+\mathit{i}\right)$ must be set equal to the $\left(\mathit{i},\mathit{j}\right)$th element of ${\theta }_{\mathit{l}}$, for $\mathit{l}=1,2,\dots ,q$, $\mathit{i}=1,2,\dots ,k$ and $\mathit{j}=1,2,\dots ,k$.
Constraint: the elements of THETA must be within the invertibility region.
9:     VAR(LDVAR,K) – REAL (KIND=nag_wp) arrayInput
On entry: ${\mathbf{VAR}}\left(\mathit{i},\mathit{j}\right)$ must contain the ($\mathit{i},\mathit{j}$)th element of $\Sigma$, for $\mathit{i}=1,2,\dots ,{\mathbf{K}}$ and $\mathit{j}=1,2,\dots ,{\mathbf{K}}$. Only the lower triangle is required.
Constraint: the elements of VAR must be such that $\Sigma$ is positive semidefinite.
10:   LDVAR – INTEGERInput
On entry: the first dimension of the array VAR as declared in the (sub)program from which G05PJF is called.
Constraint: ${\mathbf{LDVAR}}\ge {\mathbf{K}}$.
11:   R(LR) – REAL (KIND=nag_wp) arrayCommunication Array
On entry: if ${\mathbf{MODE}}=1$ or $3$, the array R as output from the previous call to G05PJF must be input without any change.
If ${\mathbf{MODE}}=0$ or $2$, the contents of R need not be set.
On exit: information required for any subsequent calls to the routine with ${\mathbf{MODE}}=1$ or $3$. See Section 8.
12:   LR – INTEGERInput
On entry: the dimension of the array R as declared in the (sub)program from which G05PJF is called.
Constraints:
• if ${\mathbf{K}}\ge 6$, ${\mathbf{LR}}\ge \left(5{\mathit{r}}^{2}+1\right)×{{\mathbf{K}}}^{2}+\left(4\mathit{r}+3\right)×{\mathbf{K}}+4$;
• if ${\mathbf{K}}<6$, ${\mathbf{LR}}\ge \left({\left({\mathbf{IP}}+{\mathbf{IQ}}\right)}^{2}+1\right)×{{\mathbf{K}}}^{2}+\phantom{\rule{0ex}{0ex}}\left(4×\left({\mathbf{IP}}+{\mathbf{IQ}}\right)+3\right)×{\mathbf{K}}+\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left\{{\mathbf{K}}\mathit{r}\left({\mathbf{K}}\mathit{r}+2\right),{{\mathbf{K}}}^{2}{\left({\mathbf{IP}}+{\mathbf{IQ}}\right)}^{2}+\mathit{l}\left(\mathit{l}+3\right)+{{\mathbf{K}}}^{2}\left({\mathbf{IQ}}+1\right)\right\}+4$.
Where $\mathit{r}=\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{IP}},{\mathbf{IQ}}\right)$ and if ${\mathbf{IP}}=0$, $\mathit{l}={\mathbf{K}}\left({\mathbf{K}}+1\right)/2$, or if ${\mathbf{IP}}\ge 1$, $\mathit{l}={\mathbf{K}}\left({\mathbf{K}}+1\right)/2+\left({\mathbf{IP}}-1\right){{\mathbf{K}}}^{2}$.
See Section 8 for some examples of the required size of the array R.
13:   STATE($*$) – INTEGER arrayCommunication Array
Note: the actual argument supplied must be the array STATE supplied to the initialization routines G05KFF or G05KGF.
On entry: contains information on the selected base generator and its current state.
On exit: contains updated information on the state of the generator.
14:   X(LDX,N) – REAL (KIND=nag_wp) arrayOutput
On exit: ${\mathbf{X}}\left(\mathit{i},\mathit{t}\right)$ will contain a realisation of the $\mathit{i}$th component of ${X}_{\mathit{t}}$, for $\mathit{i}=1,2,\dots ,k$ and $\mathit{t}=1,2,\dots ,n$.
15:   LDX – INTEGERInput
On entry: the first dimension of the array X as declared in the (sub)program from which G05PJF is called.
Constraint: ${\mathbf{LDX}}\ge {\mathbf{K}}$.
16:   IFAIL – INTEGERInput/Output
On entry: IFAIL must be set to $0$, $-1\text{​ or ​}1$. If you are unfamiliar with this parameter you should refer to Section 3.3 in the Essential Introduction for details.
For environments where it might be inappropriate to halt program execution when an error is detected, the value $-1\text{​ or ​}1$ is recommended. If the output of error messages is undesirable, then the value $1$ is recommended. Otherwise, if you are not familiar with this parameter, the recommended value is $0$. When the value $-\mathbf{1}\text{​ or ​}\mathbf{1}$ is used it is essential to test the value of IFAIL on exit.
On exit: ${\mathbf{IFAIL}}={\mathbf{0}}$ unless the routine detects an error or a warning has been flagged (see Section 6).

## 6  Error Indicators and Warnings

If on entry ${\mathbf{IFAIL}}={\mathbf{0}}$ or $-{\mathbf{1}}$, explanatory error messages are output on the current error message unit (as defined by X04AAF).
Errors or warnings detected by the routine:
${\mathbf{IFAIL}}=1$
On entry, ${\mathbf{MODE}}\ne 0$, $1$, $2$ or $3$.
${\mathbf{IFAIL}}=2$
On entry, ${\mathbf{N}}<0$.
${\mathbf{IFAIL}}=3$
On entry, ${\mathbf{K}}<1$.
${\mathbf{IFAIL}}=5$
On entry, ${\mathbf{IP}}<0$.
${\mathbf{IFAIL}}=6$
The autoregressive parameter matrices, stored in PHI, are such that the model is non-stationary.
${\mathbf{IFAIL}}=7$
On entry, ${\mathbf{IQ}}<0$.
${\mathbf{IFAIL}}=8$
On entry, the moving average parameter matrices, stored in THETA, are such that the model is non-invertible.
${\mathbf{IFAIL}}=9$
The covariance matrix $\Sigma$, stored in VAR, is not positive semidefinite.
${\mathbf{IFAIL}}=10$
On entry, ${\mathbf{LDVAR}}<{\mathbf{K}}$.
${\mathbf{IFAIL}}=11$
Either R has been corrupted or the value of K is not the same as when R was set up in a previous call to G05PJF with ${\mathbf{MODE}}=0$ or $2$.
${\mathbf{IFAIL}}=12$
On entry, LR is too small.
${\mathbf{IFAIL}}=13$
 On entry, STATE vector was not initialized or has been corrupted.
${\mathbf{IFAIL}}=15$
On entry, ${\mathbf{LDX}}<{\mathbf{K}}$.
${\mathbf{IFAIL}}=20$
This is an unlikely exit brought about by an excessive number of iterations being needed by the NAG Library routine used to evaluate the eigenvalues of $A\left(\varphi \right)$ or $B\left(\theta \right)$.
${\mathbf{IFAIL}}=21$
G05PJF has not been able to calculate all the required elements of the array R. This is likely to be because the AR parameters are very close to the boundary of the stationarity region.
${\mathbf{IFAIL}}=22$
This is an unlikely exit brought about by an excessive number of iterations being needed by the NAG Library routine used to evaluate the eigenvalues of the covariance matrix.
${\mathbf{IFAIL}}=23$
G05PJF has not been able to calculate all the required elements of the array R. This is an unlikely exit brought about by an excessive number of iterations being needed by the NAG Library routine used to evaluate eigenvalues to be stored in the array R. If this error occurs please contact NAG.

## 7  Accuracy

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

## 8  Further Comments

Note that, in reference to ${\mathbf{IFAIL}}={\mathbf{8}}$, G05PJF 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 parameter LR, for $k=1,2$ or $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 G13DXF 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$.

## 9  Example

This program generates two realisations, each of length $48$, from the bivariate AR(1) model
 $Xt-μ=ϕ1Xt-1-μ+εt$
with
 $ϕ1= 0.80 0.07 0.00 0.58 ,$
 $μ= 5.00 9.00 ,$
and
 $Σ= 2.97 0 0.64 5.38 .$
The pseudorandom number generator is initialized by a call to G05KFF. Then, in the first call to G05PJF, ${\mathbf{MODE}}=2$ in order to set up the reference vector before generating the first realisation. In the subsequent call ${\mathbf{MODE}}=3$ and a new recent history is generated and used to generate the second realisation.

### 9.1  Program Text

Program Text (g05pjfe.f90)

### 9.2  Program Data

Program Data (g05pjfe.d)

### 9.3  Program Results

Program Results (g05pjfe.r)