naginterfaces.library.tsa.multi_​varma_​estimate

naginterfaces.library.tsa.multi_varma_estimate(ip, iq, mean, par, qq, w, parhld, exact, iprint, cgetol, ishow, maxcal=None, io_manager=None)

multi_varma_estimate fits a vector autoregressive moving average (VARMA) model to an observed vector of time series using the method of Maximum Likelihood (ML). Standard errors of parameter estimates are computed along with their appropriate correlation matrix. The function also calculates estimates of the residual series.

For full information please refer to the NAG Library document for g13dd

https://support.nag.com/numeric/nl/nagdoc_30.2/flhtml/g13/g13ddf.html

Parameters

ipint

$p$, the number of AR parameter matrices.

iqint

$q$, the number of MA parameter matrices.

meanbool

$mean=True$, if components of $\mu$ have been estimated and $mean=False$, if all elements of $\mu$ are to be taken as zero.

parfloat, array-like, shape $\left(\text{npar}\right)$

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 $ip>0$, $par[(\text{l}-1)\times k\times k+(\text{i}-1)\times k+\text{j}-1]$ must be set equal to an initial estimate of the $(\text{i},\text{j})$th element of ${\varphi }_{\text{l}}$, for $\text{j}=1,2,\dots ,k$, for $\text{i}=1,2,\dots ,k$, for $\text{l}=1,2,\dots ,p$;

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 initial estimate of the $(i,j)$th element of ${\theta }_l$, $l=1,2,\dots ,q$ and $i,j=1,2,\dots ,k$;

