# Syntax

C# |
---|

public static void g05pj( int mode, int n, int k, double[] xmean, int ip, double[] phi, int iq, double[] theta, double[,] var, double[] r, G05..::..G05State g05state, double[,] x, out int ifail ) |

Visual Basic |
---|

Public Shared Sub g05pj ( _ mode As Integer, _ n As Integer, _ k As Integer, _ xmean As Double(), _ ip As Integer, _ phi As Double(), _ iq As Integer, _ theta As Double(), _ var As Double(,), _ r As Double(), _ g05state As G05..::..G05State, _ x As Double(,), _ <OutAttribute> ByRef ifail As Integer _ ) |

Visual C++ |
---|

public: static void g05pj( int mode, int n, int k, array<double>^ xmean, int ip, array<double>^ phi, int iq, array<double>^ theta, array<double,2>^ var, array<double>^ r, G05..::..G05State^ g05state, array<double,2>^ x, [OutAttribute] int% ifail ) |

F# |
---|

static member g05pj : mode : int * n : int * k : int * xmean : float[] * ip : int * phi : float[] * iq : int * theta : float[] * var : float[,] * r : float[] * g05state : G05..::..G05State * x : float[,] * ifail : int byref -> unit |

#### Parameters

- mode
- Type: System..::..Int32
*On entry*: a code for selecting the operation to be performed by the method.- ${\mathbf{mode}}=0$
- Set up reference vector and compute a realization of the recent history.
- ${\mathbf{mode}}=1$
- Generate terms in the time series using reference vector set up in a prior call to g05pj.
- ${\mathbf{mode}}=2$
- Combine the operations of ${\mathbf{mode}}=0$ and $1$.
- ${\mathbf{mode}}=3$
- 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.

*Constraint*: ${\mathbf{mode}}=0$, $1$, $2$ or $3$.

- n
- Type: System..::..Int32
*On entry*: $n$, the number of observations to be generated.*Constraint*: ${\mathbf{n}}\ge 0$.

- k
- Type: System..::..Int32
*On entry*: $k$, the dimension of the multivariate time series.*Constraint*: ${\mathbf{k}}\ge 1$.

- xmean
- Type: array<System..::..Double>[]()[][]An array of size [k]
*On entry*: $\mu $, the vector of means of the multivariate time series.

- ip
- Type: System..::..Int32
*On entry*: $p$, the number of autoregressive parameter matrices.*Constraint*: ${\mathbf{ip}}\ge 0$.

