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_tsa_multi_varma_forecast (g13dj)

## Purpose

nag_tsa_multi_varma_forecast (g13dj) computes forecasts of a multivariate time series. It is assumed that a vector ARMA model has already been fitted to the appropriately differenced/transformed time series using nag_tsa_multi_varma_estimate (g13dd). The standard deviations of the forecast errors are also returned. A reference vector is set up so that, should future series values become available, the forecasts and their standard errors may be updated by calling nag_tsa_multi_varma_update (g13dk).

## Syntax

[qq, predz, sefz, ref, ifail] = g13dj(z, tr, id, delta, ip, iq, mean, par, qq, v, lmax, lref, 'k', k, 'n', n, 'lpar', lpar)
[qq, predz, sefz, ref, ifail] = nag_tsa_multi_varma_forecast(z, tr, id, delta, ip, iq, mean, par, qq, v, lmax, lref, 'k', k, 'n', n, 'lpar', lpar)

## Description

Let the vector Zt = (z1t,z2t,,zkt)T ${Z}_{\mathit{t}}={\left({z}_{1\mathit{t}},{z}_{2\mathit{t}},\dots ,{z}_{k\mathit{t}}\right)}^{\mathrm{T}}$, for t = 1,2,,n$\mathit{t}=1,2,\dots ,n$, denote a k$k$-dimensional time series for which forecasts of Zn + 1,Zn + 2,,Zn + lmax${Z}_{n+1},{Z}_{n+2},\dots ,{Z}_{n+{l}_{\mathrm{max}}}$ are required. Let Wt = (w1t,w2t,,wkt)T ${W}_{t}={\left({w}_{1t},{w}_{2t},\dots ,{w}_{kt}\right)}^{\mathrm{T}}$ be defined as follows:
 wit = δi(B)zit * ,  i = 1,2, … ,k, $wit=δi(B)zit*, i=1,2,…,k,$
where δi(B)${\delta }_{i}\left(B\right)$ is the differencing operator applied to the i$i$th series and where zit * ${z}_{it}^{*}$ is equal to either zit${z}_{it}$, sqrt(zit)$\sqrt{{z}_{it}}$ or loge(zit)${\mathrm{log}}_{\mathrm{e}}\left({z}_{it}\right)$ depending on whether or not a transformation was required to stabilize the variance before fitting the model.
If the order of differencing required for the i$i$th series is di${\mathit{d}}_{i}$, then the differencing operator for the i$i$th series is defined by δi(B) = 1δi1Bδi2B2δidiBdi${\delta }_{i}\left(B\right)=1-{\delta }_{i1}B-{\delta }_{i2}{B}^{2}-\cdots -{\delta }_{i{\mathit{d}}_{i}}{B}^{{\mathit{d}}_{i}}$ where B$B$ is the backward shift operator; that is, BZt = Zt1$B{Z}_{t}={Z}_{t-1}$. The differencing parameters δij${\delta }_{\mathit{i}\mathit{j}}$, for i = 1,2,,k$\mathit{i}=1,2,\dots ,k$ and j = 1,2,,di$\mathit{j}=1,2,\dots ,{\mathit{d}}_{\mathit{i}}$, must be supplied by you. If the i$i$th series does not require differencing, then di = 0${\mathit{d}}_{i}=0$.
Wt${W}_{t}$ is assumed to follow a multivariate ARMA model of the form:
 Wt − μ = φ1(Wt − 1 − μ) + φ2(Wt − 2 − μ) + ⋯ + φp(Wt − p − μ) + εt − θ1εt − 1 − ⋯ − θqεt − q, $Wt-μ=ϕ1(Wt-1-μ)+ϕ2(Wt-2-μ)+⋯+ϕp(Wt-p-μ)+εt-θ1εt-1-⋯-θqεt-q,$ (1)
where εt = (ε1t,ε2t,,εkt)T ${\epsilon }_{\mathit{t}}={\left({\epsilon }_{1\mathit{t}},{\epsilon }_{2\mathit{t}},\dots ,{\epsilon }_{k\mathit{t}}\right)}^{\mathrm{T}}$, for t = 1,2,,n$\mathit{t}=1,2,\dots ,n$, is a vector of k$k$ residual series assumed to be Normally distributed with zero mean and positive definite covariance matrix Σ$\Sigma$. The components of εt${\epsilon }_{t}$ are assumed to be uncorrelated at non-simultaneous lags. The φi${\varphi }_{i}$ and θj${\theta }_{j}$ are k$k$ by k$k$ matrices of parameters. The matrices φi${\varphi }_{\mathit{i}}$, for i = 1,2,,p$\mathit{i}=1,2,\dots ,p$, are the autoregressive (AR) parameter matrices, and the matrices θi${\theta }_{\mathit{i}}$, for i = 1,2,,q$\mathit{i}=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$. The ARMA model (1) must be both stationary and invertible; see nag_tsa_uni_arma_roots (g13dx) for a method of checking these conditions.
The ARMA model (1) may be rewritten as
 φ(B)(δ(B)Zt * − μ) = θ(B)εt, $ϕ(B)(δ(B)Zt*-μ)=θ(B)εt,$
