Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

Chapter Contents
Chapter Introduction
NAG Toolbox

NAG Toolbox: nag_rand_times_mv_varma (g05pj)

Purpose

nag_rand_times_mv_varma (g05pj) 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 nag_rand_times_mv_varma (g05pj).

Syntax

[r, state, x, ifail] = g05pj(mode, n, xmean, ip, phi, iq, theta, var, r, state, 'k', k, 'lr', lr)
[r, state, x, ifail] = nag_rand_times_mv_varma(mode, n, xmean, ip, phi, iq, theta, var, r, state, 'k', k, 'lr', lr)

Description

Let the vector Xt = (x1t,x2t,,xkt)T ${X}_{t}={\left({x}_{1t},{x}_{2t},\dots ,{x}_{kt}\right)}^{\mathrm{T}}$, denote a k$k$-dimensional time series which is assumed to follow a vector autoregressive moving average (VARMA) model of the form:
 Xt − μ = φ1(Xt − 1 − μ) + φ2(Xt − 2 − μ) + ⋯ + φp(Xt − p − μ) + εt − θ1εt − 1 − θ2εt − 2 − ⋯ − θqεt − q
$Xt-μ= ϕ1(Xt-1-μ)+ϕ2(Xt-2-μ)+⋯+ϕp(Xt-p-μ)+ εt-θ1εt-1-θ2εt-2-⋯-θqεt-q$
(1)
where εt = (ε1t,ε2t,,εkt)T ${\epsilon }_{t}={\left({\epsilon }_{1t},{\epsilon }_{2t},\dots ,{\epsilon }_{kt}\right)}^{\mathrm{T}}$, is a vector of k$k$ residual series assumed to be Normally distributed with zero mean and covariance matrix Σ$\Sigma$. The components of εt${\epsilon }_{t}$ are assumed to be uncorrelated at non-simultaneous lags. The φi${\varphi }_{i}$'s and θj${\theta }_{j}$'s are k$k$ by k$k$ matrices of parameters. {φi}$\left\{{\varphi }_{i}\right\}$, for i = 1,2,,p$\mathit{i}=1,2,\dots ,p$, are called the autoregressive (AR) parameter matrices, and {θj}$\left\{{\theta }_{j}\right\}$, for j = 1,2,,q$\mathit{j}=1,2,\dots ,q$, the moving average (MA) parameter matrices. The parameters in the model are thus the p$p$ k$k$ by k$k$ φ$\varphi$-matrices, the q$q$ k$k$ by k$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
$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$I$ denotes the k$k$ by 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(φ)$A\left(\varphi \right)$ lie inside the unit circle and invertible if the eigenvalues of B(θ)$B\left(\theta \right)$ lie inside the unit circle.
For k6$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 function computes a realisation of the pre-observed vectors X0,X1,,X1p${X}_{0},{X}_{-1},\dots ,{X}_{1-p}$, ε0,ε1,,ε1q${\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$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 nag_rand_times_mv_varma (g05pj). See the description of the parameter mode in Section [Parameters] for details.
The function returns a realisation of X1,X2,,Xn${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_times_mv_varma (g05pj) may be called again to generate a realisation of Xn + 1,Xn + 2,${X}_{n+1},{X}_{n+2},\dots$, etc. See the description of the parameter mode in Section [Parameters] for details.
Further computational details are given in Shea (1988). Note, however, that nag_rand_times_mv_varma (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 functions nag_rand_init_repeat (g05kf) (for a repeatable sequence if computed sequentially) or nag_rand_init_nonrepeat (g05kg) (for a non-repeatable sequence) must be called prior to the first call to nag_rand_times_mv_varma (g05pj).

References

Barone P (1987) A method for generating independent realisations of a multivariate normal stationary and invertible ARMA(p,q)$\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

Parameters

Compulsory Input Parameters

1:     mode – int64int32nag_int scalar
A code for selecting the operation to be performed by the function.
mode = 0${\mathbf{mode}}=0$
Set up reference vector and compute a realisation of the recent history.
mode = 1${\mathbf{mode}}=1$
Generate terms in the time series using reference vector set up in a prior call to nag_rand_times_mv_varma (g05pj).
mode = 2${\mathbf{mode}}=2$
Combine the operations of mode = 0${\mathbf{mode}}=0$ and 1$1$.
mode = 3${\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 mode = 1${\mathbf{mode}}=1$ or 3$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 nag_rand_times_mv_varma (g05pj).
Constraint: mode = 0${\mathbf{mode}}=0$, 1$1$, 2$2$ or 3$3$.
2:     n – int64int32nag_int scalar
n$n$, the number of observations to be generated.
Constraint: n0${\mathbf{n}}\ge 0$.
3:     xmean(k) – double array
k, the dimension of the array, must satisfy the constraint k1${\mathbf{k}}\ge 1$.
μ$\mu$, the vector of means of the multivariate time series.
4:     ip – int64int32nag_int scalar
p$p$, the number of autoregressive parameter matrices.
Constraint: ip0${\mathbf{ip}}\ge 0$.
5:     phi(k × k × ip${\mathbf{k}}×{\mathbf{k}}×{\mathbf{ip}}$) – double array
Must contain the elements of the ip × k × k${\mathbf{ip}}×{\mathbf{k}}×{\mathbf{k}}$ autoregressive parameter matrices of the model, φ1,φ2,,φp${\varphi }_{1},{\varphi }_{2},\dots ,{\varphi }_{p}$. If phi is considered as a three-dimensional array, dimensioned as phi(k,k,ip)${\mathbf{phi}}\left({\mathbf{k}},{\mathbf{k}},{\mathbf{ip}}\right)$, then the (i,j)$\left(i,j\right)$th element of φl${\varphi }_{\mathit{l}}$ would be stored in phi(i,j,l)${\mathbf{phi}}\left(i,j,\mathit{l}\right)$; that is, phi((l1) × k × k + (j1) × k + i)${\mathbf{phi}}\left(\left(\mathit{l}-1\right)×k×k+\left(j-1\right)×k+i\right)$ must be set equal to the (i,j)$\left(i,j\right)$th element of φl${\varphi }_{\mathit{l}}$, for l = 1,2,,p$\mathit{l}=1,2,\dots ,p$, i = 1,2,,k$i=1,2,\dots ,k$ and j = 1,2,,k$j=1,2,\dots ,k$.
Constraint: the elements of phi must satisfy the stationarity condition.
6:     iq – int64int32nag_int scalar
q$q$, the number of moving average parameter matrices.
Constraint: iq0${\mathbf{iq}}\ge 0$.
7:     theta(k × k × iq${\mathbf{k}}×{\mathbf{k}}×{\mathbf{iq}}$) – double array
Must contain the elements of the iq × k × k${\mathbf{iq}}×{\mathbf{k}}×{\mathbf{k}}$ moving average parameter matrices of the model, θ1,θ2,,θq${\theta }_{1},{\theta }_{2},\dots ,{\theta }_{q}$. If theta is considered as a three-dimensional array, dimensioned as theta(k,k,iq), then the (i,j)$\left(i,j\right)$th element of θl${\theta }_{\mathit{l}}$ would be stored in theta(i,j,l)${\mathbf{theta}}\left(\mathit{i},\mathit{j},\mathit{l}\right)$; that is, theta((l1) × k × k + (j1) × k + i)${\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 (i,j)$\left(\mathit{i},\mathit{j}\right)$th element of θl${\theta }_{\mathit{l}}$, for l = 1,2,,q$\mathit{l}=1,2,\dots ,q$, i = 1,2,,k$\mathit{i}=1,2,\dots ,k$ and j = 1,2,,k$\mathit{j}=1,2,\dots ,k$.
Constraint: the elements of theta must be within the invertibility region.
8:     var(ldvar,k) – double array
ldvar, the first dimension of the array, must satisfy the constraint ldvark$\mathit{ldvar}\ge {\mathbf{k}}$.
var(i,j)${\mathbf{var}}\left(\mathit{i},\mathit{j}\right)$ must contain the (i,j$\mathit{i},\mathit{j}$)th element of Σ$\Sigma$, for i = 1,2,,k$\mathit{i}=1,2,\dots ,{\mathbf{k}}$ and j = 1,2,,k$\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.
9:     r(lr) – double array
lr, the dimension of the array, must satisfy the constraint
• if k6${\mathbf{k}}\ge 6$, lr(5r2 + 1) × k2 + (4r + 3) × k + 4${\mathbf{lr}}\ge \left(5{\mathit{r}}^{2}+1\right)×{{\mathbf{k}}}^{2}+\left(4\mathit{r}+3\right)×{\mathbf{k}}+4$;
• if k < 6${\mathbf{k}}<6$, lr((ip + iq)2 + 1) × k2 +   (4 × (ip + iq) + 3) × k + max {kr(kr + 2),k2(ip + iq)2 + l(l + 3) + k2(iq + 1)} + 4${\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$.
If mode = 1${\mathbf{mode}}=1$ or 3$3$, the array r as output from the previous call to nag_rand_times_mv_varma (g05pj) must be input without any change.
If mode = 0${\mathbf{mode}}=0$ or 2$2$, the contents of r need not be set.
10:   state( : $:$) – int64int32nag_int array
Note: the actual argument supplied must be the array state supplied to the initialization routines nag_rand_init_repeat (g05kf) or nag_rand_init_nonrepeat (g05kg).
Contains information on the selected base generator and its current state.

Optional Input Parameters

1:     k – int64int32nag_int scalar
Default: The dimension of the array xmean and the first dimension of the array var and the second dimension of the array var. (An error is raised if these dimensions are not equal.)
k$k$, the dimension of the multivariate time series.
Constraint: k1${\mathbf{k}}\ge 1$.
2:     lr – int64int32nag_int scalar
Default: The dimension of the array r.
The dimension of the array r as declared in the (sub)program from which nag_rand_times_mv_varma (g05pj) is called.
Constraints:
• if k6${\mathbf{k}}\ge 6$, lr(5r2 + 1) × k2 + (4r + 3) × k + 4${\mathbf{lr}}\ge \left(5{\mathit{r}}^{2}+1\right)×{{\mathbf{k}}}^{2}+\left(4\mathit{r}+3\right)×{\mathbf{k}}+4$;
• if k < 6${\mathbf{k}}<6$, lr((ip + iq)2 + 1) × k2 +   (4 × (ip + iq) + 3) × k + max {kr(kr + 2),k2(ip + iq)2 + l(l + 3) + k2(iq + 1)} + 4${\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 r = max (ip,iq)$\mathit{r}=\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{ip}},{\mathbf{iq}}\right)$ and if ip = 0${\mathbf{ip}}=0$, l = k(k + 1) / 2$\mathit{l}={\mathbf{k}}\left({\mathbf{k}}+1\right)/2$, or if ip1${\mathbf{ip}}\ge 1$, l = k(k + 1) / 2 + (ip1)k2$\mathit{l}={\mathbf{k}}\left({\mathbf{k}}+1\right)/2+\left({\mathbf{ip}}-1\right){{\mathbf{k}}}^{2}$.
See Section [Further Comments] for some examples of the required size of the array r.

ldvar ldx

Output Parameters

1:     r(lr) – double array
Information required for any subsequent calls to the function with mode = 1${\mathbf{mode}}=1$ or 3$3$. See Section [Further Comments].
2:     state( : $:$) – int64int32nag_int array
Note: the actual argument supplied must be the array state supplied to the initialization routines nag_rand_init_repeat (g05kf) or nag_rand_init_nonrepeat (g05kg).
Contains updated information on the state of the generator.
3:     x(ldx,n) – double array
ldxk$\mathit{ldx}\ge {\mathbf{k}}$.
x(i,t)${\mathbf{x}}\left(\mathit{i},\mathit{t}\right)$ will contain a realisation of the i$\mathit{i}$th component of Xt${X}_{\mathit{t}}$, for i = 1,2,,k$\mathit{i}=1,2,\dots ,k$ and t = 1,2,,n$\mathit{t}=1,2,\dots ,n$.
4:     ifail – int64int32nag_int scalar
${\mathrm{ifail}}={\mathbf{0}}$ unless the function detects an error (see [Error Indicators and Warnings]).

Error Indicators and Warnings

Errors or warnings detected by the function:

Cases prefixed with W are classified as warnings and do not generate an error of type NAG:error_n. See nag_issue_warnings.

ifail = 1${\mathbf{ifail}}=1$
On entry, mode0${\mathbf{mode}}\ne 0$, 1$1$, 2$2$ or 3$3$.
ifail = 2${\mathbf{ifail}}=2$
On entry, n < 0${\mathbf{n}}<0$.
ifail = 3${\mathbf{ifail}}=3$
On entry, k < 1${\mathbf{k}}<1$.
ifail = 5${\mathbf{ifail}}=5$
On entry, ip < 0${\mathbf{ip}}<0$.
ifail = 6${\mathbf{ifail}}=6$
The autoregressive parameter matrices, stored in phi, are such that the model is non-stationary.
ifail = 7${\mathbf{ifail}}=7$
On entry, iq < 0${\mathbf{iq}}<0$.
ifail = 8${\mathbf{ifail}}=8$
On entry, the moving average parameter matrices, stored in theta, are such that the model is non-invertible.
ifail = 9${\mathbf{ifail}}=9$
The covariance matrix Σ$\Sigma$, stored in var, is not positive semidefinite.
ifail = 10${\mathbf{ifail}}=10$
On entry, ldvar < k$\mathit{ldvar}<{\mathbf{k}}$.
ifail = 11${\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 nag_rand_times_mv_varma (g05pj) with mode = 0${\mathbf{mode}}=0$ or 2$2$.
ifail = 12${\mathbf{ifail}}=12$
On entry, lr is too small.
ifail = 13${\mathbf{ifail}}=13$
 On entry, state vector was not initialized or has been corrupted.
ifail = 15${\mathbf{ifail}}=15$
On entry, ldx < k$\mathit{ldx}<{\mathbf{k}}$.
ifail = 20${\mathbf{ifail}}=20$
This is an unlikely exit brought about by an excessive number of iterations being needed by the NAG Toolbox function used to evaluate the eigenvalues of A(φ)$A\left(\varphi \right)$ or B(θ)$B\left(\theta \right)$.
ifail = 21${\mathbf{ifail}}=21$
nag_rand_times_mv_varma (g05pj) 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.
ifail = 22${\mathbf{ifail}}=22$
This is an unlikely exit brought about by an excessive number of iterations being needed by the NAG Toolbox function used to evaluate the eigenvalues of the covariance matrix.
W ifail = 23${\mathbf{ifail}}=23$
nag_rand_times_mv_varma (g05pj) 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 Toolbox function used to evaluate eigenvalues to be stored in the array r. If this error occurs please contact NAG.

Accuracy

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

Note that, in reference to ${\mathbf{ifail}}={\mathbf{8}}$, nag_rand_times_mv_varma (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$k=1,2$ or 3$3$, and for various values of p$p$ and q$q$.
 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_uni_arma_roots (g13dx) may be used to check whether a VARMA model is stationary and invertible.
The time taken depends on the values of p$p$, q$q$ and especially n$n$ and k$k$.

Example

```function nag_rand_times_mv_varma_example
% Initialize the seed
seed = [int64(1762543)];
% genid and subid identify the base generator
genid = int64(1);
subid =  int64(1);

n = int64(48);
xmean = [5; 9];
ip = int64(1);
iq = int64(0);
phi = [0.8; 0; 0.07; 0.58];
theta = [];
var = [2.97, 0; 0.64, 5.38];
r = zeros(600, 1);

% Initialize the generator to a repeatable sequence
[state, ifail] = nag_rand_init_repeat(genid, subid, seed);

% Generate the first realisation
% Use mode = 2 to set up R and generate values
mode = int64(2);
[r, state, x, ifail] = nag_rand_times_mv_varma(mode, n, xmean, ip, phi, iq, theta, var, r, state);
% Display the results
if ifail == 0
fprintf('\n Realisation Number 1\n\n');
for i=1:2
fprintf(' Series number %d\n', i);
fprintf(' ---------------\n');
for j=1:6
fprintf(' %8.3f', x(i, (j-1)*8+1:j*8));
fprintf('\n');
end
fprintf('\n');
end
fprintf('\n');
else
fprintf(' nag_rand_times_mv_varma exited with ifail = %d\n', ifail);
end

% Generate a second realisation
% Use mode = 3 to reset the series and generate values
mode = int64(3);
[r, state, x, ifail] = nag_rand_times_mv_varma(mode, n, xmean, ip, phi, iq, theta, var, r, state);
% Display the results
if ifail == 0
fprintf('\n Realisation Number 2\n\n');
for i=1:2
fprintf(' Series number %d\n', i);
fprintf(' ---------------\n');
for j=1:6
fprintf(' %8.3f', x(i, (j-1)*8+1:j*8));
fprintf('\n');
end
fprintf('\n');
end
fprintf('\n');
else
fprintf(' nag_rand_times_mv_varma exited with ifail = %d\n', ifail);
end
```
```

Realisation Number 1

Series number 1
---------------
4.833    2.813    3.224    3.825    1.023    1.415    2.184    3.005
5.547    4.832    4.705    5.484    9.407   10.335    8.495    7.478
6.373    6.692    6.698    6.976    6.200    4.458    2.520    3.517
3.054    5.439    5.699    7.136    5.750    8.497    9.563   11.604
9.020   10.063    7.976    5.927    4.992    4.222    3.982    7.107
3.554    7.045    7.025    4.106    5.106    5.954    8.026    7.212

Series number 2
---------------
8.458    9.140   10.866   10.975    9.245    5.054    5.023   12.486
10.534   10.590   11.376    8.793   14.445   13.237   11.030    8.405
7.187    8.291    5.920    9.390   10.055    6.222    7.751   10.604
12.441   10.664   10.960    8.022   10.073   12.870   12.665   14.064
11.867   12.894   10.546   12.754    8.594    9.042   12.029   12.557
9.746    5.487    5.500    8.629    9.723    8.632    6.383   12.484

Realisation Number 2

Series number 1
---------------
5.396    4.811    2.685    5.824    2.449    3.563    5.663    6.209
3.130    4.308    4.333    4.903    1.770    1.278    1.340   -0.527
1.745    3.211    4.478    5.170    5.365    4.852    6.080    6.464
2.765    2.148    6.641    7.224   10.316    7.102    5.604    3.934
4.839    3.698    5.210    5.384    7.652    7.315    7.332    7.561
7.537    7.788    6.868    7.575    6.108    6.188    8.132   10.310

Series number 2
---------------
11.345   10.070   13.654   12.409   11.329   13.054   12.465    9.867
10.263   13.394   10.553   10.331    7.814    8.747   10.025   11.167
10.626    9.366    9.607    9.662   10.492   10.766   11.512   10.813
10.799    8.780    9.221   14.245   11.575   10.620    8.282    5.447
9.935    9.386   11.627   10.066   11.394    7.951    7.907   12.616
15.246    9.962   13.216   11.350   11.227    6.021    6.968   12.428

```
```function g05pj_example
% Initialize the seed
seed = [int64(1762543)];
% genid and subid identify the base generator
genid = int64(1);
subid =  int64(1);

n = int64(48);
xmean = [5; 9];
ip = int64(1);
iq = int64(0);
phi = [0.8; 0; 0.07; 0.58];
theta = [];
var = [2.97, 0; 0.64, 5.38];
r = zeros(600, 1);

% Initialize the generator to a repeatable sequence
[state, ifail] = g05kf(genid, subid, seed);

% Generate the first realisation
% Use mode = 2 to set up R and generate values
mode = int64(2);
[r, state, x, ifail] = g05pj(mode, n, xmean, ip, phi, iq, theta, var, r, state);
% Display the results
if ifail == 0
fprintf('\n Realisation Number 1\n\n');
for i=1:2
fprintf(' Series number %d\n', i);
fprintf(' ---------------\n');
for j=1:6
fprintf(' %8.3f', x(i, (j-1)*8+1:j*8));
fprintf('\n');
end
fprintf('\n');
end
fprintf('\n');
else
fprintf(' g05pj exited with ifail = %d\n', ifail);
end

% Generate a second realisation
% Use mode = 3 to reset the series and generate values
mode = int64(3);
[r, state, x, ifail] = g05pj(mode, n, xmean, ip, phi, iq, theta, var, r, state);
% Display the results
if ifail == 0
fprintf('\n Realisation Number 2\n\n');
for i=1:2
fprintf(' Series number %d\n', i);
fprintf(' ---------------\n');
for j=1:6
fprintf(' %8.3f', x(i, (j-1)*8+1:j*8));
fprintf('\n');
end
fprintf('\n');
end
fprintf('\n');
else
fprintf(' g05pj exited with ifail = %d\n', ifail);
end
```
```

Realisation Number 1

Series number 1
---------------
4.833    2.813    3.224    3.825    1.023    1.415    2.184    3.005
5.547    4.832    4.705    5.484    9.407   10.335    8.495    7.478
6.373    6.692    6.698    6.976    6.200    4.458    2.520    3.517
3.054    5.439    5.699    7.136    5.750    8.497    9.563   11.604
9.020   10.063    7.976    5.927    4.992    4.222    3.982    7.107
3.554    7.045    7.025    4.106    5.106    5.954    8.026    7.212

Series number 2
---------------
8.458    9.140   10.866   10.975    9.245    5.054    5.023   12.486
10.534   10.590   11.376    8.793   14.445   13.237   11.030    8.405
7.187    8.291    5.920    9.390   10.055    6.222    7.751   10.604
12.441   10.664   10.960    8.022   10.073   12.870   12.665   14.064
11.867   12.894   10.546   12.754    8.594    9.042   12.029   12.557
9.746    5.487    5.500    8.629    9.723    8.632    6.383   12.484

Realisation Number 2

Series number 1
---------------
5.396    4.811    2.685    5.824    2.449    3.563    5.663    6.209
3.130    4.308    4.333    4.903    1.770    1.278    1.340   -0.527
1.745    3.211    4.478    5.170    5.365    4.852    6.080    6.464
2.765    2.148    6.641    7.224   10.316    7.102    5.604    3.934
4.839    3.698    5.210    5.384    7.652    7.315    7.332    7.561
7.537    7.788    6.868    7.575    6.108    6.188    8.132   10.310

Series number 2
---------------
11.345   10.070   13.654   12.409   11.329   13.054   12.465    9.867
10.263   13.394   10.553   10.331    7.814    8.747   10.025   11.167
10.626    9.366    9.607    9.662   10.492   10.766   11.512   10.813
10.799    8.780    9.221   14.245   11.575   10.620    8.282    5.447
9.935    9.386   11.627   10.066   11.394    7.951    7.907   12.616
15.246    9.962   13.216   11.350   11.227    6.021    6.968   12.428

```