# NAG FL Interfaceg13djf (multi_​varma_​forecast)

## ▸▿ Contents

Settings help

FL Name Style:

FL Specification Language:

## 1Purpose

g13djf 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 g13ddf. 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 g13dkf.

## 2Specification

Fortran Interface
 Subroutine g13djf ( k, n, z, kmax, tr, id, ip, iq, mean, par, lpar, qq, v, lmax, sefz, ref, lref, work,
 Integer, Intent (In) :: k, n, kmax, id(k), ip, iq, lpar, lmax, lref, lwork, liwork Integer, Intent (Inout) :: ifail Integer, Intent (Out) :: iwork(liwork) Real (Kind=nag_wp), Intent (In) :: z(kmax,n), delta(kmax,*), par(lpar), v(kmax,*) Real (Kind=nag_wp), Intent (Inout) :: qq(kmax,k), predz(kmax,lmax), sefz(kmax,lmax) Real (Kind=nag_wp), Intent (Out) :: ref(lref), work(lwork) Character (1), Intent (In) :: tr(k), mean
#include <nag.h>
 void g13djf_ (const Integer *k, const Integer *n, const double z[], const Integer *kmax, const char tr[], const Integer id[], const double delta[], const Integer *ip, const Integer *iq, const char *mean, const double par[], const Integer *lpar, double qq[], const double v[], const Integer *lmax, double predz[], double sefz[], double ref[], const Integer *lref, double work[], const Integer *lwork, Integer iwork[], const Integer *liwork, Integer *ifail, const Charlen length_tr, const Charlen length_mean)
The routine may be called by the names g13djf or nagf_tsa_multi_varma_forecast.

## 3Description

Let the vector ${Z}_{\mathit{t}}={\left({z}_{1\mathit{t}},{z}_{2\mathit{t}},\dots ,{z}_{k\mathit{t}}\right)}^{\mathrm{T}}$, for $\mathit{t}=1,2,\dots ,n$, denote a $k$-dimensional time series for which forecasts of ${Z}_{n+1},{Z}_{n+2},\dots ,{Z}_{n+{l}_{\mathrm{max}}}$ are required. Let ${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,$
where ${\delta }_{i}\left(B\right)$ is the differencing operator applied to the $i$th series and where ${z}_{it}^{*}$ is equal to either ${z}_{it}$, $\sqrt{{z}_{it}}$ or ${\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$th series is ${\mathit{d}}_{i}$, then the differencing operator for the $i$th series is defined by ${\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$ is the backward shift operator; that is, $B{Z}_{t}={Z}_{t-1}$. The differencing parameters ${\delta }_{\mathit{i}\mathit{j}}$, for $\mathit{i}=1,2,\dots ,k$ and $\mathit{j}=1,2,\dots ,{\mathit{d}}_{\mathit{i}}$, must be supplied by you. If the $i$th series does not require differencing, then ${\mathit{d}}_{i}=0$.
${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,$ (1)
where ${\epsilon }_{\mathit{t}}={\left({\epsilon }_{1\mathit{t}},{\epsilon }_{2\mathit{t}},\dots ,{\epsilon }_{k\mathit{t}}\right)}^{\mathrm{T}}$, for $\mathit{t}=1,2,\dots ,n$, is a vector of $k$ residual series assumed to be Normally distributed with zero mean and positive definite covariance matrix $\Sigma$. The components of ${\epsilon }_{t}$ are assumed to be uncorrelated at non-simultaneous lags. The ${\varphi }_{i}$ and ${\theta }_{j}$ are $k×k$ matrices of parameters. The matrices ${\varphi }_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,p$, are the autoregressive (AR) parameter matrices, and the matrices ${\theta }_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,q$, the moving average (MA) parameter matrices. The parameters in the model are thus the $p$ ($k×k$) $\varphi$-matrices, the $q$ ($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 g13dxf for a method of checking these conditions.
The ARMA model (1) may be rewritten as
 $ϕ(B)(δ(B)Zt*-μ)=θ(B)εt,$
where $\varphi \left(B\right)$ and $\theta \left(B\right)$ are the autoregressive and moving average polynomials and $\delta \left(B\right)$ denotes the $k×k$ diagonal matrix whose $i$th diagonal elements is ${\delta }_{i}\left(B\right)$ and ${Z}_{t}^{*}={\left({z}_{1t}^{*},{z}_{2t}^{*}\dots {z}_{kt}^{*}\right)}^{\mathrm{T}}$.
This may be rewritten as
 $ϕ(B)δ(B)Zt*=ϕ(B)μ+θ(B)εt$
or
 $Zt*=τ+ψ (B)εt=τ+εt+ψ1εt- 1+ψ2εt- 2+⋯$
where $\psi \left(B\right)={\delta }^{-1}\left(B\right){\varphi }^{-1}\left(B\right)\theta \left(B\right)$ and $\tau ={\delta }^{-1}\left(B\right)\mu$ is a vector of length $k$.
Forecasts are computed using a multivariate version of the procedure described in Box and Jenkins (1976). If ${\stackrel{^}{Z}}_{n}^{*}\left(l\right)$ denotes the forecast of ${Z}_{n+l}^{*}$, then ${\stackrel{^}{Z}}_{n}^{*}\left(l\right)$ is taken to be that linear function of ${Z}_{n}^{*},{Z}_{n-1}^{*},\dots \text{}$ which minimizes the elements of $E\left\{{e}_{n}\left(l\right){e}_{n}^{\prime }\left(l\right)\right\}$ where ${e}_{n}\left(l\right)={Z}_{n+l}^{*}-{\stackrel{^}{Z}}_{n}^{*}\left(l\right)$ is the forecast error. ${\stackrel{^}{Z}}_{n}^{*}\left(l\right)$ is referred to as the linear minimum mean square error forecast of ${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+⋯.$
The forecast error at $t$ for lead $l$ is then
 $en(l)=Zn+l*-Z^n*(l)=εn+l+ψ1εn+l-1+ψ2εn+l-2+⋯+ψl-1εn+1.$
Let $\mathit{d}=\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathit{d}}_{\mathit{i}}\right)$, for $\mathit{i}=1,2,\dots ,k$. Unless $q=0$ the routine requires estimates of ${\epsilon }_{\mathit{t}}$, for $\mathit{t}=\mathit{d}+1,\dots ,n$, which are obtainable from g13ddf. The terms ${\epsilon }_{\mathit{t}}$ are assumed to be zero, for $\mathit{t}=n+1,\dots ,n+{l}_{\mathrm{max}}$. You may use g13dkf to update these ${l}_{\mathrm{max}}$ forecasts should further observations, ${Z}_{n+1},{Z}_{n+2},\dots \text{}$, become available. Note that when ${l}_{\mathrm{max}}$ or more further observations are available then g13djf must be used to produce new forecasts for ${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, ${Z}_{t}$; see Granger and Newbold (1976).

## 4References

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

## 5Arguments

The quantities k, n, kmax, ip, iq, par, npar, qq and v from g13ddf are suitable for input to g13djf.
1: $\mathbf{k}$Integer Input
On entry: $k$, the dimension of the multivariate time series.
Constraint: ${\mathbf{k}}\ge 1$.
2: $\mathbf{n}$Integer Input
On entry: $n$, the number of observations in the series, ${Z}_{t}$, prior to differencing.
Constraint: ${\mathbf{n}}\ge 3$.
The total number of observations must exceed the total number of parameters in the model; that is
• if ${\mathbf{mean}}=\text{'Z'}$, ${\mathbf{n}}×{\mathbf{k}}>\left({\mathbf{ip}}+{\mathbf{iq}}\right)×{\mathbf{k}}×{\mathbf{k}}+{\mathbf{k}}×\left({\mathbf{k}}+1\right)/2$;
• if ${\mathbf{mean}}=\text{'M'}$, ${\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 arguments ip, iq and mean).
3: $\mathbf{z}\left({\mathbf{kmax}},{\mathbf{n}}\right)$Real (Kind=nag_wp) array Input
On entry: ${\mathbf{z}}\left(\mathit{i},\mathit{t}\right)$ must contain, ${z}_{\mathit{i}\mathit{t}}$, the $\mathit{i}$th component of ${Z}_{\mathit{t}}$, for $\mathit{i}=1,2,\dots ,k$ and $\mathit{t}=1,2,\dots ,n$.
Constraints:
• if ${\mathbf{tr}}\left(i\right)=\text{'L'}$, ${\mathbf{z}}\left(i,t\right)>0.0$;
• if ${\mathbf{tr}}\left(i\right)=\text{'S'}$, ${\mathbf{z}}\left(\mathit{i},\mathit{t}\right)\ge 0.0$, for $\mathit{i}=1,2,\dots ,k$ and $\mathit{t}=1,2,\dots ,n$.
4: $\mathbf{kmax}$Integer Input
On entry: the first dimension of the arrays z, delta, qq, v, predz and sefz as declared in the (sub)program from which g13djf is called.
Constraint: ${\mathbf{kmax}}\ge {\mathbf{k}}$.
5: $\mathbf{tr}\left({\mathbf{k}}\right)$Character(1) array Input
On entry: ${\mathbf{tr}}\left(\mathit{i}\right)$ indicates whether the $\mathit{i}$th time series is to be transformed, for $\mathit{i}=1,2,\dots ,k$.
${\mathbf{tr}}\left(i\right)=\text{'N'}$
No transformation is used.
${\mathbf{tr}}\left(i\right)=\text{'L'}$
A log transformation is used.
${\mathbf{tr}}\left(i\right)=\text{'S'}$
A square root transformation is used.
Constraint: ${\mathbf{tr}}\left(\mathit{i}\right)=\text{'N'}$, $\text{'L'}$ or $\text{'S'}$, for $\mathit{i}=1,2,\dots ,k$.
6: $\mathbf{id}\left({\mathbf{k}}\right)$Integer array Input
On entry: ${\mathbf{id}}\left(i\right)$ must specify, ${\mathit{d}}_{i}$, the order of differencing required for the $i$th series.
Constraint: $0\le {\mathbf{id}}\left(\mathit{i}\right)<{\mathbf{n}}-\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{ip}},{\mathbf{iq}}\right)$, for $\mathit{i}=1,2,\dots ,k$.
7: $\mathbf{delta}\left({\mathbf{kmax}},*\right)$Real (Kind=nag_wp) array Input
Note: the second dimension of the array delta must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,\mathit{d}\right)$, where $\mathit{d}=\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{id}}\left(i\right)\right)$.
On entry: if ${\mathbf{id}}\left(i\right)>0$, then ${\mathbf{delta}}\left(\mathit{i},\mathit{j}\right)$ must be set equal to ${\delta }_{\mathit{i}\mathit{j}}$, for $\mathit{j}=1,2,\dots ,{\mathit{d}}_{\mathit{i}}$ and $\mathit{i}=1,2,\dots ,k$.
If $\mathit{d}=0$, delta is not referenced.
8: $\mathbf{ip}$Integer Input
On entry: $p$, the number of AR parameter matrices.
Constraint: ${\mathbf{ip}}\ge 0$.
9: $\mathbf{iq}$Integer Input
On entry: $q$, the number of MA parameter matrices.
Constraint: ${\mathbf{iq}}\ge 0$.
10: $\mathbf{mean}$Character(1) Input
On entry: ${\mathbf{mean}}=\text{'M'}$, if components of $\mu$ have been estimated and ${\mathbf{mean}}=\text{'Z'}$, if all elements of $\mu$ are to be taken as zero.
Constraint: ${\mathbf{mean}}=\text{'M'}$ or $\text{'Z'}$.
11: $\mathbf{par}\left({\mathbf{lpar}}\right)$Real (Kind=nag_wp) array Input
On entry: must contain the parameter estimates read in row by row in the order ${\varphi }_{1},{\varphi }_{2},\dots ,{\varphi }_{p}$, ${\theta }_{1},{\theta }_{2},\dots ,{\theta }_{q}$, $\mu$.
Thus,
• if ${\mathbf{ip}}>0$, ${\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 $\left(\mathit{i},\mathit{j}\right)$th element of ${\varphi }_{\mathit{l}}$, for $\mathit{l}=1,2,\dots ,p$, $\mathit{i}=1,2,\dots ,k$ and $\mathit{j}=1,2,\dots ,k$;
• if ${\mathbf{iq}}>0$, ${\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 $\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$;
• if ${\mathbf{mean}}=\text{'M'}$, ${\mathbf{par}}\left(\left(p+q\right)×k×k+\mathit{i}\right)$ must be set equal to an estimate of the $\mathit{i}$th component of $\mu$, for $\mathit{i}=1,2,\dots ,k$.
Constraint: the first ${\mathbf{ip}}×{\mathbf{k}}×{\mathbf{k}}$ elements of par must satisfy the stationarity condition and the next ${\mathbf{iq}}×{\mathbf{k}}×{\mathbf{k}}$ elements of par must satisfy the invertibility condition.
12: $\mathbf{lpar}$Integer Input
On entry: the dimension of the array par as declared in the (sub)program from which g13djf is called.
Constraints:
• if ${\mathbf{mean}}=\text{'Z'}$, ${\mathbf{lpar}}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,\left({\mathbf{ip}}+{\mathbf{iq}}\right)×{\mathbf{k}}×{\mathbf{k}}\right)$;
• if ${\mathbf{mean}}=\text{'M'}$, ${\mathbf{lpar}}\ge \left({\mathbf{ip}}+{\mathbf{iq}}\right)×{\mathbf{k}}×{\mathbf{k}}+{\mathbf{k}}$.
13: $\mathbf{qq}\left({\mathbf{kmax}},{\mathbf{k}}\right)$Real (Kind=nag_wp) array Input/Output
On entry: ${\mathbf{qq}}\left(i,j\right)$ must contain an estimate of the $\left(i,j\right)$th element of $\Sigma$. The lower triangle only is needed.
Constraint: ${\mathbf{qq}}$ must be positive definite.
On exit: if ${\mathbf{ifail}}\ne {\mathbf{1}}$, then the upper triangle is set equal to the lower triangle.
14: $\mathbf{v}\left({\mathbf{kmax}},*\right)$Real (Kind=nag_wp) array Input
Note: the second dimension of the array v must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}-\mathit{d}\right)$, where $\mathit{d}=\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{id}}\left(i\right)\right)$.
On entry: ${\mathbf{v}}\left(\mathit{i},\mathit{t}\right)$ must contain an estimate of the $\mathit{i}$th component of ${\epsilon }_{\mathit{t}+\mathit{d}}$, for $\mathit{i}=1,2,\dots ,k$ and $\mathit{t}=1,2,\dots ,n-\mathit{d}$.
If $q=0$, v is not used.
15: $\mathbf{lmax}$Integer Input
On entry: the number, ${l}_{\mathrm{max}}$, of forecasts required.
Constraint: ${\mathbf{lmax}}\ge 1$.
16: $\mathbf{predz}\left({\mathbf{kmax}},{\mathbf{lmax}}\right)$Real (Kind=nag_wp) array Output
On exit: ${\mathbf{predz}}\left(\mathit{i},\mathit{l}\right)$ contains the forecast of ${z}_{\mathit{i},n+\mathit{l}}$, for $\mathit{i}=1,2,\dots ,k$ and $\mathit{l}=1,2,\dots ,{\mathit{l}}_{\mathrm{max}}$.
17: $\mathbf{sefz}\left({\mathbf{kmax}},{\mathbf{lmax}}\right)$Real (Kind=nag_wp) array Output
On exit: ${\mathbf{sefz}}\left(\mathit{i},\mathit{l}\right)$ contains an estimate of the standard error of the forecast of ${z}_{\mathit{i},n+\mathit{l}}$, for $\mathit{i}=1,2,\dots ,k$ and $\mathit{l}=1,2,\dots ,{\mathit{l}}_{\mathrm{max}}$.
18: $\mathbf{ref}\left({\mathbf{lref}}\right)$Real (Kind=nag_wp) array Output
On exit: the reference vector which may be used to update forecasts using g13dkf. The first $\left({\mathbf{lmax}}-1\right)×{\mathbf{k}}×{\mathbf{k}}$ elements contain the $\psi$ weight matrices, ${\psi }_{1},{\psi }_{2},\dots ,{\psi }_{{l}_{\mathrm{max}}-1}$. The next ${\mathbf{k}}×{\mathbf{lmax}}$ elements contain the forecasts of the transformed series ${\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.
19: $\mathbf{lref}$Integer Input
On entry: the dimension of the array ref as declared in the (sub)program from which g13djf is called.
Constraint: ${\mathbf{lref}}\ge \left({\mathbf{lmax}}-1\right)×{\mathbf{k}}×{\mathbf{k}}+2×{\mathbf{k}}×{\mathbf{lmax}}+{\mathbf{k}}$.
20: $\mathbf{work}\left({\mathbf{lwork}}\right)$Real (Kind=nag_wp) array Workspace
21: $\mathbf{lwork}$Integer Input
On entry: the dimension of the array work as declared in the (sub)program from which g13djf is called.
Constraint: if $r=\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{ip}},{\mathbf{iq}}\right)$ and $\mathit{d}=\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{id}}\left(\mathit{i}\right)\right)$, for $\mathit{i}=1,2,\dots ,k$, ${\mathbf{lwork}}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left\{{\mathbf{k}}r\left({\mathbf{k}}r+2\right),\left({\mathbf{ip}}+\mathit{d}+2\right){{\mathbf{k}}}^{2}+\left({\mathbf{n}}+{\mathbf{lmax}}\right){\mathbf{k}}\right\}$.
22: $\mathbf{iwork}\left({\mathbf{liwork}}\right)$Integer array Workspace
23: $\mathbf{liwork}$Integer Input
On entry: the dimension of the array iwork as declared in the (sub)program from which g13djf is called.
Constraint: ${\mathbf{liwork}}\ge {\mathbf{k}}×\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{ip}},{\mathbf{iq}}\right)$.
24: $\mathbf{ifail}$Integer Input/Output
On entry: ifail must be set to $0$, $-1$ or $1$ to set behaviour on detection of an error; these values have no effect when no error is detected.
A value of $0$ causes the printing of an error message and program execution will be halted; otherwise program execution continues. A value of $-1$ means that an error message is printed while a value of $1$ means that it is not.
If halting is not appropriate, the value $-1$ or $1$ is recommended. If message printing is undesirable, then the value $1$ is recommended. Otherwise, the value $0$ is recommended. When the value $-\mathbf{1}$ 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).

