NAG Library Routine Document

g13ddf (multi_varma_estimate)

1

Purpose

2

Specification

3

Description

4

References

5

Arguments

1:     $\mathbf{k}$ – IntegerInput

On entry: $k$, the number of observed time series.

Constraint: $\mathbf{k}\ge 1$.

2:     $\mathbf{n}$ – IntegerInput

On entry: $n$, the number of observations in each time series.

3:     $\mathbf{ip}$ – IntegerInput

On entry: $p$, the number of AR parameter matrices.

Constraint: $\mathbf{ip}\ge 0$.

4:     $\mathbf{iq}$ – IntegerInput

On entry: $q$, the number of MA parameter matrices.

Constraint: $\mathbf{iq}\ge 0$.

$\mathbf{ip}=\mathbf{iq}=0$ is not permitted.

5:     $\mathbf{mean}$ – LogicalInput

On entry: $\mathbf{mean}=\mathrm{.TRUE.}$, if components of $\mu$ have been estimated and $\mathbf{mean}=\mathrm{.FALSE.}$, if all elements of $\mu$ are to be taken as zero.

Constraint: $\mathbf{mean}=\mathrm{.TRUE.}$ or $\mathrm{.FALSE.}$.

6:     $\mathbf{par}\left(\mathbf{npar}\right)$ – Real (Kind=nag_wp) arrayInput/Output

On entry: initial 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)\times k\times k+\left(\mathit{i}-1\right)\times k+\mathit{j}\right)$ must be set equal to an initial 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\times k\times k+\left(l-1\right)\times k\times k+\left(i-1\right)\times k+j\right)$ must be set equal to an initial estimate of the $\left(i,j\right)$th element of ${\theta }_l$, $l=1,2,\dots ,q$ and $i,j=1,2,\dots ,k$;
• if $\mathbf{mean}=\mathrm{.TRUE.}$, $\mathbf{par}\left(\left(p+q\right)\times k\times k+i\right)$ should be set equal to an initial estimate of the $i$th component of $\mu$ ($\mu \left(i\right)$). (If you set $\mathbf{par}\left(\left(p+q\right)\times k\times k+i\right)$ to $0.0$ then g13ddf will calculate the mean of the $i$th series and use this as an initial estimate of $\mu \left(i\right)$.)

The first $p\times k\times k$ elements of par must satisfy the stationarity condition and the next $q\times k\times k$ elements of par must satisfy the invertibility condition.

If in doubt set all elements of par to $0.0$.

On exit: if $\mathbf{ifail}=\mathbf{0}$ or $\mathbf{ifail}\ge \mathbf{4}$ then all the elements of par will be overwritten by the latest estimates of the corresponding ARMA parameters.

7:     $\mathbf{npar}$ – IntegerInput

On entry: the dimension of the arrays par, parhld and g and the second dimension of the array cm as declared in the (sub)program from which g13ddf is called. npar is the number of initial parameter estimates.

Constraints:

• if $\mathbf{mean}=\mathrm{.FALSE.}$, npar must be set equal to $\left(p+q\right)\times k\times k$;
• if $\mathbf{mean}=\mathrm{.TRUE.}$, npar must be set equal to $\left(p+q\right)\times k\times k+k$.

The total number of observations $\left(n\times k\right)$ must exceed the total number of parameters in the model ($\mathbf{npar}+k\left(k+1\right)/2$).

8:     $\mathbf{qq}\left(\mathbf{kmax},\mathbf{k}\right)$ – Real (Kind=nag_wp) arrayInput/Output

On entry: $\mathbf{qq}\left(\mathit{i},\mathit{j}\right)$ must be set equal to an initial estimate of the $\left(\mathit{i},\mathit{j}\right)$th element of $\Sigma$. The lower triangle only is needed. qq must be positive definite. It is strongly recommended that on entry the elements of qq are of the same order of magnitude as at the solution point. If you set $\mathbf{qq}\left(\mathit{i},\mathit{j}\right)=0.0$, for $\mathit{i}=1,2,\dots ,k$ and $\mathit{j}=1,2,\dots ,i$, then g13ddf will calculate the covariance matrix between the $k$ time series and use this as an initial estimate of $\Sigma$.

