An Algorithm for the Fisher Information Matrix of a VARMAX Process

: In this paper, an algorithm for Mathematica is proposed for the computation of the asymptotic Fisher information matrix for a multivariate time series, more precisely for a controlled vector autoregressive moving average stationary process, or VARMAX process. Meanwhile, we present brieﬂy several algorithms published in the literature and discuss the sufﬁcient condition of invertibility of that matrix based on the eigenvalues of the process operators. The results are illustrated by numerical computations.


Introduction
In this paper, we present an algorithm for the symbolic software package Mathematica to compute the asymptotic Fisher information matrix for a multivariate time series, more precisely a controlled vector autoregressive moving average stationary process, or VAR-MAX process.There are several ways to obtain the Fisher information matrix for time series models.The most used procedure is computing the Hessian of the Gaussian log-likelihood when a model is fitted to the time series, usually after non-linear optimization (since the model is non-linear with respect to the parameters except for the pure VAR case), e.g., [1].Another much less used method consists of determining the exact (Gaussian) Fisher information, see [2][3][4][5][6][7][8].This is nice when the data are available but, as [9] indicates, it can be useful to know the matrix at a tentative parameter point before obtaining the data, in order to determine the series length necessary to reach some specific accuracy for the estimators of the parameters.Then, an algorithm for the computation of the so-called asymptotic Fisher information matrix at a given parameter point is what is needed.It was provided in the form of an algebraic expression for a handful of small-order (1 or 2) ARMA models by Box and Jenkins in 1970, see [10].Otherwise, there is a vast literature for scalar models, i.e., simple ARMA models, e.g., [11,12], seasonal ARMA models [13], ARMAX models [14], single input single output models [15], and for vector models, i.e., VARMA models [16,17] and VARMAX models [18,19].
We introduce a new algorithm well-suited to symbolic computation software packages (here Mathematica) and we present an illustrative example.Then, we discuss the invertibility of the Fisher information matrices and provide other examples.
The proposed algorithm differs from the works of Zadrovny [3,20], see also Zadrozny and Mittnick [21], who developed their methods for VARMAX models, in the sense that they are based on an approximation of the exact Fisher information.Here, we follow Whittle's approach [22] of expressing the asymptotic information matrix by computing a circular integral of a rational function, related to the autoregressive and moving average operators.That was the way used by Newton for VARMA models, see [16], with the difference that we make the effective computation and that we consider VARMAX models.The more fundamental source is [19], who obtained an explicit expression for the Fisher information matrix of a VARMAX model under the form of a sum of two integrals while [18] contains expressions for each element of that matrix.As a matter of fact, we first need to generalize [19], which is developed for the special case where the number of explanatory variables equals the number of variables of interest and, moreover, the coefficient matrix at lag 0 is the identity matrix.These restrictions were explained by the revealed objective of the paper which was to exploit tensor Sylvester matrices and try to extend to VARMAX models the resultant property of the Fisher information matrix wrongly claimed for VARMA models by [17].Indeed, in [17] (p.1986), it is written "[a] representation [based on tensor Sylvester matrices] is more appropriate . . . to prove a possible resultant property of the Fisher information matrix of a VARMAX process".Such a tensor Sylvester matrix is associated with two matrix polynomials and it becomes singular if, and only if, the two matrix polynomials have at least one common eigenvalue, see [23].
In [17], it is said that the Fisher information matrix of a VARMA process becomes invertible if, and only if, the tensor Sylvester matrix is invertible, in other words, if, and only if, the autoregressive and moving average matrix polynomials of the VARMA process have no common eigenvalue.This is called the Sylvester resultant property.In [17], the Sylvester resultant property is shown for a scalar ARMAX process but it is no longer true, in general, for other than scalar ARMA or ARMAX processes.Indeed, in [24] Mélard indicated the error in the claim and explained why the example in [17] was wrong.In the present paper, we correct the assertion stated in [17] for VARMA processes, and we extend the "if" part to a class of VARMAX processes.As will be shown, the "only if" part of that result is, however, wrong when the dimension of the process is larger than 1.We will also compare our algorithms with other algorithms proposed for special cases in the literature.Although the results are less powerful, they are what is needed in practice: a sufficient condition of invertibility of the Fisher information matrix, indicating a possible lack of identifiability so that parameter estimation is risky.A necessary condition that should involve more information about the matrix polynomials is not as useful.
Consider the vector stochastic difference equation representation of a linear system of order (p, r, q) of the process {y t , t ∈ Z}, Z is the set of integers where y t , x t , and ε t are, respectively, the n-dimensional stochastic observed output, the mdimensional observed input, and the n-dimensional unobserved errors, and where A j ∈ R n×n , j = 1, . . ., p, C j ∈ R n×m , j = 0, . . ., r − 1, and B j ∈ R n×n , j = 1, . . ., q, are associated parameter matrices.We additionally assume A 0 = B 0 = I n , the n × n identity matrix.We also assume that C 0 is an invertible matrix.In the examples, we will sometimes take C 0 = I n fixed instead of being a matrix of parameters so that the maximum lag r − 1 is replaced by r.Note that the absence of lag induced by ( 1) is purely conventional.For example, a lagged effect of x t on y t can be produced by defining x t = x t−1 .More precisely, we suppose that x t is stochastic and that (y t , x t ) is a Gaussian stochastic process.The error {ε t , t ∈ Z} is a sequence of independent zero-mean n-dimensional Gaussian random variables with a strictly positive definite covariance matrix Σ.We denote this by Σ 0. We denote the transposition by T , the complex conjugate transposition by * , and the mathematical expectation is E. We assume E{x t ε T t } = 0, for all t, t .We denote z as the backward shift operator, for example, zx t = x t−1 .Then (1) can be written as where are the associated matrix polynomials, where z ∈ C (the duplicate use of z as an operator and as a complex variable is usual in the signal and time series literature, e.g., [25][26][27]).The assumptions det(A(z)) = 0 for all |z| ≤ 1 (causality) and det(B(z)) = 0 for all |z| ≤ 1 (invertibility) are imposed so that the eigenvalues of the matrix polynomials A(z) and B(z) will be outside the unit circle, see [25].Remind that the eigenvalues of an n × n matrix polynomial D(z) are the roots of the equation det(D(z)) = 0 (see [28], p. 14).We will also use the reciprocal polynomial of D(z), with degree d, say.It is is full-rank, there are dn eigenvalues for D(z).The eigenvalues of D(z) are the inverse of those of D(z), provided the missing eigenvalues of the former are treated as being at ∞.
Remark 1.Note that there are restrictions in the VARMAX model being considered.The dimensions of the unobserved errors ε t and the observed output y t are the same.This is often the case in the literature, although [29], for example, consider a VARMA model where the dimension of the unobserved errors is smaller than that of the observed output.Except when we will mention the tensor Sylvester matrices, the dimensions of the observed input x t and the observed output y t need not be the same, like [7,18] but contrary to [19].
We store the VARMAX (p, r, q) coefficients in an l = (n 2 (p + q) + nmr) × 1 vector ϑ defined as follows: The vec operator transforms a matrix into a vector by stacking the columns of the matrix one underneath the other, according to vec X = col(col(X ij ) n i=1 ) n j=1 , see, e.g., [30].The observed input variable x t is assumed to be a stationary m-dimensional Gaussian VARMA process with white noise process η t satisfying E{η t η T t } = Ω 0 and where a(z) and b(z) are, respectively, the autoregressive and moving average matrix polynomials, such that a(0) = b(0) = I m and det(a(z)) = 0 far all |z| ≤ 1 and det(b(z)) = 0 for all |z| ≤ 1.The spectral density of process x t is defined as, see, e.g., [26] where i is the standard imaginary unit, f is the frequency, the spectral density matrix R x (e i f ) is Hermitian, and we further have, R x (e i f ) 0 and Therefore, there is at least one solution of (1) which is stationary.