- phi
- Type: array<System..::..Double>[]()[][]An array of size [${\mathbf{k}}\times {\mathbf{k}}\times {\mathbf{ip}}$]
*On entry*: must contain the elements of the ${\mathbf{ip}}\times {\mathbf{k}}\times {\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}}[{\mathbf{k}}-1,{\mathbf{k}}-1,{\mathbf{ip}}-1]$, then the $\left(i,j\right)$th element of ${\varphi}_{\mathit{l}}$ would be stored in ${\mathbf{phi}}[i-1,j-1,\mathit{l}-1]$; that is, ${\mathbf{phi}}\left[\left(\mathit{l}-1\right)\times k\times k+\left(j-1\right)\times k+i-1\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.

- iq
- Type: System..::..Int32
*On entry*: $q$, the number of moving average parameter matrices.*Constraint*: ${\mathbf{iq}}\ge 0$.

- theta
- Type: array<System..::..Double>[]()[][]An array of size [${\mathbf{k}}\times {\mathbf{k}}\times {\mathbf{iq}}$]
*On entry*: must contain the elements of the ${\mathbf{iq}}\times {\mathbf{k}}\times {\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}}[\mathit{i}-1,\mathit{j}-1,\mathit{l}-1]$; that is, ${\mathbf{theta}}\left[\left(\mathit{l}-1\right)\times k\times k+\left(\mathit{j}-1\right)\times k+\mathit{i}-1\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.

- var
- Type: array<System..::..Double,2>[,](,)[,][,]An array of size [dim1, k]
**Note:**dim1 must satisfy the constraint: $\mathrm{dim1}\ge {\mathbf{k}}$*On entry*: ${\mathbf{var}}[\mathit{i}-1,\mathit{j}-1]$ 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.

- r
- Type: array<System..::..Double>[]()[][]An array of size [lr]
*On entry*: if ${\mathbf{mode}}=1$ or $3$, the array r as output from the previous call to g05pj must be input without any change.*On exit*: information required for any subsequent calls to the method with ${\mathbf{mode}}=1$ or $3$. See [Further Comments].

- g05state
- Type: NagLibrary..::..G05..::..G05StateAn Object of type G05.G05State.

- x
- Type: array<System..::..Double,2>[,](,)[,][,]An array of size [dim1, n]
**Note:**dim1 must satisfy the constraint: $\mathrm{dim1}\ge {\mathbf{k}}$*On exit*: ${\mathbf{x}}[\mathit{i}-1,\mathit{t}-1]$ will contain a realization of the $\mathit{i}$th component of ${X}_{\mathit{t}}$, for $\mathit{i}=1,2,\dots ,k$ and $\mathit{t}=1,2,\dots ,n$.

- ifail
- Type: System..::..Int32%
*On exit*: ${\mathbf{ifail}}={0}$ unless the method detects an error or a warning has been flagged (see [Error Indicators and Warnings]).

# 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:

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

where $I$ denotes the $k$ by $k$ identity matrix.

$$\begin{array}{cc}{X}_{t}-\mu =& {\varphi}_{1}\left({X}_{t-1}-\mu \right)+{\varphi}_{2}\left({X}_{t-2}-\mu \right)+\cdots +{\varphi}_{p}\left({X}_{t-p}-\mu \right)+\\ & {\epsilon}_{t}-{\theta}_{1}{\epsilon}_{t-1}-{\theta}_{2}{\epsilon}_{t-2}-\cdots -{\theta}_{q}{\epsilon}_{t-q}\end{array}$$ | (1) |

$$A\left(\varphi \right)={\left[\begin{array}{ccccccc}{\varphi}_{1}& I& 0& .& .& .& 0\\ {\varphi}_{2}& 0& I& 0& .& .& 0\\ .& & & .& & & \\ .& & & & .& & \\ .& & & & & .& \\ {\varphi}_{p-1}& 0& .& .& .& 0& I\\ {\varphi}_{p}& 0& .& .& .& 0& 0\end{array}\right]}_{pk\times pk}\text{\hspace{1em} and \hspace{1em}}B\left(\theta \right)={\left[\begin{array}{ccccccc}{\theta}_{1}& I& 0& .& .& .& 0\\ {\theta}_{2}& 0& I& 0& .& .& 0\\ .& & & .& & & \\ .& & & & .& & \\ .& & & & & .& \\ {\theta}_{q-1}& 0& .& .& .& 0& I\\ {\theta}_{q}& 0& .& .& .& 0& 0\end{array}\right]}_{qk\times qk}$$ |

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 realization of the state vector at time zero computed. For all other cases the method computes a realization 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 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 g05pj. See the description of the parameter mode in [Parameters] for details.

The method 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 g05pj may be called again to generate a realization of ${X}_{n+1},{X}_{n+2},\dots $, etc. See the description of the parameter mode in [Parameters] for details.

Further computational details are given in Shea (1988). Note, however, that g05pj 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 methods (G05KFF not in this release) (for a repeatable sequence if computed sequentially) or (G05KGF not in this release) (for a non-repeatable sequence) must be called prior to the first call to g05pj.

# References

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–130Shea 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# Error Indicators and Warnings

Errors or warnings detected by the method:

Some error messages may refer to parameters that are dropped from this interface
(LDVAR, LDX) In these
cases, an error in another parameter has usually caused an incorrect value to be inferred.

- ${\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$

- ${\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 method used to evaluate the eigenvalues of $A\left(\varphi \right)$ or $B\left(\theta \right)$.

- ${\mathbf{ifail}}=21$

- ${\mathbf{ifail}}=22$
- This is an unlikely exit brought about by an excessive number of iterations being needed by the NAG Library method used to evaluate the eigenvalues of the covariance matrix.

- ${\mathbf{ifail}}=23$

- ${\mathbf{ifail}}=-9000$
- An error occured, see message report.
- ${\mathbf{ifail}}=-6000$
- Invalid Parameters $\u2329\mathit{\text{value}}\u232a$
- ${\mathbf{ifail}}=-4000$
- Invalid dimension for array $\u2329\mathit{\text{value}}\u232a$
- ${\mathbf{ifail}}=-8000$
- Negative dimension for array $\u2329\mathit{\text{value}}\u232a$
- ${\mathbf{ifail}}=-6000$
- Invalid Parameters $\u2329\mathit{\text{value}}\u232a$

# Accuracy

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

# Parallelism and Performance

None.

# Further Comments

Note that, in reference to ${\mathbf{ifail}}={8}$, g05pj 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 g13dx 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$.