On exit: if $\mathbf{ifail}=\mathbf{0}$ or $\mathbf{ifail}\ge \mathbf{4}$ then $\mathbf{qq}\left(i,j\right)$ will contain the latest estimate of the $\left(i,j\right)$th element of $\Sigma$. The lower triangle only is returned.

9:     $\mathbf{kmax}$ – IntegerInput

On entry: the first dimension of the arrays qq, w and v as declared in the (sub)program from which g13ddf is called.

Constraint: $\mathbf{kmax}\ge \mathbf{k}$.

10:   $\mathbf{w}\left(\mathbf{kmax},\mathbf{n}\right)$ – Real (Kind=nag_wp) arrayInput

On entry: $\mathbf{w}\left(\mathit{i},\mathit{t}\right)$ must be set equal to the $\mathit{i}$th component of $W_{\mathit{t}}$, for $\mathit{i}=1,2,\dots ,k$ and $\mathit{t}=1,2,\dots ,n$.

11:   $\mathbf{parhld}\left(\mathbf{npar}\right)$ – Logical arrayInput

On entry: $\mathbf{parhld}\left(\mathit{i}\right)$ must be set to .TRUE. if $\mathbf{par}\left(\mathit{i}\right)$ is to be held constant at its input value and .FALSE. if $\mathbf{par}\left(\mathit{i}\right)$ is a free parameter, for $\mathit{i}=1,2,\dots ,\mathbf{npar}$.

If in doubt try setting all elements of parhld to .FALSE..

12:   $\mathbf{exact}$ – LogicalInput

On entry: must be set equal to .TRUE. if you wish g13ddf to compute exact maximum likelihood estimates. exact must be set equal to .FALSE. if only conditional likelihood estimates are required.

13:   $\mathbf{iprint}$ – IntegerInput

On entry: the frequency with which the automatic monitoring routine is to be called.

$\mathbf{iprint}>0$

The ML search procedure is monitored once every iprint iterations and just before exit from the search routine.

$\mathbf{iprint}=0$

The search routine is monitored once at the final point.

$\mathbf{iprint}<0$

The search routine is not monitored at all.

14:   $\mathbf{cgetol}$ – Real (Kind=nag_wp)Input

On entry: the accuracy to which the solution in par and qq is required.

If cgetol is set to ${10}^{-l}$ and on exit $\mathbf{ifail}=\mathbf{0}$ or $\mathbf{ifail}\ge \mathbf{6}$, then all the elements in par and qq should be accurate to approximately $l$ decimal places. For most practical purposes the value ${10}^{-4}$ should suffice. You should be wary of setting cgetol too small since the convergence criteria may then have become too strict for the machine to handle.

If cgetol has been set to a value which is less than the machine precision, $\epsilon$, then g13ddf will use the value $10.0\times \sqrt{\epsilon }$ instead.

15:   $\mathbf{maxcal}$ – IntegerInput

On entry: the maximum number of likelihood evaluations to be permitted by the search procedure.

Suggested value: $\mathbf{maxcal}=40\times \mathbf{npar}\times \left(\mathbf{npar}+5\right)$.

Constraint: $\mathbf{maxcal}\ge 1$.

16:   $\mathbf{ishow}$ – IntegerInput

On entry: specifies which of the following two quantities are to be printed.

(i)	table of maximum likelihood estimates and their standard errors (as returned in the output arrays par, qq and cm);
(ii)	table of residual series (as returned in the output array v).

$\mathbf{ishow}=0$

None of the above are printed.

$\mathbf{ishow}=1$

(i) only is printed.

$\mathbf{ishow}=2$

(i) and (ii) are printed.

Constraint: $0\le \mathbf{ishow}\le 2$.

17:   $\mathbf{niter}$ – IntegerOutput

On exit: if $\mathbf{ifail}=\mathbf{0}$ or $\mathbf{ifail}\ge \mathbf{4}$ then niter contains the number of iterations performed by the search routine.

18:   $\mathbf{rlogl}$ – Real (Kind=nag_wp)Output

On exit: if $\mathbf{ifail}=\mathbf{0}$ or $\mathbf{ifail}\ge \mathbf{4}$ then rlogl contains the value of the log-likelihood function corresponding to the final point held in par and qq.