The Fisher Information Matrix of VARMAX Processes
Using (2) for all t ∈ Z enables us to determine the residual ε t (ϑ) which depends on the parameter vector, ϑ.Under the above assumptions, the asymptotic Fisher information matrix, F (ϑ), is defined by the following (l × l) matrix, which does not depend on t: where the (n × l) matrix ∂(•)/∂ϑ T is the derivative with respect to ϑ T .That derivative is computed in [19], using the Kronecker product ⊗, to obtain: For each positive integer k, denote u k (z) = (1, z, z 2 , . . ., z k−1 ) T .Let us define: ) and the square matrices of dimension (p + q)n + rm using the Hermitian spectral density matrix, R x (z), defined in (4).Then, as shown in [19] (for the special case where m = n), the Fisher information matrix is given by: The integrals are counter-clockwise.It is easy to check that the result in (Proposition 2.1, [19]) is valid even if m = n and C(0) is not the identity matrix.The main changes are in (7) but also in the middle block of (Equations ( 18)-( 20), [19]) where I n 2 has to be replaced by I nm .Note that the subsequent sections of [19] are not valid if m = n.

New Algorithm and Comparison with Other Algorithms
We propose an algorithm based on the method in the previous section and an implementation as a notebook for Mathematica.To make use of exact computations, the coefficients of the VARMAX model as well as Σ and R x are entered as fractions, if possible.The algorithm is provided in Figure 1 using some notations that are detailed here.
The arguments (z) are written [z] for the polynomials.Given the polynomial matrices A(z), B(z), and C(z) (hence denoted A[z], B[z], and C[z], respectively, in Figure 1), the vectors u p (z), u q (z), and u r (z) (denoted up[z], uq[z], and ur[z], respectively), the zero matrices O rm×n and O qn×m (denoted Orm$n and Oqn$m, respectively), the identity matrices I n and I m (denoted In and Imm, respectively, noting that Im is a locked symbol), plus Σ(z) (denoted Sigma[z]) and R x (z) (denoted Rx[z]).Then, the block matrices G(z) (Gcal[z]) and K(z) (Kcal[z]), are computed using the ArrayFlatten command of Mathematica, plus KroneckerProduct and Inverse, then σ(z) (sigma[z]), P (z) (Pcal[z]), and Q(z) (Qcal[z]) using also the Transpose function.Since Mathematica does not treat well circular integrals, we transform them into line integrals between 0 and 2π by using the change of variable z = e i f , 0 ≤ f ≤ 2π, hence dz/z = id f .We compute the two transformed integrals of matrices using the command Integrate and sum them up.Contrary to the computations in (Appendix B, [7]), where individual formulas (taken from [18]) were used for each element and computed using the Cauchy rule, here the (two) integrals of the whole matrices are computed symbolically.To simplify the expressions as much as possible as rational functions of z, we use several times the commands Together and also Simplify.The algorithm is shown in Figure 1 as a Mathematica notebook.The example is related to a VARMAX model used by [31] and also used by [32].The model is defined by: The coefficients are entered as fractions, hence Part of the results are thus shown in Figures 2-4.It took about 7.2 min to produce the output.Making the computation using the numerical program described in [32] on the same computer took 30 s but for one million evaluations.Supplementary Materials Supplementary S1 contains the corresponding notebook (Figure1234.nb)that can be read using the free Wolfram's CDF Player and its PDF copy (Figure1234.pdf).We will show and discuss other examples in the next Section.Note that we tried to implement the algorithm in two open-source packages for symbolic computation, Octave and Maxima.It is feasible because both programs allow us to integrate matrices but it appears that there is an error in the computation of the integrals being involved.We have informed the respective developers of the bug but the errors are not corrected yet.
Before showing other examples, let us look at other algorithms published in the literature.As said already, most of them are for the exact Fisher information, not the asymptotic Fisher information and we will not discuss them here, focusing instead on the asymptotic information matrix.Note that none of the existing algorithms is developed for VARMAX models, but only for VARMA or univariate models.(10).Fcal1 and Fcal2 are, respectively, the first and second terms in the right-hand side of (9) and Fcal is their sum, so equal to F (ϑ).
The first algorithm was due to Newton [16] for VARMA models.It is based on [22] who expressed the information matrix by a trigonometric integral involving the derivatives of the spectral density matrix of the process with respect to the parameters.The spectral density at frequency f , f ∈ [−π, π], of a VARMAX process, e.g., (1) without the term ∑ r−1 j=0 C j x t−j , is given by (1/2π)A −1 (e i f )B(e i f )ΣB * (e i f )A * −1 (e i f ).Newton [16] does not indicate how to evaluate these matrix integrals that he obtained.It can be seen that the integrand of each element of these matrix integrals is a rational function of e i f .Furthermore, the condition of stationarity implies that the denominator has a factor with zeroes in the unit circle and zeroes outside of it.Peterka and Vidinčev [33] have proposed an algorithm for evaluating circular integrals of a rational function around the unit circle and [34] has proposed an implementation that consists of a series of reduction in the degrees of the involved polynomials.That method was used by [12] for the evaluation of the Fisher information matrix of ARMA models.
A completely different early and effective method but restricted to ARMA models was proposed by Godolphin and Unwin [11].It makes use of simple matrix operations with the inversion of two matrices, one of order p and one of order q.
Another approach was proposed by [35] which consists of replacing the computation of the integral by the evaluation of the cross-covariance of two ARMA processes based on the same white noise process, and this problem can be transformed in the computation of the autocovariance of a single AR process.There are excellent algorithms for that computation due to [36,37].They require a number of operations which is quadratic in the AR order, much less than previous algorithms where a matrix of that order had to be inverted.The first algorithm was used by [13] for the evaluation of the Fisher information matrix of seasonal ARMA models.Ref. [38] have complained that several high-order AR models are sometimes involved and have proposed an alternative approach.Nevertheless, the algorithm making use of the autocovariances of an AR process using [37] could be developed for SISO models [15] and even for seasonal SISO models [32].For ARMAX models, ref. [14] did not insist on algorithms but more on checking the invertibility of the Fisher information matrix, the subject of the next section.For ARMAX or SISO models, the Fisher information matrix is a sum of two integrals.The approach has however reached its limitations because the construction of the polynomials involved in the procedure becomes complex.For instance, the algorithm in [32] made use of a table with letters identifying the regular polynomials (a to f) and the seasonal polynomials (A to F) among others (the constant and the regression coefficient) involved as numerators or denominators in the different integrals.For a VARMAX model, there are 3 polynomials, hence the information matrix is a 3 × 3 block matrix with, according to symmetry, 6 different blocks, including 1 corresponding to the parameters in the A j where a sum of two integrals appears.A simple inspection of the expressions for the different blocks in [18] reveals that obtaining the polynomials is a complex task.

