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 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 g13dkc.

2 Specification

#include <nag.h>

void

g13djc (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)

The function may be called by the names: g13djc, nag_tsa_multi_varma_forecast or nag_tsa_varma_forecast.

3 Description

Let the vector

Z_{t} = {(z_{1 t}, z_{2 t}, \dots, z_{k t})}^{T}

, for

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_{\max}}

are required. Let

W_{t} = {(w_{1 t}, w_{2 t}, \dots, w_{k t})}^{T}

be defined as follows:

w_{i t} = δ_{i} (B) z_{i t}^{*}, i = 1, 2, \dots, k,

where

δ_{i} (B)

is the differencing operator applied to the

i

th series and where

z_{i t}^{*}

is equal to either

z_{i t}

\sqrt{z_{i t}}

\log_{e} (z_{i t})

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

d_{i}

, then the differencing operator for the

i

th series is defined by

δ_{i} (B) = 1 - δ_{i 1} B - δ_{i 2} B^{2} - \dots - δ_{i d_{i}} B^{d_{i}}

where

B

is the backward shift operator; that is,

B Z_{t} = Z_{t - 1}

. The differencing parameters

δ_{i j}

, for

i = 1, 2, \dots, k

and

j = 1, 2, \dots, d_{i}

, must be supplied by you. If the

i

th series does not require differencing, then

d_{i} = 0

W_{t}

is assumed to follow a multivariate ARMA model of the form:

W_{t} - μ = ϕ_{1} (W_{t - 1} - μ) + ϕ_{2} (W_{t - 2} - μ) + \dots + ϕ_{p} (W_{t - p} - μ) + ε_{t} - θ_{1} ε_{t - 1} - \dots - θ_{q} ε_{t - q},

(1)

where

ε_{t} = {(ε_{1 t}, ε_{2 t}, \dots, ε_{k t})}^{T}

, for

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

Σ

. The components of

ε_{t}

are assumed to be uncorrelated at non-simultaneous lags. The

ϕ_{i}

and

θ_{j}

are

k

k

matrices of parameters. The matrices

ϕ_{i}

, for

i = 1, 2, \dots, p

, are the autoregressive (AR) parameter matrices, and the matrices

θ_{i}

, for

i = 1, 2, \dots, q

, the moving average (MA) parameter matrices. The parameters in the model are thus the

p

(

k

k

)

ϕ

-matrices, the

q

(

k

k

)

θ

-matrices, the mean vector

μ

and the residual error covariance matrix

Σ

. The ARMA model (1) must be both stationary and invertible; see g13dxc for a method of checking these conditions.

The ARMA model (1) may be rewritten as

ϕ (B) (δ (B) Z_{t}^{*} - μ) = θ (B) ε_{t},

where

ϕ (B)

and

θ (B)

are the autoregressive and moving average polynomials and

δ (B)

denotes the

k

k

diagonal matrix whose

i

th diagonal elements is

δ_{i} (B)

and

Z_{t}^{*} = {(z_{1 t}^{*}, z_{2 t}^{*} \dots z_{k t}^{*})}^{T}

This may be rewritten as

ϕ (B) δ (B) Z_{t}^{*} = ϕ (B) μ + θ (B) ε_{t}

Z_{t}^{*} = τ + ψ (B) ε_{t} = τ + ε_{t} + ψ_{1} ε_{t - 1} + ψ_{2} ε_{t - 2} + \dots

where

ψ (B) = δ^{- 1} (B) ϕ^{- 1} (B) θ (B)

and

τ = δ^{- 1} (B) μ

is a vector of length

k

Forecasts are computed using a multivariate version of the procedure described in Box and Jenkins (1976). If

{\hat{Z}}_{n}^{*} (l)

denotes the forecast of

Z_{n + l}^{*}

, then

{\hat{Z}}_{n}^{*} (l)

is taken to be that linear function of

Z_{n}^{*}, Z_{n - 1}^{*}, \dots

which minimizes the elements of