if $mean=True$, $par[(p+q)\times k\times k+i-1]$ should be set equal to an initial estimate of the $i$th component of $\mu$ ($\mu \left(i\right)$). (If you set $par[(p+q)\times k\times k+i-1]$ to $0.0$ then multi_varma_estimate 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$.

qqfloat, array-like, shape $(k,k)$

$qq[\text{i}-1,\text{j}-1]$ must be set equal to an initial estimate of the $(\text{i},\text{j})$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 $qq[\text{i}-1,\text{j}-1]=0.0$, for $\text{j}=1,2,\dots ,i$, for $\text{i}=1,2,\dots ,k$, then multi_varma_estimate will calculate the covariance matrix between the $k$ time series and use this as an initial estimate of $\Sigma$.

wfloat, array-like, shape $(k,n)$

$w[\text{i}-1,\text{t}-1]$ must be set equal to the $\text{i}$th component of $W_{\text{t}}$, for $\text{t}=1,2,\dots ,n$, for $\text{i}=1,2,\dots ,k$.

parhldbool, array-like, shape $\left(\text{npar}\right)$

$parhld[\text{i}-1]$ must be set to $True$ if $par[\text{i}-1]$ is to be held constant at its input value and $False$ if $par[\text{i}-1]$ is a free parameter, for $\text{i}=1,2,\dots ,\text{npar}$.

If in doubt try setting all elements of $parhld$ to $False$.

exactbool

Must be set equal to $True$ if you wish multi_varma_estimate to compute exact maximum likelihood estimates. $exact$ must be set equal to $False$ if only conditional likelihood estimates are required.

iprintint

The frequency with which the automatic monitoring function is to be called.

$iprint>0$

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

$iprint=0$

The search function is monitored once at the final point.

$iprint<0$

The search function is not monitored at all.

cgetolfloat

The accuracy to which the solution in $par$ and $qq$ is required.

If $cgetol$ is set to ${10}^{-l}$ and on exit no exception or warning is raised or $errno$ >= 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, $ϵ$, then multi_varma_estimate will use the value $10.0\times \sqrt{ϵ}$ instead.

ishowint

Specifies which of the following two quantities are to be printed.

1. table of maximum likelihood estimates and their standard errors (as returned in the output arrays $par$, $qq$ and $cm$);

2. table of residual series (as returned in the output array $v$).

$ishow=0$

None of the above are printed.

$ishow=1$

(i) only is printed.

$ishow=2$

(i) and (ii) are printed.

maxcalNone or int, optional

Note: if this argument is None then a default value will be used, determined as follows: $40\times \text{npar}\times (\text{npar}+5)$.

The maximum number of likelihood evaluations to be permitted by the search procedure.

io_managerFileObjManager, optional

Manager for I/O in this routine.

Returns

parfloat, ndarray, shape $\left(\text{npar}\right)$

If no exception or warning is raised or $errno$ >= 4 then all the elements of $par$ will be overwritten by the latest estimates of the corresponding ARMA parameters.

qqfloat, ndarray, shape $(k,k)$

If no exception or warning is raised or $errno$ >= 4 then $qq[i-1,j-1]$ will contain the latest estimate of the $(i,j)$th element of $\Sigma$. The lower triangle only is returned.

niterint

If no exception or warning is raised or $errno$ >= 4 then $niter$ contains the number of iterations performed by the search function.

rloglfloat

If no exception or warning is raised or $errno$ >= 4 then $rlogl$ contains the value of the log-likelihood function corresponding to the final point held in $par$ and $qq$.

vfloat, ndarray, shape $(k,n)$

If no exception or warning is raised or $errno$ >= 4 then $v[\text{i}-1,\text{t}-1]$ will contain an estimate of the $\text{i}$th component of ${ϵ}_{\text{t}}$, for $\text{t}=1,2,\dots ,n$, for $\text{i}=1,2,\dots ,k$, corresponding to the final point held in $par$ and $qq$.

gfloat, ndarray, shape $\left(\text{npar}\right)$

If no exception or warning is raised or $errno$ >= 4 then $g[i-1]$ 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.

cmfloat, ndarray, shape $(\text{npar},\text{npar})$

If no exception or warning is raised or $errno$ >= 4 then $cm[i-1,j-1]$ 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 \text{npar}$, $1\le j\le \text{npar}$. If $i=j$, then $cm[i-1,j-1]$ will contain the estimated standard error of $par[i-1]$. If the $l$th component of $par$ has been held constant, i.e., $parhld[l-1]$ 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.

Raises

NagValueError

(errno $1$)

On entry, $k=⟨value⟩$.

Constraint: $k\ge 1$.

(errno $1$)

On entry, $\text{npar}=⟨value⟩$.

Constraint: $\text{npar}\ge 0$.

(errno $1$)

On entry, $n=⟨value⟩$, $k=⟨value⟩$ and $\text{npar}=⟨value⟩$.

Constraint: $n\times k>\text{npar}+k\times (k+1)/2$.

(errno $1$)

On entry, $ip=⟨value⟩$.

Constraint: $ip\ge 0$.

(errno $1$)

On entry, $iq=⟨value⟩$.

Constraint: $iq\ge 0$.

(errno $1$)

On entry, $ip=0$ and $iq=0$.

(errno $1$)

On entry, $\text{npar}=⟨value⟩$.

Constraint: $\text{npar}=⟨value⟩$.

(errno $1$)

On entry, $maxcal=⟨value⟩$.

Constraint: $maxcal\ge 1$.

(errno $1$)

On entry, $ishow=⟨value⟩$.

Constraint: $0\le ishow\le 2$.

(errno $2$)

The initial estimate of $\Sigma$ is not positive definite.

(errno $2$)

The initial AR parameter estimates are outside the stationarity region.

(errno $2$)

The initial MA parameter estimates are outside the invertibility region.

(errno $2$)

The starting point is too close to the boundary of the admissibility region.

(errno $3$)

The function cannot compute a sufficiently accurate estimate of the gradient vector at the user-supplied starting point.

(errno $4$)

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

See the function document for further details.

Warns

NagAlgorithmicWarning

(errno $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.

See the function document for further details.

(errno $6$)

The ML solution is so close to the boundary of either the stationarity region or the invertibility region that multi_varma_estimate 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.

(errno $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.

(errno $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.

Notes

Let $W_{\text{t}}={(w_{1\text{t}},w_{2\text{t}},\dots ,w_{\text{k}\text{t}})}^T$, for $\text{t}=1,2,\dots ,n$, denote a vector of $k$ time series which is assumed to follow a multivariate ARMA model of the form

$$\begin{array}{c}\begin{array}{cc}{W}_t-\mu =& {\varphi }_1(W_{t-1}-\mu )+{\varphi }_2(W_{t-2}-\mu )+\cdots +{\varphi }_p(W_{t-p}-\mu )\\ & +{ϵ}_t-{\theta }_1{ϵ}_{t-1}-{\theta }_2{ϵ}_{t-2}-\cdots -{\theta }_q{ϵ}_{t-q}\end{array}\end{array}$$

where ${ϵ}_{\text{t}}={({ϵ}_{1\text{t}},{ϵ}_{2\text{t}},\dots ,{ϵ}_{k\text{t}})}^T$, for $\text{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 ${ϵ}_t$ are assumed to be uncorrelated at non-simultaneous lags. The ${\varphi }_i$ and ${\theta }_j$ are $k\times k$ matrices of parameters. $\left\{{\varphi }_{\text{i}}\right\}$, for $\text{i}=1,2,\dots ,p$, are called the autoregressive (AR) parameter matrices, and $\left\{{\theta }_{\text{i}}\right\}$, for $\text{i}=1,2,\dots ,q$, the moving average (MA) parameter matrices. The parameters in the model are thus the $p$ ($k\times k$) $\varphi$-matrices, the $q$ ($k\times k$) $\theta$-matrices, the mean vector, $\mu$, and the residual error covariance matrix $\Sigma$. Let

$$\begin{array}{c}A\left(\varphi \right)={\left[\begin{array}{ccccccc}{\varphi }_1& I& 0& .& .& .& 0\\ {\varphi }_2& 0& I& 0& .& .& 0\\ .& & & .& & & \\ .& & & & .& & \\ .& & & & & .& \\ {\varphi }_{p-1}& 0& .& .& .& 0& I\\ {\varphi }_p& 0& .& .& .& 0& 0\end{array}\right]}_{pk\times pk}\text{and}B\left(\theta \right)={\left[\begin{array}{ccccccc}{\theta }_1& I& 0& .& .& .& 0\\ {\theta }_2& 0& I& 0& .& .& 0\\ .& & & .& & & \\ .& & & & .& & \\ .& & & & & .& \\ {\theta }_{q-1}& 0& .& .& .& .& I\\ {\theta }_q& 0& .& .& .& .& 0\end{array}\right]}_{qk\times qk}\end{array}$$

where $I$ denotes the $k\times k$ identity matrix.

The ARMA model (1) is said to be stationary if the eigenvalues of $A\left(\varphi \right)$ lie inside the unit circle. Similarly, the ARMA model (1) is said to be invertible if the eigenvalues of $B\left(\theta \right)$ lie inside the unit circle.

The method of computing the exact likelihood function (using a Kalman filter algorithm) is discussed in Shea (1987). A quasi-Newton algorithm (see Gill and Murray (1972)) is then used to search for the maximum of the log-likelihood function. Stationarity and invertibility are enforced on the model using the reparameterisation discussed in Ansley and Kohn (1986). Conditional on the maximum likelihood estimates being equal to their true values the estimates of the residual series are uncorrelated with zero mean and constant variance $\Sigma$.

You have the option of setting an argument ($exact$ to $False$) so that multi_varma_estimate calculates conditional maximum likelihood estimates (conditional on $W_0=W_{-1}=\cdots =W_{1-p}={ϵ}_0={ϵ}_{-1}=\cdots ={ϵ}_{1-q}=0$). This may be useful if the exact maximum likelihood estimates are close to the boundary of the invertibility region.

You also have the option (see Parameters) of requesting multi_varma_estimate to constrain elements of the $\varphi$ and $\theta$ matrices and $\mu$ vector to have pre-specified values.

References

Ansley, C F and Kohn, R, 1986, A note on reparameterising a vector autoregressive moving average model to enforce stationarity, J. Statist. Comput. Simulation (24), 99–106

Gill, P E and Murray, W, 1972, Quasi-Newton methods for unconstrained optimization, J. Inst. Math. Appl. (9), 91–108

Shea, B L, 1987, Estimation of multivariate time series, J. Time Ser. Anal. (8), 95–110