g13dd 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 method also calculates estimates of the residual series.
Public Shared Sub g13dd ( _ k As Integer, _ n As Integer, _ ip As Integer, _ iq As Integer, _ mean As Boolean, _ par As Double(), _ qq As Double(,), _ w As Double(,), _ parhld As Boolean(), _ exact As Boolean, _ iprint As Integer, _ cgetol As Double, _ maxcal As Integer, _ ishow As Integer, _ <OutAttribute> ByRef niter As Integer, _ <OutAttribute> ByRef rlogl As Double, _ v As Double(,), _ g As Double(), _ cm As Double(,), _ <OutAttribute> ByRef ifail As Integer _ )
public: static void g13dd( int k, int n, int ip, int iq, bool mean, array<double>^ par, array<double,2>^ qq, array<double,2>^ w, array<bool>^ parhld, bool exact, int iprint, double cgetol, int maxcal, int ishow, [OutAttribute] int% niter, [OutAttribute] double% rlogl, array<double,2>^ v, array<double>^ g, array<double,2>^ cm, [OutAttribute] int% ifail )
static member g13dd : k : int * n : int * ip : int * iq : int * mean : bool * par : float * qq : float[,] * w : float[,] * parhld : bool * exact : bool * iprint : int * cgetol : float * maxcal : int * ishow : int * niter : int byref * rlogl : float byref * v : float[,] * g : float * cm : float[,] * ifail : int byref -> unit
- Type: System..::..Int32On entry: , the number of observed time series.Constraint: .
- Type: System..::..Int32On entry: , the number of observations in each time series.
- Type: System..::..Int32On entry: , the number of AR parameter matrices.Constraint: .
- Type: System..::..Int32On entry: , the number of MA parameter matrices.Constraint: .is not permitted.
- Type: System..::..BooleanOn entry: , if components of have been estimated and , if all elements of are to be taken as zero.Constraint: or .
- Type: array<System..::..Double>()On entry: initial parameter estimates read in row by row in the order , .Thus,
The first elements of par must satisfy the stationarity condition and the next elements of par must satisfy the invertibility condition.If in doubt set all elements of par to .On exit: if or then all the elements of par will be overwritten by the latest estimates of the corresponding ARMA parameters.
- if , must be set equal to an initial estimate of the th element of , for , and ;
- if , must be set equal to an initial estimate of the th element of , and ;
- if , should be set equal to an initial estimate of the th component of (). (If you set to then g13dd will calculate the mean of the th series and use this as an initial estimate of .)
- Type: array<System..::..Double,2>[,](,)[,][,]An array of size [dim1, k]Note: dim1 must satisfy the constraint:On entry: must be set equal to an initial estimate of the th element of . 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 , for and , then g13dd will calculate the covariance matrix between the time series and use this as an initial estimate of .On exit: if or then will contain the latest estimate of the th element of . The lower triangle only is returned.
- Type: array<System..::..Double,2>[,](,)[,][,]An array of size [dim1, n]Note: dim1 must satisfy the constraint:On entry: must be set equal to the th component of , for and .
- Type: array<System..::..Boolean>()On entry: must be set to true if is to be held constant at its input value and false if is a free parameter, for .If in doubt try setting all elements of parhld to false.
- Type: System..::..Boolean
- Type: System..::..Int32On entry: the frequency with which the automatic monitoring method is to be called.
- The ML search procedure is monitored once every iprint iterations and just before exit from the search method.
- The search method is monitored once at the final point.
- The search method is not monitored at all.
- Type: System..::..DoubleOn entry: the accuracy to which the solution in par and qq is required.If cgetol is set to and on exit or , then all the elements in par and qq should be accurate to approximately decimal places. For most practical purposes the value 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.
- Type: System..::..Int32On entry: the maximum number of likelihood evaluations to be permitted by the search procedure.Suggested value: .Constraint: .
- Type: System..::..Int32On 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).Constraint: .
- Type: System..::..Int32%On exit: if or then niter contains the number of iterations performed by the search method.
- Type: System..::..Double%
- Type: array<System..::..Double,2>[,](,)[,][,]An array of size [dim1, n]Note: dim1 must satisfy the constraint:
- Type: array<System..::..Double>()
- Type: array<System..::..Double,2>[,](,)[,][,]An array of size [dim1, npar]Note: dim1 must satisfy the constraint:On exit: if or then will contain an estimate of the correlation coefficient between the th and th elements in the par array for , . If , then will contain the estimated standard error of . If the th component of par has been held constant, i.e., was set to true, then the 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.
Let , for , denote a vector of time series which is assumed to follow a multivariate ARMA model of the form
where , for , is a vector of residual series assumed to be Normally distributed with zero mean and positive definite covariance matrix . The components of are assumed to be uncorrelated at non-simultaneous lags. The and are by matrices of parameters. , for , are called the autoregressive (AR) parameter matrices, and , for , the moving average (MA) parameter matrices. The parameters in the model are thus the ( by ) -matrices, the ( by ) -matrices, the mean vector, , and the residual error covariance matrix . Let
where denotes the by identity matrix.
The ARMA model (1) is said to be stationary if the eigenvalues of lie inside the unit circle. Similarly, the ARMA model (1) is said to be invertible if the eigenvalues of 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 .
You have the option of setting a parameter (exact to false) so that g13dd calculates conditional maximum likelihood estimates (conditional on ). This may be useful if the exact maximum likelihood estimates are close to the boundary of the invertibility region.
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
Note: g13dd may return useful information for one or more of the following detected errors or warnings.
Errors or warnings detected by the method:
Some error messages may refer to parameters that are dropped from this interface (KMAX, LDCM) In these cases, an error in another parameter has usually caused an incorrect value to be inferred.
On entry, , or , or , or , or , where if or if , or , or , or , or ,
- On entry, either the initial estimate of is not positive definite, or the initial estimates of the AR parameters are such that the model is non-stationary, or the initial estimates of the MA parameters are such that the model is non-invertible. To proceed, you must try a different starting point.
- The method 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 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.
- There have been maxcal log-likelihood evaluations made in the method. 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 (as described in [Parameters]) are computed at the final point held in par and qq, except that if g or cm cannot be computed, in which case they are set to zero.
- 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 , 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 (as described in [Parameters]) are computed at the final point held in par and qq, except that if g or cm cannot be computed, in which case they are set to zero.
- This is an unlikely exit, which could occur in e04xaa, which computes an estimate of the second derivative matrix and the gradient vector at the solution point. 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. In this case the elements of cm will be set to zero on exit as will the elements of g. All other output quantities will be correct.
- The second derivative matrix at the solution point is not positive definite. In this case the elements of cm will be set to zero on exit. All other output quantities will be correct.
- Internal memory allocation failed.
On exit from g13dd, if or and cgetol has been set to , then all the parameters should be accurate to approximately decimal places. If cgetol was set equal to a value less than the machine precision, , then all the parameters should be accurate to approximately .
If on exit (i.e., maxcal likelihood evaluations have been made but the convergence conditions of the search method have not been satisfied), then the elements in par and qq may still be good approximations to the ML estimates. Inspection of the elements of g may help you determine whether this is likely.
Let and . Local workspace arrays of fixed lengths are allocated internally by g13dd. The total size of these arrays amounts to integer elements and real elements.
The number of iterations required depends upon the number of parameters in the model and the distance of the user-supplied starting point from the solution.
If the solution lies on the boundary of the admissibility region (stationarity and invertibility region) then g13dd may get into difficulty and exit with . If this exit occurs you are advised to either try a different starting point or a different setting for exact. If this still continues to occur then you are urged to try fitting a more parsimonious model.
You are advised to try and avoid fitting models with an excessive number of parameters since over-parameterisation can cause the maximization problem to become ill-conditioned.
The standardized estimates of the residual series (denoted by ) can easily be calculated by forming the Cholesky decomposition of , e.g., and setting . (F07FDF not in this release) may be used to calculate the array g. The components of which are now uncorrelated at all lags can sometimes be more easily interpreted.
If your time series model provides a good fit to the data then the residual series should be approximately white noise, i.e., exhibit no serial cross-correlation. An examination of the residual cross-correlation matrices should confirm whether this is likely to be so. You are advised to call (G13DSF not in this release) to provide information for diagnostic checking. (G13DSF not in this release) returns the residual cross-correlation matrices along with their asymptotic standard errors. (G13DSF not in this release) also computes a portmanteau statistic and its asymptotic significance level for testing model adequacy. If or on exit from g13dd then the quantities output k, n, v, kmax, ip, iq, par, parhld, and qq will be suitable for input to (G13DSF not in this release).
This example shows how to fit a bivariate AR(1) model to two series each of length . will be estimated and will be constrained to be zero.