where φ(B)$\varphi \left(B\right)$ and θ(B)$\theta \left(B\right)$ are the autoregressive and moving average polynomials and δ(B)$\delta \left(B\right)$ denotes the k$k$ by k$k$ diagonal matrix whose i$i$th diagonal elements is δi(B)${\delta }_{i}\left(B\right)$ and Zt * = (z1t * , z2t * zkt * )T ${Z}_{t}^{*}={\left({z}_{1t}^{*},{z}_{2t}^{*}\dots {z}_{kt}^{*}\right)}^{\mathrm{T}}$.
This may be rewritten as
 φ(B)δ(B)Zt * = φ(B)μ + θ(B)εt $ϕ(B)δ(B)Zt*=ϕ(B)μ+θ(B)εt$
or
 Zt * = τ + ψ (B)εt = τ + εt + ψ1εt − 1 + ψ2εt − 2 + ⋯ $Zt*=τ+ψ (B)εt=τ+εt+ψ1εt- 1+ψ2εt- 2+⋯$
where ψ(B) = δ1(B)φ1(B)θ(B)$\psi \left(B\right)={\delta }^{-1}\left(B\right){\varphi }^{-1}\left(B\right)\theta \left(B\right)$ and τ = δ1(B)μ$\tau ={\delta }^{-1}\left(B\right)\mu$ is a vector of length k$k$.
Forecasts are computed using a multivariate version of the procedure described in Box and Jenkins (1976). If n * (l)${\stackrel{^}{Z}}_{n}^{*}\left(l\right)$ denotes the forecast of Zn + l * ${Z}_{n+l}^{*}$, then n * (l)${\stackrel{^}{Z}}_{n}^{*}\left(l\right)$ is taken to be that linear function of Zn * ,Zn1 * ,${Z}_{n}^{*},{Z}_{n-1}^{*},\dots \text{}$ which minimizes the elements of E{en(l)en(l)}$E\left\{{e}_{n}\left(l\right){e}_{n}^{\prime }\left(l\right)\right\}$ where en(l) = Zn + l * n * (l)${e}_{n}\left(l\right)={Z}_{n+l}^{*}-{\stackrel{^}{Z}}_{n}^{*}\left(l\right)$ is the forecast error. n * (l)${\stackrel{^}{Z}}_{n}^{*}\left(l\right)$ is referred to as the linear minimum mean square error forecast of Zn + l * ${Z}_{n+l}^{*}$.
The linear predictor which minimizes the mean square error may be expressed as
 Ẑn * (l) = τ + ψlεn + ψl + 1εn − 1 + ψl + 2εn − 2 + ⋯ . $Z^n*(l)=τ+ψlεn+ψl+1εn-1+ψl+2εn-2+⋯.$
The forecast error at t$t$ for lead l$l$ is then
 en(l) = Zn + l * − Ẑn * (l) = εn + l + ψ1εn + l − 1 + ψ2εn + l − 2 + ⋯ + ψl − 1εn + 1. $en(l)=Zn+l*-Z^n*(l)=εn+l+ψ1εn+l-1+ψ2εn+l-2+⋯+ψl-1εn+1.$
Let d = max (di)$\mathit{d}=\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathit{d}}_{\mathit{i}}\right)$, for i = 1,2,,k$\mathit{i}=1,2,\dots ,k$. Unless q = 0$q=0$ the function requires estimates of εt${\epsilon }_{\mathit{t}}$, for t = d + 1,,n$\mathit{t}=\mathit{d}+1,\dots ,n$, which are obtainable from nag_tsa_multi_varma_estimate (g13dd). The terms εt${\epsilon }_{\mathit{t}}$ are assumed to be zero, for t = n + 1,,n + lmax$\mathit{t}=n+1,\dots ,n+{l}_{\mathrm{max}}$. You may use nag_tsa_multi_varma_update (g13dk) to update these lmax${l}_{\mathrm{max}}$ forecasts should further observations, Zn + 1,Zn + 2,${Z}_{n+1},{Z}_{n+2},\dots \text{}$, become available. Note that when lmax${l}_{\mathrm{max}}$ or more further observations are available then nag_tsa_multi_varma_forecast (g13dj) must be used to produce new forecasts for Zn + lmax + 1,Zn + lmax + 2,${Z}_{n+{l}_{\mathrm{max}}+1},{Z}_{n+{l}_{\mathrm{max}}+2},\dots \text{}$, should they be required.
When a transformation has been used the forecasts and their standard errors are suitably modified to give results in terms of the original series, Zt${Z}_{t}$; see Granger and Newbold (1976).

## References

Box G E P and Jenkins G M (1976) Time Series Analysis: Forecasting and Control (Revised Edition) Holden–Day
Granger C W J and Newbold P (1976) Forecasting transformed series J. Roy. Statist. Soc. Ser. B 38 189–203
Wei W W S (1990) Time Series Analysis: Univariate and Multivariate Methods Addison–Wesley

