g13 Chapter Contents
g13 Chapter Introduction
NAG C Library Manual

# NAG Library Function Documentnag_tsa_varma_forecast (g13djc)

## 1  Purpose

nag_tsa_varma_forecast (g13djc) 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_varma_estimate (g13ddc). 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_varma_update (g13dkc).

## 2  Specification

 #include #include
 void nag_tsa_varma_forecast (Integer k, Integer n, const double z[], Integer kmax, const Integer tr[], const Integer id[], const double delta[], Integer ip, Integer iq, Nag_IncludeMean mean, const double par[], Integer lpar, double qq[], const double v[], Integer lmax, double predz[], double sefz[], double ref[], Integer lref, NagError *fail)

## 3  Description

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=δiBzit*, 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}}_{e}{z}_{it}$ 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-μ=ϕ1Wt-1-μ+ϕ2Wt-2-μ+⋯+ϕpWt-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$ by $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$ by $k$) $\varphi$-matrices, the $q$ ($k$ by $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_arma_roots (g13dxc) for a method of checking these conditions.
The ARMA model (1) may be rewritten as
 $ϕBδBZt*-μ=θ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$ by $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δBZt*=ϕ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
 $enl=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 function requires estimates of ${\epsilon }_{\mathit{t}}$, for $\mathit{t}=\mathit{d}+1,\dots ,n$, which are obtainable from nag_tsa_varma_estimate (g13ddc). The terms ${\epsilon }_{\mathit{t}}$ are assumed to be zero, for $\mathit{t}=n+1,\dots ,n+{l}_{\mathrm{max}}$. You may use nag_tsa_varma_update (g13dkc) 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 nag_tsa_varma_forecast (g13djc) 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).

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

## 5  Arguments