E \{e_{n} (l) e_{n}^{'} (l)\}

where

e_{n} (l) = Z_{n + l}^{*} - {\hat{Z}}_{n}^{*} (l)

is the forecast error.

{\hat{Z}}_{n}^{*} (l)

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

{\hat{Z}}_{n}^{*} (l) = τ + ψ_{l} ε_{n} + ψ_{l + 1} ε_{n - 1} + ψ_{l + 2} ε_{n - 2} + \dots .

The forecast error at

t

for lead

l

is then

e_{n} (l) = Z_{n + l}^{*} - {\hat{Z}}_{n}^{*} (l) = ε_{n + l} + ψ_{1} ε_{n + l - 1} + ψ_{2} ε_{n + l - 2} + \dots + ψ_{l - 1} ε_{n + 1} .

Let

d = \max (d_{i})

, for

i = 1, 2, \dots, k

. Unless

q = 0

the function requires estimates of

ε_{t}

, for

t = d + 1, \dots, n

, which are obtainable from g13ddc. The terms

ε_{t}

are assumed to be zero, for

t = n + 1, \dots, n + l_{\max}

. You may use g13dkc to update these

l_{\max}

forecasts should further observations,

Z_{n + 1}, Z_{n + 2}, \dots

, become available. Note that when

l_{\max}

or more further observations are available then g13djc must be used to produce new forecasts for

Z_{n + l_{\max} + 1}, Z_{n + l_{\max} + 2}, \dots

, 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 quantities k, n, kmax, ip, iq, par, npar, qq and v from g13ddc are suitable for input to g13djc.

1: $k$ – Integer Input

On entry:

k

, the dimension of the multivariate time series.

Constraint:

k \geq 1

2: $n$ – Integer Input

On entry:

n

, the number of observations in the series,

Z_{t}

, prior to differencing.

Constraint:

n \geq 3

The total number of observations must exceed the total number of parameters in the model; that is

if $mean = Nag_MeanZero$ , $n \times k > (ip + iq) \times k \times k + k \times (k + 1) / 2$ ;
if $mean = Nag_MeanInclude$ , $n \times k > (ip + iq) \times k \times k + k + k \times (k + 1) / 2$ ,

(see the arguments ip, iq and mean).

3: $z [kmax \times n]$ – const double Input

On entry:

z [(t - 1) \times kmax + i - 1]

must contain the

i

th series at time

t

, for

t = 1, 2, \dots, n

and

i = 1, 2, \dots, k

4: $kmax$ – Integer Input

On entry: the stride separating row elements in the two-dimensional data stored in the arrays z, delta, qq, v, predz, sefz.

Constraint:

kmax \geq k

5: $tr [k]$ – const Integer Input

On entry:

tr [i - 1]

indicates whether the

i

th series is to be transformed, for

i = 1, 2, \dots, k

$tr [i - 1] = - 1$: A square root transformation is used.
$tr [i - 1] = 0$: No transformation is used.
$tr [i - 1] = 1$: A log transformation is used.

Constraint:

tr [i - 1] = - 1

0

1

, for

i = 1, 2, \dots, k

6: $id [k]$ – const Integer Input

On entry:

id [i - 1]

must specify,

d_{i}

, the order of differencing required for the

i

th series.

Constraint:

0 \leq id [i - 1] < n - \max (ip, iq)

, for

i = 1, 2, \dots, k

7: $delta [\dim]$ – const double Input

Note: the dimension, dim, of the array delta must be at least

kmax \times d

, where

d = \max (id [i - 1])

On entry: if

id [i - 1] > 0

, then

delta [(j - 1) \times kmax + i - 1]

must be set equal to

δ_{i j}

, for

j = 1, 2, \dots, d_{i}

and

i = 1, 2, \dots, k

d = 0

, delta is not referenced.

8: $ip$ – Integer Input

On entry:

p

, the number of AR parameter matrices.

Constraint:

ip \geq 0

9: $iq$ – Integer Input

On entry:

q

, the number of MA parameter matrices.

Constraint:

iq \geq 0

10: $mean$ – Nag_IncludeMean Input

On entry:

mean = Nag_MeanInclude

, if components of

μ

have been estimated and

mean = Nag_MeanZero

, if all elements of

μ

are to be taken as zero.

Constraint:

mean = Nag_MeanInclude

Nag_MeanZero

11: $par [lpar]$ – const double Input

On entry: must contain the parameter estimates read in row by row in the order

ϕ_{1}, ϕ_{2}, \dots, ϕ_{p}

θ_{1}, θ_{2}, \dots, θ_{q}

μ

Thus,

if $ip > 0$ , $par [(l - 1) \times k \times k + (i - 1) \times k + j - 1]$ must be set equal to an estimate of the $(i, j)$ th element of $ϕ_{l}$ , for $l = 1, 2, \dots, p$ , $i = 1, 2, \dots, k$ and $j = 1, 2, \dots, k$ ;
if $iq > 0$ , $par [p \times k \times k + (l - 1) \times k \times k + (i - 1) \times k + j - 1]$ must be set equal to an estimate of the $(i, j)$ th element of $θ_{l}$ , for $l = 1, 2, \dots, q$ , $i = 1, 2, \dots, k$ and $j = 1, 2, \dots, k$ ;
if $mean = Nag_MeanInclude$ , $par [(p + q) \times k \times k + i - 1]$ must be set equal to an estimate of the $i$ th component of $μ$ , for $i = 1, 2, \dots, k$ .

Constraint: the first

ip \times k \times k

elements of par must satisfy the stationarity condition and the next

iq \times k \times k

elements of par must satisfy the invertibility condition.

12: $lpar$ – Integer Input

On entry: the dimension of the array par.

Constraints:

if $mean = Nag_MeanZero$ , $lpar \geq \max (1, (ip + iq) \times k \times k)$ ;
if $mean = Nag_MeanInclude$ , $lpar \geq (ip + iq) \times k \times k + k$ .

13: $qq [kmax \times k]$ – double Input/Output

On entry:

qq [(j - 1) \times kmax + i - 1]

must contain an estimate of the

(i, j)

th element of

Σ

. The lower triangle only is needed.

Constraint:

qq

must be positive definite.

On exit: if

fail . code =

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 [\dim]$ – const double Input

Note: the dimension, dim, of the array v must be at least

kmax \times (n - d)

, where

d = \max (id [i - 1])

On entry:

v [(t - 1) \times kmax + i - 1]

must contain an estimate of the

i

th component of

ε_{t + d}

, for

i = 1, 2, \dots, k

and

t = 1, 2, \dots, n - d

q = 0

, v is not used.

15: $lmax$ – Integer Input

On entry: the number,

l_{\max}

, of forecasts required.

Constraint:

lmax \geq 1

16: $predz [kmax \times lmax]$ – double Output

On exit:

predz [(l - 1) \times kmax + i - 1]

contains the forecast of

z_{i, n + l}

, for

i = 1, 2, \dots, k

and

l = 1, 2, \dots, l_{\max}

17: $sefz [kmax \times lmax]$ – double Output

On exit:

sefz [(l - 1) \times kmax + i - 1]

contains an estimate of the standard error of the forecast of

z_{i, n + l}

, for

i = 1, 2, \dots, k

and

l = 1, 2, \dots, l_{\max}

18: $ref [lref]$ – double Output

On exit: the reference vector which may be used to update forecasts using g13dkc. The first

(lmax - 1) \times k \times k

elements contain the

ψ

weight matrices,

ψ_{1}, ψ_{2}, \dots, ψ_{l_{\max} - 1}

. The next

k \times lmax

elements contain the forecasts of the transformed series

{\hat{Z}}_{n + 1}^{*}, {\hat{Z}}_{n + 2}^{*}, \dots, {\hat{Z}}_{n + l_{\max}}^{*}

and the next

k \times 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: $lref$ – Integer Input

On entry: the dimension of the array ref.

Constraint:

lref \geq (lmax - 1) \times k \times k + 2 \times k \times lmax + k

20: $fail$ – NagError * Input/Output

The NAG error argument (see Section 7 in the Introduction to the NAG Library CL Interface).

6 Error Indicators and Warnings

NE_ALLOC_FAIL: Dynamic memory allocation failed.
See Section 3.1.2 in the Introduction to the NAG Library CL Interface for further information.
NE_BAD_PARAM: On entry, argument $〈value〉$ had an illegal value.
NE_EIGENVALUES: An excessive number of iterations were needed 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. To proceed you must supply different parameter estimates in the arrays par and qq.
NE_G13D_MA: 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.
NE_INT: On entry, $ip = 〈value〉$ .
Constraint: $ip \geq 0$ .

On entry, $iq = 〈value〉$ .
Constraint: $iq \geq 0$ .

On entry, $k = 〈value〉$ .
Constraint: $k \geq 1$ .

On entry, $lmax = 〈value〉$ .
Constraint: $lmax \geq 1$ .

On entry, lpar is too small: $lpar = 〈value〉$ and the minimum size $required = 〈value〉$ .

On entry, $lref = 〈value〉$ and the minimum size $required = 〈value〉$ .
Constraint: $lref \geq (lmax - 1) \times k \times k + 2 \times k \times lmax + k$ .

On entry, $n = 〈value〉$ .
Constraint: $n \geq 3$ .
NE_INT_2: On entry, $kmax = 〈value〉$ and $k = 〈value〉$ .
Constraint: $kmax \geq k$ .
NE_INT_ARRAY: On entry, $i = 〈value〉$ , $id [i - 1] = 〈value〉$ and $n - \max (ip, iq) = 〈value〉$ .
Constraint: $0 \leq id [i - 1] < n - \max (ip, iq)$ .

On entry, $tr [〈value〉] = 〈value〉$ .
Constraint: $tr [i] = - 1$ , $0$ or $1$ .
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.
See Section 7.5 in the Introduction to the NAG Library CL Interface for further information.
NE_NEARLY_POS_DEF: 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.
NE_NO_LICENCE: Your licence key may have expired or may not have been installed correctly.
See Section 8 in the Introduction to the NAG Library CL Interface for further information.
NE_NOT_POS_DEF: On entry, the covariance matrix qq is not positive definite. To proceed you must supply different parameter estimates in the arrays par and qq.
NE_OBSERV_LT_P: On entry, number of observations $= 〈value〉$ and number of parameters $= 〈value〉$ .
Constraint: number of $observations \geq number$ of parameters.
NE_OVERFLOW_LIKELY: The forecasts will overflow if computed. You should check whether the transformations requested in the array tr are sensible.
NE_TRANSFORMATION: 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.

7 Accuracy

The matrix computations are believed to be stable.

8 Parallelism and Performance

g13djc is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.

g13djc 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 function. Please also consult the Users' Note for your implementation for any additional implementation-specific information.

9 Further Comments

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:

\begin{matrix} w_{1 t} = \nabla^{2} z_{1 t} = {(1 - B)}^{2} z_{1 t} = (1 - 2 B + B^{2}) Z_{1 t},   and \\ w_{2 t} = \nabla z_{2 t} = (1 - B) z_{2 t} . \end{matrix}

Then

d_{1} = 2, d_{2} = 1

d = \max (d_{1}, d_{2}) = 2

, and

delta = [\begin{array}{r} δ_{11} & δ_{12} \\ δ_{21} \end{array}] = [\begin{array}{r} 2 & - 1 \\ 1 \end{array}] .

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.

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 g13dxc, before calling this function, 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 g13dkc.

10 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. g13ddc is first called to fit an AR(1) model to the series. The mean vector

μ

is to be estimated and

ϕ_{1} (2, 1)

constrained to be zero.

g13dj: FL CL CPP AD

NAG CL Interfaceg13djc (multi_​varma_​forecast)

▸▿ Contents

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

6 Error Indicators and Warnings

7 Accuracy

8 Parallelism and Performance

9 Further Comments

10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results

NAG CL Interface
g13djc (multi_varma_forecast)