## Parameters

The output quantities k, n, kmax, ip, iq, par, npar, qq and v from nag_tsa_multi_varma_estimate (g13dd) are suitable for input to nag_tsa_multi_varma_forecast (g13dj).

### Compulsory Input Parameters

1:     z(kmax,n) – double array
kmax, the first dimension of the array, must satisfy the constraint kmaxk$\mathit{kmax}\ge {\mathbf{k}}$.
z(i,t)${\mathbf{z}}\left(\mathit{i},\mathit{t}\right)$ must contain, zit${z}_{\mathit{i}\mathit{t}}$, the i$\mathit{i}$th component of Zt${Z}_{\mathit{t}}$, for i = 1,2,,k$\mathit{i}=1,2,\dots ,k$ and t = 1,2,,n$\mathit{t}=1,2,\dots ,n$.
Constraints:
• if tr(i) = 'L'${\mathbf{tr}}\left(i\right)=\text{'L'}$, z(i,t) > 0.0${\mathbf{z}}\left(i,t\right)>0.0$;
• if tr(i) = 'S'${\mathbf{tr}}\left(i\right)=\text{'S'}$, z(i,t)0.0${\mathbf{z}}\left(\mathit{i},\mathit{t}\right)\ge 0.0$, for i = 1,2,,k$\mathit{i}=1,2,\dots ,k$ and t = 1,2,,n$\mathit{t}=1,2,\dots ,n$.
2:     tr(k) – cell array of strings
k, the dimension of the array, must satisfy the constraint k1${\mathbf{k}}\ge 1$.
tr(i)${\mathbf{tr}}\left(\mathit{i}\right)$ indicates whether the i$\mathit{i}$th time series is to be transformed, for i = 1,2,,k$\mathit{i}=1,2,\dots ,k$.
tr(i) = 'N'${\mathbf{tr}}\left(i\right)=\text{'N'}$
No transformation is used.
tr(i) = 'L'${\mathbf{tr}}\left(i\right)=\text{'L'}$
A log transformation is used.
tr(i) = 'S'${\mathbf{tr}}\left(i\right)=\text{'S'}$
A square root transformation is used.
Constraint: tr(i) = 'N'${\mathbf{tr}}\left(\mathit{i}\right)=\text{'N'}$, 'L'$\text{'L'}$ or 'S'$\text{'S'}$, for i = 1,2,,k$\mathit{i}=1,2,\dots ,k$.
3:     id(k) – int64int32nag_int array
k, the dimension of the array, must satisfy the constraint k1${\mathbf{k}}\ge 1$.
id(i)${\mathbf{id}}\left(i\right)$ must specify, di${\mathit{d}}_{i}$, the order of differencing required for the i$i$th series.
Constraint: 0id(i) < nmax (ip,iq)$0\le {\mathbf{id}}\left(\mathit{i}\right)<{\mathbf{n}}-\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{ip}},{\mathbf{iq}}\right)$, for i = 1,2,,k$\mathit{i}=1,2,\dots ,k$.
4:     delta(kmax, : $:$) – double array
The first dimension of the array delta must be at least k${\mathbf{k}}$
The second dimension of the array must be at least max (1,d)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,\mathit{d}\right)$, where d = max (id(i))$\mathit{d}=\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{id}}\left(i\right)\right)$
If id(i) > 0${\mathbf{id}}\left(i\right)>0$, then delta(i,j)${\mathbf{delta}}\left(\mathit{i},\mathit{j}\right)$ must be set equal to δij${\delta }_{\mathit{i}\mathit{j}}$, for j = 1,2,,di$\mathit{j}=1,2,\dots ,{\mathit{d}}_{\mathit{i}}$ and i = 1,2,,k$\mathit{i}=1,2,\dots ,k$.
If d = 0$\mathit{d}=0$, delta is not referenced.
5:     ip – int64int32nag_int scalar
p$p$, the number of AR parameter matrices.
Constraint: ip0${\mathbf{ip}}\ge 0$.
6:     iq – int64int32nag_int scalar
q$q$, the number of MA parameter matrices.
Constraint: iq0${\mathbf{iq}}\ge 0$.
7:     mean – string (length ≥ 1)
mean = 'M'${\mathbf{mean}}=\text{'M'}$, if components of μ$\mu$ have been estimated and mean = 'Z'${\mathbf{mean}}=\text{'Z'}$, if all elements of μ$\mu$ are to be taken as zero.
Constraint: mean = 'M'${\mathbf{mean}}=\text{'M'}$ or 'Z'$\text{'Z'}$.
8:     par(lpar) – double array
lpar, the dimension of the array, must satisfy the constraint
• if mean = 'Z'${\mathbf{mean}}=\text{'Z'}$, lparmax (1,(ip + iq) × k × k)${\mathbf{lpar}}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,\left({\mathbf{ip}}+{\mathbf{iq}}\right)×{\mathbf{k}}×{\mathbf{k}}\right)$;
• if mean = 'M'${\mathbf{mean}}=\text{'M'}$, lpar(ip + iq) × k × k + k${\mathbf{lpar}}\ge \left({\mathbf{ip}}+{\mathbf{iq}}\right)×{\mathbf{k}}×{\mathbf{k}}+{\mathbf{k}}$.
Must contain the parameter estimates read in row by row in the order φ1,φ2,,φp${\varphi }_{1},{\varphi }_{2},\dots ,{\varphi }_{p}$, θ1,θ2,,θq${\theta }_{1},{\theta }_{2},\dots ,{\theta }_{q}$, μ$\mu$.
Thus,
• if ip > 0${\mathbf{ip}}>0$, par((l1) × k × k + (i1) × k + j)${\mathbf{par}}\left(\left(\mathit{l}-1\right)×k×k+\left(\mathit{i}-1\right)×k+\mathit{j}\right)$ must be set equal to an estimate of the (i,j)$\left(\mathit{i},\mathit{j}\right)$th element of φl${\varphi }_{\mathit{l}}$, for l = 1,2,,p$\mathit{l}=1,2,\dots ,p$, i = 1,2,,k$\mathit{i}=1,2,\dots ,k$ and j = 1,2,,k$\mathit{j}=1,2,\dots ,k$;
• if iq > 0${\mathbf{iq}}>0$, par(p × k × k + (l1) × k × k + (i1) × k + j)${\mathbf{par}}\left(p×k×k+\left(\mathit{l}-1\right)×k×k+\left(\mathit{i}-1\right)×k+\mathit{j}\right)$ must be set equal to an estimate of 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$;
• if mean = 'M'${\mathbf{mean}}=\text{'M'}$, par((p + q) × k × k + i)${\mathbf{par}}\left(\left(p+q\right)×k×k+\mathit{i}\right)$ must be set equal to an estimate of the i$\mathit{i}$th component of μ$\mu$, for i = 1,2,,k$\mathit{i}=1,2,\dots ,k$.
Constraint: the first ip × k × k${\mathbf{ip}}×{\mathbf{k}}×{\mathbf{k}}$ elements of par must satisfy the stationarity condition and the next iq × k × k${\mathbf{iq}}×{\mathbf{k}}×{\mathbf{k}}$ elements of par must satisfy the invertibility condition.
9:     qq(kmax,k) – double array
kmax, the first dimension of the array, must satisfy the constraint kmaxk$\mathit{kmax}\ge {\mathbf{k}}$.
qq(i,j)${\mathbf{qq}}\left(i,j\right)$ must contain an estimate of the (i,j)$\left(i,j\right)$th element of Σ$\Sigma$. The lower triangle only is needed.
Constraint: qq${\mathbf{qq}}$ must be positive definite.
10:   v(kmax, : $:$) – double array
The first dimension of the array v must be at least k${\mathbf{k}}$
The second dimension of the array must be at least max (1,nd)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}-\mathit{d}\right)$, where d = max (id(i))$\mathit{d}=\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{id}}\left(i\right)\right)$
v(i,t)${\mathbf{v}}\left(\mathit{i},\mathit{t}\right)$ must contain an estimate of the i$\mathit{i}$th component of εt + d${\epsilon }_{\mathit{t}+\mathit{d}}$, for i = 1,2,,k$\mathit{i}=1,2,\dots ,k$ and t = 1,2,,nd$\mathit{t}=1,2,\dots ,n-\mathit{d}$.
If q = 0$q=0$, v is not used.
11:   lmax – int64int32nag_int scalar
The number, lmax${l}_{\mathrm{max}}$, of forecasts required.
Constraint: lmax1${\mathbf{lmax}}\ge 1$.
12:   lref – int64int32nag_int scalar
The dimension of the array ref as declared in the (sub)program from which nag_tsa_multi_varma_forecast (g13dj) is called.
Constraint: lref(lmax1) × k × k + 2 × k × lmax + k${\mathbf{lref}}\ge \left({\mathbf{lmax}}-1\right)×{\mathbf{k}}×{\mathbf{k}}+2×{\mathbf{k}}×{\mathbf{lmax}}+{\mathbf{k}}$.