## 6Error Indicators and Warnings

If on entry ${\mathbf{ifail}}=0$ or $-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, $i=⟨\mathit{\text{value}}⟩$, ${\mathbf{id}}\left(i\right)=⟨\mathit{\text{value}}⟩$ and ${\mathbf{n}}-\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{ip}},{\mathbf{iq}}\right)=⟨\mathit{\text{value}}⟩$.
Constraint: $0\le {\mathbf{id}}\left(\mathit{i}\right)<{\mathbf{n}}-\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{ip}},{\mathbf{iq}}\right)$.
On entry, ${\mathbf{ip}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{ip}}\ge 0$.
On entry, ${\mathbf{iq}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{iq}}\ge 0$.
On entry, ${\mathbf{k}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{k}}\ge 1$.
On entry, ${\mathbf{kmax}}=⟨\mathit{\text{value}}⟩$ and ${\mathbf{k}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{kmax}}\ge {\mathbf{k}}$.
On entry, ${\mathbf{liwork}}=⟨\mathit{\text{value}}⟩$ and the minimum size $\text{required}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{liwork}}\ge {\mathbf{k}}×\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{ip}},{\mathbf{iq}}\right)$.
On entry, ${\mathbf{lmax}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{lmax}}\ge 1$.
On entry, lpar is too small: ${\mathbf{lpar}}=⟨\mathit{\text{value}}⟩$ and the minimum size $\text{required}=⟨\mathit{\text{value}}⟩$.
On entry, ${\mathbf{lref}}=⟨\mathit{\text{value}}⟩$ and the minimum size $\text{required}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{lref}}\ge \left({\mathbf{lmax}}-1\right)×{\mathbf{k}}×{\mathbf{k}}+2×{\mathbf{k}}×{\mathbf{lmax}}+{\mathbf{k}}$.
On entry, lwork is too small: ${\mathbf{lwork}}=⟨\mathit{\text{value}}⟩$ and the minimum size $\text{required}=⟨\mathit{\text{value}}⟩$.
On entry, mean is an invalid character.
Constraint: ${\mathbf{mean}}=\text{'M'}$ or $\text{'Z'}$.
On entry, ${\mathbf{n}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{n}}\ge 3$.
On entry, number of observations $\text{}=⟨\mathit{\text{value}}⟩$ and number of parameters $\text{}=⟨\mathit{\text{value}}⟩$.
Constraint: number of $\text{observations}\ge \text{number}$ of parameters.
${\mathbf{ifail}}=2$
On entry, $i=⟨\mathit{\text{value}}⟩$ and ${\mathbf{tr}}\left(i\right)$ is invalid.
Constraint: ${\mathbf{tr}}\left(i\right)=\text{'N'}$, $\text{'L'}$ or $\text{'S'}$.
${\mathbf{ifail}}=3$
On entry, one (or more) of the transformations requested is invalid. Check that you are not trying to log or square-root a series, some of whose values are negative.
${\mathbf{ifail}}=4$
On entry, the AR parameter matrices are outside the stationarity region. To proceed you must supply different parameter estimates in the arrays par and qq.
On entry, the covariance matrix qq is not positive definite. To proceed you must supply different parameter estimates in the arrays par and qq.
On entry, the MA parameter matrices are outside the invertibility region. To proceed you must supply different parameter estimates in the arrays par and qq.
${\mathbf{ifail}}=5$
An excessive number of iterations were needed to evaluate the eigenvalues of the matrices used to test for stationarity and invertibility.
${\mathbf{ifail}}=6$
The covariance matrix may be 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.
${\mathbf{ifail}}=7$
The forecasts will overflow if computed. You should check whether the transformations requested in the array tr are sensible.
${\mathbf{ifail}}=-99$
See Section 7 in the Introduction to the NAG Library FL Interface for further information.
${\mathbf{ifail}}=-399$
Your licence key may have expired or may not have been installed correctly.
See Section 8 in the Introduction to the NAG Library FL Interface for further information.
${\mathbf{ifail}}=-999$
Dynamic memory allocation failed.
See Section 9 in the Introduction to the NAG Library FL Interface for further information.