The output quantities k, n, kmax, ip, iq, par, npar, qq and v from nag_tsa_varma_estimate (g13ddc) are suitable for input to nag_tsa_varma_forecast (g13djc).
1:     kIntegerInput
On entry: $k$, the dimension of the multivariate time series.
Constraint: ${\mathbf{k}}\ge 1$.
2:     nIntegerInput
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}}=\mathrm{Nag_MeanZero}$, ${\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}}=\mathrm{Nag_MeanInclude}$, ${\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:     z[${\mathbf{kmax}}×{\mathbf{n}}$]const doubleInput
On entry: ${\mathbf{z}}\left[\left(\mathit{t}-1\right){\mathbf{kmax}}+\mathit{i}-1\right]$ must contain the $\mathit{i}$th series at time $\mathit{t}$, for $\mathit{t}=1,2,\dots ,n$ and $\mathit{i}=1,2,\dots ,k$.
4:     kmaxIntegerInput
On entry: the first dimension of the arrays z, qq, predz and sefz.
Constraint: ${\mathbf{kmax}}\ge {\mathbf{k}}$.
5:     tr[k]const IntegerInput
On entry: ${\mathbf{tr}}\left[\mathit{i}-1\right]$ indicates whether the $\mathit{i}$th series is to be transformed, for $\mathit{i}=1,2,\dots ,k$.
${\mathbf{tr}}\left[i-1\right]=-1$
A square root transformation is used.
${\mathbf{tr}}\left[i-1\right]=0$
No transformation is used.
${\mathbf{tr}}\left[i-1\right]=1$
A log transformation is used.
Constraint: ${\mathbf{tr}}\left[\mathit{i}-1\right]=-1$, $0$ or $1$, for $\mathit{i}=1,2,\dots ,k$.
6:     id[k]const IntegerInput
On entry: ${\mathbf{id}}\left[i-1\right]$ must specify, ${\mathit{d}}_{i}$, the order of differencing required for the $i$th series.
Constraint: $0\le {\mathbf{id}}\left[\mathit{i}-1\right]<{\mathbf{n}}-\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{ip}},{\mathbf{iq}}\right)$, for $\mathit{i}=1,2,\dots ,k$.
7:     delta[$\mathit{dim}$]const doubleInput
Note: the dimension, dim, of the array delta must be at least $\mathit{d}×{\mathbf{kmax}}$, where $\mathit{d}=\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{id}}\left[i-1\right]\right)$.
On entry: if ${\mathbf{id}}\left[i-1\right]>0$, then ${\mathbf{delta}}\left[×\left(\mathit{j}-1\right)+\mathit{i}-1\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:     ipIntegerInput
On entry: $p$, the number of AR parameter matrices.
Constraint: ${\mathbf{ip}}\ge 0$.
9:     iqIntegerInput
On entry: $q$, the number of MA parameter matrices.
Constraint: ${\mathbf{iq}}\ge 0$.
10:   meanNag_IncludeMeanInput
On entry: ${\mathbf{mean}}=\mathrm{Nag_MeanInclude}$, if components of $\mu$ have been estimated and ${\mathbf{mean}}=\mathrm{Nag_MeanZero}$, if all elements of $\mu$ are to be taken as zero.
Constraint: ${\mathbf{mean}}=\mathrm{Nag_MeanInclude}$ or $\mathrm{Nag_MeanZero}$.
11:   par[lpar]const doubleInput
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}-1\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}-1\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}}=\mathrm{Nag_MeanInclude}$, ${\mathbf{par}}\left[\left(p+q\right)×k×k+\mathit{i}-1\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:   lparIntegerInput
On entry: the dimension of the array par.
Constraints:
• if ${\mathbf{mean}}=\mathrm{Nag_MeanZero}$, ${\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}}=\mathrm{Nag_MeanInclude}$, ${\mathbf{lpar}}\ge \left({\mathbf{ip}}+{\mathbf{iq}}\right)×{\mathbf{k}}×{\mathbf{k}}+{\mathbf{k}}$.
13:   qq[${\mathbf{kmax}}×{\mathbf{k}}$]doubleInput/Output
On entry: ${\mathbf{qq}}\left[{\mathbf{kmax}}×\left(j-1\right)+i-1\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 NE_EIGENVALUES, NE_G13D_AR, NE_G13D_MA, NE_NEARLY_POS_DEF, NE_NOT_POS_DEF, NE_OVERFLOW_LIKELY or NE_TRANSFORMATION, then the upper triangle is set equal to the lower triangle.
14:   v[$\mathit{dim}$]const doubleInput
Note: the dimension, dim, of the array v must be at least ${\mathbf{kmax}}×\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-1\right]\right)$.
On entry: ${\mathbf{v}}\left[×\left(\mathit{t}-1\right)+\mathit{i}-1\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:   lmaxIntegerInput
On entry: the number, ${l}_{\mathrm{max}}$, of forecasts required.
Constraint: ${\mathbf{lmax}}\ge 1$.
16:   predz[${\mathbf{kmax}}×{\mathbf{lmax}}$]doubleOutput
On exit: ${\mathbf{predz}}\left[{\mathbf{kmax}}×\left(\mathit{l}-1\right)+\mathit{i}-1\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:   sefz[${\mathbf{kmax}}×{\mathbf{lmax}}$]doubleOutput
On exit: ${\mathbf{sefz}}\left[{\mathbf{kmax}}×\left(\mathit{l}-1\right)+\mathit{i}-1\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:   ref[lref]doubleOutput
On exit: the reference vector which may be used to update forecasts using nag_tsa_varma_update (g13dkc). 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:   lrefIntegerInput
On entry: the dimension of the array ref.
Constraint: ${\mathbf{lref}}\ge \left({\mathbf{lmax}}-1\right)×{\mathbf{k}}×{\mathbf{k}}+2×{\mathbf{k}}×{\mathbf{lmax}}+{\mathbf{k}}$.
20:   failNagError *Input/Output
The NAG error argument (see Section 3.6 in the Essential Introduction).

## 6  Error Indicators and Warnings

NE_ALLOC_FAIL
Dynamic memory allocation failed.
On entry, argument $〈\mathit{\text{value}}〉$ had an illegal value.
NE_EIGENVALUES
An excessive number of iterations were needed by nag_real_eigenvalues (f02afc) to evaluate the eigenvalues of the matrices used to test for stationarity and invertibility.
NE_G13D_AR
On entry, the AR parameter matrices are outside the stationarity region.
NE_G13D_MA
On entry, the MA parameter matrices are outside the invertibility region.
NE_INT
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{lmax}}=〈\mathit{\text{value}}〉$.
Constraint: ${\mathbf{lmax}}\ge 1$.
On entry, lpar is too small: ${\mathbf{lpar}}=〈\mathit{\text{value}}〉$ but must be at least $〈\mathit{\text{value}}〉$.
On entry, lref is too small: ${\mathbf{lref}}=〈\mathit{\text{value}}〉$ but must be at least $〈\mathit{\text{value}}〉$.
On entry, ${\mathbf{n}}=〈\mathit{\text{value}}〉$.
Constraint: ${\mathbf{n}}\ge 3$.
NE_INT_2
On entry, ${\mathbf{kmax}}=〈\mathit{\text{value}}〉$ and ${\mathbf{k}}=〈\mathit{\text{value}}〉$.
Constraint: ${\mathbf{kmax}}\ge {\mathbf{k}}$.
NE_INT_ARRAY
On entry, ${\mathbf{id}}\left[〈\mathit{\text{value}}〉\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}-1\right]<{\mathbf{n}}-\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{ip}},{\mathbf{iq}}\right)$.
On entry, ${\mathbf{tr}}\left[\mathit{i}-1\right]=〈\mathit{\text{value}}〉$ and ${\mathbf{k}}=〈\mathit{\text{value}}〉$.
Constraint: ${\mathbf{tr}}\left[\mathit{i}-1\right]=-1$, $0$ or $1$, for $\mathit{i}=1,2,\dots ,k$.
NE_INTERNAL_ERROR
An internal error has occurred in this function. Check the function call and any array sizes. If the call is correct then please contact NAG for assistance.
NE_NEARLY_POS_DEF
The covariance matrix may be nearly non positive definite.
NE_NOT_POS_DEF
On entry, the covariance matrix qq is not positive definite.
NE_OBSERV_LT_P
On entry, the total number of observations is less than the total number of parameters (including the covariance matrix). Number of observations $\text{}=〈\mathit{\text{value}}〉$ and number of parameters $\text{}=〈\mathit{\text{value}}〉$.
NE_OVERFLOW_LIKELY
The forecasts will overflow if computed.
NE_TRANSFORMATION
On entry, one (or more) of the transformations requested is invalid.

## 7  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$, 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+B2Z1t, and w2t=∇z2t=1-Bz2t.$
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 function to produce the forecasts.
nag_tsa_varma_forecast (g13djc) 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_arma_roots (g13dxc), 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_varma_update (g13dkc).

## 9  Example

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. nag_tsa_varma_estimate (g13ddc) 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.

### 9.1  Program Text

Program Text (g13djce.c)

### 9.2  Program Data

Program Data (g13djce.d)

### 9.3  Program Results

Program Results (g13djce.r)