### Optional Input Parameters

1:     k – int64int32nag_int scalar
Default: The dimension of the arrays tr, id and the first dimension of the arrays z, delta, qq, v and the second dimension of the array qq. (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:     n – int64int32nag_int scalar
Default: The second dimension of the array z.
n$n$, the number of observations in the series, Zt${Z}_{t}$, prior to differencing.
Constraint: n3${\mathbf{n}}\ge 3$.
The total number of observations must exceed the total number of parameters in the model; that is
• if mean = 'Z'${\mathbf{mean}}=\text{'Z'}$, n × k > (ip + iq) × k × k + k × (k + 1) / 2${\mathbf{n}}×{\mathbf{k}}>\left({\mathbf{ip}}+{\mathbf{iq}}\right)×{\mathbf{k}}×{\mathbf{k}}+{\mathbf{k}}×\left({\mathbf{k}}+1\right)/2$;
• if mean = 'M'${\mathbf{mean}}=\text{'M'}$, n × k > (ip + iq) × k × k + k + k × (k + 1) / 2${\mathbf{n}}×{\mathbf{k}}>\left({\mathbf{ip}}+{\mathbf{iq}}\right)×{\mathbf{k}}×{\mathbf{k}}+{\mathbf{k}}+{\mathbf{k}}×\left({\mathbf{k}}+1\right)/2$,
(see the parameters ip, iq and mean).
3:     lpar – int64int32nag_int scalar
Default: The dimension of the array par.
The dimension of the array par as declared in the (sub)program from which nag_tsa_multi_varma_forecast (g13dj) is called.
Constraints:
• if mean = 'Z'${\mathbf{mean}}=\text{'Z'}$, lparmax (1,(ip + iq) × k × k)${\mathbf{lpar}}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,\left({\mathbf{ip}}+{\mathbf{iq}}\right)×{\mathbf{k}}×{\mathbf{k}}\right)$;
• if mean = 'M'${\mathbf{mean}}=\text{'M'}$, lpar(ip + iq) × k × k + k${\mathbf{lpar}}\ge \left({\mathbf{ip}}+{\mathbf{iq}}\right)×{\mathbf{k}}×{\mathbf{k}}+{\mathbf{k}}$.