## 7Accuracy

The matrix computations are believed to be stable.

## 8Parallelism and Performance

g13djf is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
g13djf 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 routine. Please also consult the Users' Note for your implementation for any additional implementation-specific information.

The same differencing operator does not have to be applied to all the series. For example, suppose we have $k=2$, and wish to apply the second order differencing operator ${\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.$
Then ${\mathit{d}}_{1}=2,{\mathit{d}}_{2}=1$, $\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 ] .$
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 routine to produce the forecasts.
g13djf should not be used when the moving average parameters lie close to the boundary of the invertibility region. The routine does test for both invertibility and stationarity but if in doubt, you may use g13dxf, before calling this routine, to check that the VARMA model being used is invertible.
On a successful exit, the quantities k, lmax, kmax, ref and lref will be suitable for input to g13dkf.

## 10Example

This example computes forecasts of the next five values in two series each of length $48$. No transformation is to be used and no differencing is to be applied to either of the series. g13ddf is first called to fit an AR(1) model to the series. The mean vector $\mu$ is to be estimated and ${\varphi }_{1}\left(2,1\right)$ constrained to be zero.

### 10.1Program Text

Program Text (g13djfe.f90)

### 10.2Program Data

Program Data (g13djfe.d)

### 10.3Program Results

Program Results (g13djfe.r)