19:   $\mathbf{v}\left(\mathbf{kmax},\mathbf{n}\right)$ – Real (Kind=nag_wp) arrayOutput

On exit: if $\mathbf{ifail}=\mathbf{0}$ or $\mathbf{ifail}\ge \mathbf{4}$ then $\mathbf{v}\left(\mathit{i},\mathit{t}\right)$ will contain an estimate of the $\mathit{i}$th component of ${\epsilon }_{\mathit{t}}$, for $\mathit{i}=1,2,\dots ,k$ and $\mathit{t}=1,2,\dots ,n$, corresponding to the final point held in par and qq.

20:   $\mathbf{g}\left(\mathbf{npar}\right)$ – Real (Kind=nag_wp) arrayOutput

On exit: if $\mathbf{ifail}=\mathbf{0}$ or $\mathbf{ifail}\ge \mathbf{4}$ then $\mathbf{g}\left(i\right)$ will contain the estimated first derivative of the log-likelihood function with respect to the $i$th element in the array par. If the gradient cannot be computed then all the elements of g are returned as zero.

21:   $\mathbf{cm}\left(\mathbf{ldcm},\mathbf{npar}\right)$ – Real (Kind=nag_wp) arrayOutput

On exit: if $\mathbf{ifail}=\mathbf{0}$ or $\mathbf{ifail}\ge \mathbf{4}$ then $\mathbf{cm}\left(i,j\right)$ will contain an estimate of the correlation coefficient between the $i$th and $j$th elements in the par array for $1\le i\le \mathbf{npar}$, $1\le j\le \mathbf{npar}$. If $i=j$, then $\mathbf{cm}\left(i,j\right)$ will contain the estimated standard error of $\mathbf{par}\left(i\right)$. If the $l$th component of par has been held constant, i.e., $\mathbf{parhld}\left(l\right)$ was set to .TRUE., then the $l$th row and column of cm will be set to zero. If the second derivative matrix cannot be computed then all the elements of cm are returned as zero.

22:   $\mathbf{ldcm}$ – IntegerInput

On entry: the first dimension of the array cm as declared in the (sub)program from which g13ddf is called.

Constraint: $\mathbf{ldcm}\ge \mathbf{npar}$.

23:   $\mathbf{ifail}$ – IntegerInput/Output

On entry: ifail must be set to $0$, $-1\text{ or }1$. If you are unfamiliar with this argument you should refer to Section 3.4 in How to Use the NAG Library and its Documentation for details.

For environments where it might be inappropriate to halt program execution when an error is detected, the value $-1\text{ or }1$ is recommended. If the output of error messages is undesirable, then the value $1$ is recommended. Otherwise, because for this routine the values of the output arguments may be useful even if $\mathbf{ifail}\ne \mathbf{0}$ on exit, the recommended value is $-1$. When the value $-\mathbf{1}\text{ or }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).

6

Error Indicators and Warnings

$\mathbf{ifail}=1$

On entry, $\mathbf{ip}=〈\mathit{\text{value}}〉$.
Constraint: $\mathbf{ip}\ge 0$.

On entry, $\mathbf{ip}=0$ and $\mathbf{iq}=0$.

On entry, $\mathbf{iq}=〈\mathit{\text{value}}〉$.
Constraint: $\mathbf{iq}\ge 0$.

On entry, $\mathbf{ishow}=〈\mathit{\text{value}}〉$.
Constraint: $0\le \mathbf{ishow}\le 2$.

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{ldcm}=〈\mathit{\text{value}}〉$ and $\mathbf{npar}=〈\mathit{\text{value}}〉$.
Constraint: $\mathbf{ldcm}\ge \mathbf{npar}$.

On entry, $\mathbf{maxcal}=〈\mathit{\text{value}}〉$.
Constraint: $\mathbf{maxcal}\ge 1$.

On entry, $\mathbf{n}=〈\mathit{\text{value}}〉$, $\mathbf{k}=〈\mathit{\text{value}}〉$ and $\mathbf{npar}=〈\mathit{\text{value}}〉$.
Constraint: $\mathbf{n}\times \mathbf{k}>\mathbf{npar}+\mathbf{k}\times \left(\mathbf{k}+1\right)/2$.