### Input Parameters Omitted from the MATLAB Interface

kmax work lwork iwork liwork

### Output Parameters

1:     qq(kmax,k) – double array
kmaxk$\mathit{kmax}\ge {\mathbf{k}}$.
If ${\mathbf{ifail}}\ne {\mathbf{1}}$, then the upper triangle is set equal to the lower triangle.
2:     predz(kmax,lmax) – double array
kmaxk$\mathit{kmax}\ge {\mathbf{k}}$.
predz(i,l)${\mathbf{predz}}\left(\mathit{i},\mathit{l}\right)$ contains the forecast of zi,n + l${z}_{\mathit{i},n+\mathit{l}}$, for i = 1,2,,k$\mathit{i}=1,2,\dots ,k$ and l = 1,2,,lmax$\mathit{l}=1,2,\dots ,{\mathit{l}}_{\mathrm{max}}$.
3:     sefz(kmax,lmax) – double array
kmaxk$\mathit{kmax}\ge {\mathbf{k}}$.
sefz(i,l)${\mathbf{sefz}}\left(\mathit{i},\mathit{l}\right)$ contains an estimate of the standard error of the forecast of zi,n + l${z}_{\mathit{i},n+\mathit{l}}$, for i = 1,2,,k$\mathit{i}=1,2,\dots ,k$ and l = 1,2,,lmax$\mathit{l}=1,2,\dots ,{\mathit{l}}_{\mathrm{max}}$.
4:     ref(lref) – double array
The reference vector which may be used to update forecasts using nag_tsa_multi_varma_update (g13dk). The first (lmax1) × k × k$\left({\mathbf{lmax}}-1\right)×{\mathbf{k}}×{\mathbf{k}}$ elements contain the ψ$\psi$ weight matrices, ψ1,ψ2,,ψlmax1${\psi }_{1},{\psi }_{2},\dots ,{\psi }_{{l}_{\mathrm{max}}-1}$. The next ${\mathbf{k}}×{\mathbf{lmax}}$ elements contain the forecasts of the transformed series n + 1 * ,n + 2 * ,, n + lmax * ${\stackrel{^}{Z}}_{n+1}^{*},{\stackrel{^}{Z}}_{n+2}^{*},\dots ,{\stackrel{^}{Z}}_{n+{l}_{\mathrm{max}}}^{*}$ and the next ${\mathbf{k}}×{\mathbf{lmax}}$ contain the variances of the forecasts of the transformed variables. The last k elements are used to store the transformations for the series.
5:     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:
ifail = 1${\mathbf{ifail}}=1$
 On entry, k < 1${\mathbf{k}}<1$, or n < 3${\mathbf{n}}<3$, or kmax < k$\mathit{kmax}<{\mathbf{k}}$, or id(i) < 0${\mathbf{id}}\left(i\right)<0$ for some i = 1,2, … ,k$i=1,2,\dots ,k$, or id(i) ≥ n − max (ip,iq)${\mathbf{id}}\left(i\right)\ge {\mathbf{n}}-\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{ip}},{\mathbf{iq}}\right)$ for some i = 1,2, … ,k$i=1,2,\dots ,k$, or ip < 0${\mathbf{ip}}<0$, or iq < 0${\mathbf{iq}}<0$, or mean ≠ 'M'${\mathbf{mean}}\ne \text{'M'}$ or 'Z'$\text{'Z'}$, or lpar < (ip + iq) × k × k + k${\mathbf{lpar}}<\left({\mathbf{ip}}+{\mathbf{iq}}\right)×{\mathbf{k}}×{\mathbf{k}}+{\mathbf{k}}$, and mean = 'M'${\mathbf{mean}}=\text{'M'}$, or lpar < (ip + iq) × k × k${\mathbf{lpar}}<\left({\mathbf{ip}}+{\mathbf{iq}}\right)×{\mathbf{k}}×{\mathbf{k}}$ and mean = 'Z'${\mathbf{mean}}=\text{'Z'}$, or n × k ≤ (ip + iq) × k × k + k + k(k + 1) / 2${\mathbf{n}}×{\mathbf{k}}\le \left({\mathbf{ip}}+{\mathbf{iq}}\right)×{\mathbf{k}}×{\mathbf{k}}+{\mathbf{k}}+{\mathbf{k}}\left({\mathbf{k}}+1\right)/2$, and mean = 'M'${\mathbf{mean}}=\text{'M'}$, or n × k ≤ (ip + iq) × k × k + k(k + 1) / 2${\mathbf{n}}×{\mathbf{k}}\le \left({\mathbf{ip}}+{\mathbf{iq}}\right)×{\mathbf{k}}×{\mathbf{k}}+{\mathbf{k}}\left({\mathbf{k}}+1\right)/2$ and mean = 'Z'${\mathbf{mean}}=\text{'Z'}$, or lmax < 1${\mathbf{lmax}}<1$, or lref < (lmax − 1) × k × k + 2 × k × lmax + k${\mathbf{lref}}<\left({\mathbf{lmax}}-1\right)×{\mathbf{k}}×{\mathbf{k}}+2×{\mathbf{k}}×{\mathbf{lmax}}+{\mathbf{k}}$, or lwork is too small, or liwork is too small.