A Sufficient Condition of Identifiability
The conditions for the identifiability of a VARMAX model have been studied by [25], among others.They are repeated later.A clear indication of identifiability is the invertibility of the Fisher information matrix.In this section, we derive a sufficient condition of invertibility for the case where m = n and we enforce that the condition is not necessary, as could be deduced from the wrong arguments of [17].For an ARMAX model, i.e., when n = m = 1, [14] proved that a necessary and sufficient condition for invertibility of the Fisher information matrix is that the three polynomials have no common root.Indeed, as will be seen, the result of [14] cannot be extended to vector processes.
We have already been reminded of the concepts of reciprocal polynomials and eigenvalues of matrix polynomials.There remains to define what are a unimodular polynomial matrix and the common left divisor of matrix polynomials, see, e.g., [39].
A unimodular polynomial matrix U(z) is a polynomial square matrix such that its determinant is a non-zero constant, instead of being a polynomial in z.Consequently, U −1 (z) is also a polynomial matrix.In general, the inverse of a square polynomial matrix is a matrix with rational elements, not polynomial elements.
Three matrix polynomials A(z), B(z), and C(z) with the same number of rows n have a common left divisor if there exist n × n polynomial matrices F(z), A (z), B (z), and C (z), such that A(z) = F(z)A (z), B(z) = F(z)B (z), and C(z) = F(z)C (z).In matrix form, we can write (A(z) B(z) C(z)) = F(z)(A (z) B (z) C (z)).Then (A(z) B(z) C(z)) is called a right multiple of F(z).A left divisor F(z) is called the greatest common left divisor of (A(z) B(z) C(z)) if it is a right multiple of all left divisors.Multiplying a greatest common left divisor to the right by any unimodular polynomial matrix yields another greatest common left divisor.As indicated by (Section 2.2, [25]), a greatest common left divisor can be constructed by elementary column operations: interchange any two columns, multiply any column by a real number different from 0, and add a polynomial multiple of any column to any column.Also, it corresponds to the right multiplication of (A(z) B(z) C(z)) by an appropriate unimodular matrix.The concept can also be defined for rectangular matrices with the same number of rows although we will not consider that generalization.Lemma 1. Assume that A(z), B(z), and C(z) have the prescribed degrees, respectively, p, q, and r − 1, with det(A(z)) = 0 and det(B(z)) = 0, for z in the unit circle of C, and Σ is non-singular.Assume also that x t has an absolutely continuous spectrum with spectral-density non-zero on a set of positive measure on (−π, π].Then a necessary and sufficient condition of identifiability of a stationary VARMAX process satisfying (1) is (i) A(z), B(z), and C(z) have I n as greatest common left divisor, (ii) the matrix (A p B q C r−1 ) has rank n.
Proof.Since Σ is non-singular and given the assumptions on A(z) and B(z), the conditions (8a), (8b), and (8c) of [40] are satisfied for A −1 (z)B(z), and the lemma is a special case of (Theorem 2, [40]) in the case (iii) 2 .
Since identifiability is equivalent to the inversibility of the Fisher information matrix, we are interested in finding a sufficient condition for identifiability.Equivalently, we are looking for a simple necessary condition for the lack of identifiability or the singularity of the Fisher information matrix, when m = n.
We can state the following theorem.
Theorem 1.Under the assumptions of Lemma 1 and if m = n, a sufficient condition of invertibility of the Fisher information matrix, F (ϑ), associated with the VARMAX model with the matrix polynomials A(z), C(z), and B(z) of degree p, r − 1, q, respectively, is that the reciprocal matrix polynomials Ã(z), B(z), and C(z) have no common eigenvalue.
Proof.If the lack of identifiability occurs because (i) in Lemma 1 is not satisfied, that means that there exists a non-unimodular polynomial matrix F(z), i.e., with det(F(z)) different from a constant and matrices A (z), B (z), and C (z), such that (A(z) Since the three polynomial matrices are square, we can consider their determinant stored in a row vector for which where det(F(z)) is a polynomial different from a constant.Hence, the equations det(A(z)) = 0, det(B(z)) = 0, and det(C(z)) = 0 have at least one common root.Consequently, the matrix polynomials A(z), B(z), and C(z) have at least one common eigenvalue.The same is also true for the reciprocal matrix polynomials Ã(z), B(z), and C(z).Now, if the lack of identifiability occurs because (ii) in Lemma 1 is not satisfied, that means that the matrix (A p B q C r−1 ) has a rank strictly smaller than n.Consequently, row n, say, of that matrix is a linear combination of the other n − 1 rows and the matrices A p , B q , and C r−1 do not have full rank.Hence the determinants det(A(z)), det(B(z)), and det(C(z)) are polynomials of degree strictly smaller than, respectively, np, nq, and n(r − 1).Each of the reciprocal matrix polynomials Ã(z), B(z), and C(z) have at least an eigenvalue equal to 0, hence at least one common eigenvalue.
Note that [17] has implicitly assumed that the determinants of A(z) and B(z) have their maximum degrees.For a discussion of identifiability without co-primeness, see [41].For VARMA models, Ref. [24] has shown that it is not true to say, as [17] did, that the conditions in Theorem 1 are also necessary for invertibility of the Fisher information matrix.The approach of [17] used tensor Sylvester matrices which are a generalization of Sylvester matrices for the case of matrix polynomials instead of scalar polynomials.Consequently, the expectations of [19] mentioned in the second paragraph of Section 1 are not met for VARMAX models.Moreover, since there are three matrix polynomials and C(z) does not need to be a square matrix, the tensor Sylvester matrices are not useful.It is also not useful to detail the alternative representations to (6) developed in (Sections 2.4 to 2.7, [19]) using the tensor Sylvester matrices.Let us just say that these representations seem to be correct in the sense that we have checked on all our examples for which m = n that the Fisher information matrix obtained using them is identical to what is obtained using ( 7)-(9).

Numerical Experiments
The approach developed above can be used for any VARMAX model, for instance, the one produced by [42] obtained using a least-squares-based iterative estimation method.
To save space, we will use examples based on the simplest case, i.e., n = m = 2, C 0 = I 2 , and p = r = q = 1.We denote simply A = A 1 , B = B 1 , and C = C 1 .For our first three examples, we will have: Example 1.Let A, B, and C be defined by where a and b are constants.The eigenvalues of A(z), B(z), and C(z) are, respectively, the pairs (0.8, a), (0.6, b), and (0.7, a) so that, whatever a and b, there is a common eigenvalue for A(z) and C(z).Clearly, the model is identifiable except if a = b = 0.8 because then the factor 1 − 0.8z can be simplified on the first row of the system equation.In this case, because of the presence of the symbolic constants a and b, the algorithm needs to be changed by adding a command Assuming, as shown in Figure 5, with the mention that a and b are restricted to the interval [−1, 1].It is of course not possible to show here more of the matrix, see Supplementary Materials.Take b = 0.8 to simplify the discussion.Then, using Mathematica, it can be seen that the Fisher information matrix has a determinant proportional to (4 − 5a) with a strictly positive factor, and indeed it is 0 if, and only if, a = 0.8.If a = 0.8, the 2nd, 6th, and 10th rows of the Fisher information matrix contain the fractions and it is easy to check that −Row2 = −Row6 − Row10.Since the asymptotic information matrix is singular, there is no need to show it.As we said, the computations described in Section 3 were implemented as matrix computations in Mathematica with coefficients given in the fractional form (e.g., 8/10 instead of 0.8) to have exact results.The complete Mathematica notebook (Example51.nb) is provided as Supplementary Materials, with its PDF output (Example51.pdf).
Version July 20, 2023 submitted to Algorithms 11 of 15   This example is a generalization of a VARMA model considered by [24] based on an example in [43].The determinants of A(z), B(z), and C(z) have degree 1, so we need to consider the reciprocal matrix polynomials Ã(z), B(z), and C(z) which have respective roots (−0.6, 0), (−0.5, 0), and (−0.8, 0).Hence, we can consider the singularity of the Fisher information matrix.This is confirmed by This example is a generalization of the bivariate VARMA(1,1) model in (Example 1, [17]).There, the two matrix polynomials had the same two eigenvalues.In [17], it was concluded, wrongly, that the Fisher information matrix be singular although the numerical computations in Matlab lead to the smallest eigenvalue equal to 0.0067.This incoherency was explained by numerical inaccuracy.In [24], Mélard gave a second view to that example and discovered that the necessity and sufficiency of invertibility of the Fisher information matrix is only a necessary condition.
For the VARMAX model related to (11), the three matrix polynomials have the same eigenvalues ±5/ √ 11.Exact computations with Mathematica indicate, however, that the determinant of the 12 × 12 Fisher information matrix is strictly positive and that the smallest eigenvalue is 0.0919.The inverse of the asymptotic Fisher information matrix is shown in Figure 6.This example provides a counter-example that the sufficient condition in Theorem 1 is not necessary.
The total number of parameters is 4 + 6 + 6 + 4 = 20.In (Section 4.1, [7]), the exact Fisher information matrix was computed and compared to those obtained using the E4 Toolbox for Matlab, see [5,44], but (Appendix B, [7]) contains the results for the asymptotic Fisher information matrix obtained using individual formulas in [18] for each of the 10 different blocks, the computation of the integrals being performed using the algorithm in [33].Here, we do not bother with complex individual expressions and the computation of integrals.We simply use the notebook described in Figure 1 and adapt the model specifications.See the resulting information matrix in Figure 7.

Conclusions
In this paper, we have proposed for the first time an exact method to compute the asymptotic Fisher information matrix of VARMAX processes, an algorithm for symbolic computation software, and an implementation for It is based on an extension of the method proposed in [19].
Therefore, we could not compare it with other algorithms in the general case.Nevertheless, we could compare it to the special case of an ARMAX process where there exists alternative algorithms.
The open-source packages Octave and Maxima should be able to do the job since their symbolic integration procedure can work on matrices of functions, but they do not work presently.More developments would be needed to circumvent their failing integration.
We also provide a simple sufficient condition of invertibility of the asymptotic Fisher information matrix and Example 3 shows that the condition is not necessary.
We did not examine the complexity of our algorithm because the concept is not well-established for symbolic computation.Let us just say that, at this stage, it will take much more running time, as seen in the special case discussed in Section 3.For instance, the timings given for the ARMAX model of Section 3 show that the numerical method is considerably faster.Also, the program implementing the algorithm relies on the capabilities of the symbolic computation package and possible hardware and software limitations.This is the price to pay for the ability to provide exact results, at least when the matrix entries are given as fractions, not as decimal numbers.We believe, however, that the fact that, when it works at least, the notebook that we have produced the result for an arbitrary VARMAX model.
We should conclude with possible extensions of our results.There have been papers on the asymptotic information matrix for extended models like ARMA models with periodic coefficients [45,46], state-space models [46], Markov switching VARMA models [47], or VARFIMA models [48].It is possible that the approach described in this paper can be extended to these extended models.

Figure 1 .
Figure 1.Snapshot of the Mathematica notebook to treat a VARMAX model, given A[z], B[z], and C[z].The result Fcal is F (ϑ) as defined in (9).See the text for the other notations.

Figure 2 .
Figure 2. Partial snapshot of the Mathematica notebook to treat the example in(10).Gcal[z] is G(z) and Kcal[z] is G(z), both defined in(7).

Figure 3 .
Figure 3. Partial snapshot of the Mathematica notebook to treat the example in (10).Pcal[z] is P (z) and Qcal[z] Q(z), both defined in (8).To save space, only the first two elements of the first integrand in (9) are shown.

Figure 4 .
Figure 4. Partial snapshot of the Mathematica notebook to treat the example in(10).Fcal1 and Fcal2 are, respectively, the first and second terms in the right-hand side of (9) and Fcal is their sum, so equal to F (ϑ).

Figure 6 .
Figure 6.Output of the inverse of the Fisher information matrix for Example 3.

Figure 7 .
Figure 7. Output of the Fisher information matrix for Example 4.
.nb: Mathematica notebook containing Figures 1-4; File Figure 1234.pdf:pdf version of the Mathematica notebook containing Figures 1-4; File Example51.nb:Mathematica notebook corresponding to Example 5.1; File Example51.pdf:pdf version of the Mathematica notebook corresponding to Example 5.1.