On entry, $\mathbf{npar}=〈\mathit{\text{value}}〉$.
Constraint: $\mathbf{npar}=〈\mathit{\text{value}}〉$.

On entry, $\mathbf{npar}=〈\mathit{\text{value}}〉$.
Constraint: $\mathbf{npar}\ge 0$.

$\mathbf{ifail}=2$

The initial AR parameter estimates are outside the stationarity region. To proceed you must try a different starting point.

The initial estimate of $\Sigma$ is not positive definite. To proceed you must try a different starting point.

The initial MA parameter estimates are outside the invertibility region. To proceed you must try a different starting point.

The starting point is too close to the boundary of the admissibility region. To proceed you must try a different starting point.

$\mathbf{ifail}=3$

The routine cannot compute a sufficiently accurate estimate of the gradient vector at the user-supplied starting point. This usually occurs if either the initial parameter estimates are very close to the ML parameter estimates, or you have supplied a very poor estimate of $\Sigma$, or the starting point is very close to the boundary of the stationarity or invertibility region. To proceed you must try a different starting point.

$\mathbf{ifail}=4$

There have been maxcal log-likelihood evaluations made in the routine.

If steady increases in the log-likelihood function were monitored up to the point where this exit occurred, then the exit probably simply occurred because maxcal was set too small, so the calculations should be restarted from the final point held in par and qq. This type of exit may also indicate that there is no maximum to the likelihood surface. Output quantities were computed at the final point held in par and qq, except that if g or cm could not be computed, in which case they are set to zero.

$\mathbf{ifail}=5$

The conditions for a solution have not all been met, but a point at which the log-likelihood took a larger value could not be found.

Provided that the estimated first derivatives are sufficiently small, and that the estimated condition number of the second derivative (Hessian) matrix, as printed when $\mathbf{iprint}\ge 0$, is not too large, this error exit may simply mean that, although it has not been possible to satisfy the specified requirements, the algorithm has in fact found the solution as far as the accuracy of the machine permits.

Such a condition can arise, for instance, if cgetol has been set so small that rounding error in evaluating the likelihood function makes attainment of the convergence conditions impossible.

If the estimated condition number at the final point is large, it could be that the final point is a solution but that the smallest eigenvalue of the Hessian matrix is so close to zero at the solution that it is not possible to recognize it as a solution. Output quantities were computed at the final point held in par and qq, except that if g or cm could not be computed, in which case they are set to zero.

$\mathbf{ifail}=6$

The ML solution is so close to the boundary of either the stationarity region or the invertibility region that g13ddf cannot evaluate the Hessian matrix. The elements of cm are set to zero, as are the elements of g. All other output quantities are correct.

$\mathbf{ifail}=7$

An estimate of the second derivative matrix and the gradient vector at the solution point was computed. Either the Hessian matrix was found to be too ill-conditioned to be evaluated accurately or the gradient vector could not be computed to an acceptable degree of accuracy. The elements of cm are set to zero, as are the elements of g. All other output quantities are correct.

$\mathbf{ifail}=8$

The second-derivative matrix at the solution point is not positive definite. The elements of cm are set to zero. All other output quantities are correct.

$\mathbf{ifail}=-99$

An unexpected error has been triggered by this routine. Please contact NAG.

See Section 3.9 in How to Use the NAG Library and its Documentation for further information.

$\mathbf{ifail}=-399$

Your licence key may have expired or may not have been installed correctly.

See Section 3.8 in How to Use the NAG Library and its Documentation for further information.

$\mathbf{ifail}=-999$

Dynamic memory allocation failed.

See Section 3.7 in How to Use the NAG Library and its Documentation for further information.

7

Accuracy

8

Parallelism and Performance

9

Further Comments

9.1

Memory Usage

9.2

Timing

9.3

Constraining for Stationarity and Invertibility

9.4

Over-parameterisation

9.5

Standardizing the Residual Series

9.6

Assessing the Fit of the Model

10

Example

10.1

Program Text

Program Text (g13ddfe.f90)

10.2

Program Data

Program Data (g13ddfe.d)

10.3

Program Results

Program Results (g13ddfe.r)