ifail = 2${\mathbf{ifail}}=2$
 On entry, at least one of the first k$k$ elements of tr is not equal to 'N', 'L' or 'S'.
ifail = 3${\mathbf{ifail}}=3$
On entry, one or more of the transformations requested cannot be computed; that is, you may be trying to log or square-root a series, some of whose values are negative.
ifail = 4${\mathbf{ifail}}=4$
On entry, either qq is not positive definite or the autoregressive parameter matrices are extremely close to or outside the stationarity region, or the moving average parameter matrices are extremely close to or outside the invertibility region. To proceed, you must supply different parameter estimates in the arrays par and qq.
ifail = 5${\mathbf{ifail}}=5$
This is an unlikely exit brought about by an excessive number of iterations being needed to evaluate the eigenvalues of the matrices required to check for stationarity and invertibility; see nag_tsa_uni_arma_roots (g13dx). All output parameters are undefined.
ifail = 6${\mathbf{ifail}}=6$
This is an unlikely exit which could occur if qq is nearly non positive definite. In this case the standard deviations of the forecast errors may be non-positive. To proceed, you must supply different parameter estimates in the array qq.
ifail = 7${\mathbf{ifail}}=7$
This is an unlikely exit. For one of the series, overflow will occur if the forecasts are computed. You should check whether the transformations requested in the array tr are sensible. All output parameters are undefined.

## Accuracy

The matrix computations are believed to be stable.

The same differencing operator does not have to be applied to all the series. For example, suppose we have k = 2$k=2$, and wish to apply the second order differencing operator 2${\nabla }^{2}$ to the first series and the first-order differencing operator $\nabla$ to the second series:
 w1t = ∇2z1t = (1 − B)2z1t = (1 − 2B + B2)Z1t,   and w2t = ∇z2t = (1 − B)z2t.
$w1t=∇2z1t= (1-B) 2z1t=(1-2B+B2)Z1t, and w2t=∇z2t=(1-B)z2t.$
Then d1 = 2,d2 = 1${\mathit{d}}_{1}=2,{\mathit{d}}_{2}=1$, d = max (d1,d2) = 2$\mathit{d}=\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathit{d}}_{1},{\mathit{d}}_{2}\right)=2$, and
delta =
 [ δ11 δ12 δ21 ]
=
 [ 2 − 1 1 ]
.
$delta=[ δ11 δ12 δ21 ]=[ 2 -1 1 ] .$
Note:  although differencing may already have been applied prior to the model fitting stage, the differencing parameters supplied in delta are part of the model definition and are still required by this function to produce the forecasts.
nag_tsa_multi_varma_forecast (g13dj) should not be used when the moving average parameters lie close to the boundary of the invertibility region. The function does test for both invertibility and stationarity but if in doubt, you may use nag_tsa_uni_arma_roots (g13dx), before calling this function, to check that the VARMA model being used is invertible.
On a successful exit, the output quantities k, lmax, kmax, ref and lref will be suitable for input to nag_tsa_multi_varma_update (g13dk).

## Example

```function nag_tsa_multi_varma_forecast_example
z = [-1.49, -1.62, 5.2, 6.23, 6.21, 5.86, 4.09, 3.18, 2.62, 1.49, 1.17, ...
0.85, -0.35, 0.24, 2.44, 2.58, 2.04, 0.4, 2.26, 3.34, 5.09, 5, 4.78, ...
4.11, 3.45, 1.65, 1.29, 4.09, 6.32, 7.5, 3.89, 1.58, 5.21, 5.25, 4.93, ...
7.38, 5.87, 5.81, 9.68, 9.07, 7.29, 7.84, 7.55, 7.32, 7.97, 7.76, 7, 8.35;
7.34, 6.35, 6.96, 8.54, 6.62, 4.97, 4.55, 4.81, 4.75, ...
4.76, 10.88, 10.01, 11.62, 10.36, 6.4, 6.24, 7.93, 4.04, 3.73, 5.6, ...
5.35, 6.81, 8.27, 7.68, 6.65, 6.08, 10.25, 9.14, 17.75, 13.3, ...
9.63, 6.8, 4.08, 5.06, 4.94, 6.65, 7.94, 10.76, 11.89, ...
5.85, 9.01, 7.5, 10.02, 10.38, 8.15, 8.37, 10.73, 12.14];
tr = {'N'; 'N'};
id = [int64(0);0];
delta = [0;
0];
ip = int64(1);
iq = int64(0);
mean_p = 'M';
par = [0.8016071892386086;
0.0648134906597352;
0;
0.575015951133362;
4.271122828253269;
7.825342792089621];
qq = [2.964154253391392, 0.6372583252520638;
0.6372583252520638, 5.379903126133676];
v = [-3.326133715726029, -1.241508590516912, 5.746865699837245, 1.270368439927496, ...
0.3223077197693467, 0.112781765620811, -1.269713458557113, -0.7336470675276869, ...
-0.5810360328920845, -1.25824719747888, -0.673079208545849, -1.133223470827074, ...
-2.02032143339675, -0.572742526272592, 1.235974230307894, -0.1309001630044928, ...
-0.7727550109923405, -2.089421928018445, 1.337338340999243, 0.9464411511199494, ...
1.709504159208547, 0.2328949507059163, -0.009588098625822217, -0.5978622133565419, ...
-0.6825454370774304, -1.886726796800422, -0.7668901664948766, ...
2.051416165579926, 2.108859010344129, 0.9432308237617111, -3.324245626104025, ...
-2.50257816223142, 3.162556623476816, 0.4690152211351468, 0.05343371271906172, ...
2.767725632154585, -0.8170430505081534, 0.2497744022910866, 3.985096789984951, ...
0.1996377231860347, -0.6999084077936142, 1.072141758566346, 0.4391261753813116, ...
0.2782622637979761, 1.08929906068535, 0.5027884718514648, -0.103132986353569, 1.703128639510798;
-0.1864812965451277, -1.196262944870486, -0.01699715324845752, 1.212243116560191, ...
-1.61628208623052, -2.162251460054466, -1.633475140684418, -1.131968441208406, ...
-1.34147258850308, -1.296971631435079, 4.817278209053589, 0.4281805881174119, ...
2.538444465603437, 0.3526687842787246, -2.882811117293238, -0.7657479508051256, ...
1.016254601376212, -3.84552235603917, -1.918710306130392, 0.1295446387209505, ...
-1.195735189898436, 0.4080187978849044, 1.028495509230196, -0.4010277794245127, ...
-1.091768368255828, -1.069501938588466, 3.42825715355755, -0.07955936266856867, ...
9.168708343089461, -0.2321789961687832, -1.343358013625323, -2.063049472965885, ...
-3.15575433125847, -0.6117109441757265, -1.29522657628642, 0.4837753378495833, ...
0.7904980614115339, 2.868727484449496, 2.377182502253417, -4.312585522527284, ...
2.320510822318223, -1.0065395832632, 2.381734502948175, 1.292694306092105, ...
-1.144311436315907, 0.3579741347114898, 2.591470625462152, 2.644432980787418];
lmax = int64(5);
lref = int64(150);
[qqOut, predz, sefz, ref, ifail] = ...
nag_tsa_multi_varma_forecast(z, tr, id, delta, ip, iq, mean_p, par, qq, v, lmax, lref)
```
```

qqOut =

2.9642    0.6373
0.6373    5.3799

predz =

7.8204    7.2771    6.7732    6.3300    5.9521
10.3063    9.2520    8.6457    8.2970    8.0966

sefz =

1.7217    2.2266    2.5095    2.6817    2.7898
2.3195    2.6756    2.7833    2.8180    2.8294

ref =

0.8016
0
0.0648
0.5750
0.6426
0
0.0892
0.3306
0.5151
0
0.0930
0.1901
0.4129
0
0.0868
0.1093
7.8204
10.3063
7.2771
9.2520
6.7732
8.6457
6.3300
8.2970
5.9521
8.0966
2.9642
5.3799
4.9577
7.1587
6.2975
7.7469
7.1914
7.9414
7.7830
8.0057
100.0000
100.0000
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0

ifail =

0

```
```function g13dj_example
z = [-1.49, -1.62, 5.2, 6.23, 6.21, 5.86, 4.09, 3.18, 2.62, 1.49, 1.17, ...
0.85, -0.35, 0.24, 2.44, 2.58, 2.04, 0.4, 2.26, 3.34, 5.09, 5, 4.78, ...
4.11, 3.45, 1.65, 1.29, 4.09, 6.32, 7.5, 3.89, 1.58, 5.21, 5.25, 4.93, ...
7.38, 5.87, 5.81, 9.68, 9.07, 7.29, 7.84, 7.55, 7.32, 7.97, 7.76, 7, 8.35;
7.34, 6.35, 6.96, 8.54, 6.62, 4.97, 4.55, 4.81, 4.75, ...
4.76, 10.88, 10.01, 11.62, 10.36, 6.4, 6.24, 7.93, 4.04, 3.73, 5.6, ...
5.35, 6.81, 8.27, 7.68, 6.65, 6.08, 10.25, 9.14, 17.75, 13.3, ...
9.63, 6.8, 4.08, 5.06, 4.94, 6.65, 7.94, 10.76, 11.89, ...
5.85, 9.01, 7.5, 10.02, 10.38, 8.15, 8.37, 10.73, 12.14];
tr = {'N'; 'N'};
id = [int64(0);0];
delta = [0;
0];
ip = int64(1);
iq = int64(0);
mean_p = 'M';
par = [0.8016071892386086;
0.0648134906597352;
0;
0.575015951133362;
4.271122828253269;
7.825342792089621];
qq = [2.964154253391392, 0.6372583252520638;
0.6372583252520638, 5.379903126133676];
v = [-3.326133715726029, -1.241508590516912, 5.746865699837245, 1.270368439927496, ...
0.3223077197693467, 0.112781765620811, -1.269713458557113, -0.7336470675276869, ...
-0.5810360328920845, -1.25824719747888, -0.673079208545849, -1.133223470827074, ...
-2.02032143339675, -0.572742526272592, 1.235974230307894, -0.1309001630044928, ...
-0.7727550109923405, -2.089421928018445, 1.337338340999243, 0.9464411511199494, ...
1.709504159208547, 0.2328949507059163, -0.009588098625822217, -0.5978622133565419, ...
-0.6825454370774304, -1.886726796800422, -0.7668901664948766, ...
2.051416165579926, 2.108859010344129, 0.9432308237617111, -3.324245626104025, ...
-2.50257816223142, 3.162556623476816, 0.4690152211351468, 0.05343371271906172, ...
2.767725632154585, -0.8170430505081534, 0.2497744022910866, 3.985096789984951, ...
0.1996377231860347, -0.6999084077936142, 1.072141758566346, 0.4391261753813116, ...
0.2782622637979761, 1.08929906068535, 0.5027884718514648, -0.103132986353569, 1.703128639510798;
-0.1864812965451277, -1.196262944870486, -0.01699715324845752, 1.212243116560191, ...
-1.61628208623052, -2.162251460054466, -1.633475140684418, -1.131968441208406, ...
-1.34147258850308, -1.296971631435079, 4.817278209053589, 0.4281805881174119, ...
2.538444465603437, 0.3526687842787246, -2.882811117293238, -0.7657479508051256, ...
1.016254601376212, -3.84552235603917, -1.918710306130392, 0.1295446387209505, ...
-1.195735189898436, 0.4080187978849044, 1.028495509230196, -0.4010277794245127, ...
-1.091768368255828, -1.069501938588466, 3.42825715355755, -0.07955936266856867, ...
9.168708343089461, -0.2321789961687832, -1.343358013625323, -2.063049472965885, ...
-3.15575433125847, -0.6117109441757265, -1.29522657628642, 0.4837753378495833, ...
0.7904980614115339, 2.868727484449496, 2.377182502253417, -4.312585522527284, ...
2.320510822318223, -1.0065395832632, 2.381734502948175, 1.292694306092105, ...
-1.144311436315907, 0.3579741347114898, 2.591470625462152, 2.644432980787418];
lmax = int64(5);
lref = int64(150);
[qqOut, predz, sefz, ref, ifail] = ...
g13dj(z, tr, id, delta, ip, iq, mean_p, par, qq, v, lmax, lref)
```
```

qqOut =

2.9642    0.6373
0.6373    5.3799

predz =

7.8204    7.2771    6.7732    6.3300    5.9521
10.3063    9.2520    8.6457    8.2970    8.0966

sefz =

1.7217    2.2266    2.5095    2.6817    2.7898
2.3195    2.6756    2.7833    2.8180    2.8294

ref =

0.8016
0
0.0648
0.5750
0.6426
0
0.0892
0.3306
0.5151
0
0.0930
0.1901
0.4129
0
0.0868
0.1093
7.8204
10.3063
7.2771
9.2520
6.7732
8.6457
6.3300
8.2970
5.9521
8.0966
2.9642
5.3799
4.9577
7.1587
6.2975
7.7469
7.1914
7.9414
7.7830
8.0057
100.0000
100.0000
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0

ifail =

0

```