Yule-Walker Estimation for the Moving-Average Model

The standard Yule-Walker equations, as they are known for an autoregression, are generalized to involve the moments of a moving-average process indexed on any number of dimensions. Once observations become available, new moments estimators are set to imitate the theoretical equations. These estimators are not only consistent but also asymptotically normal for any number of indexes. Their variance matrix resembles a standard result from maximum Gaussian likelihood estimation. A simulation study is added to conclude on their efficiency.


Introduction
At first, processes taking place on more than one dimension were considered on surfaces or in the three-dimensional space for the sake of spatial dynamics. Next, the time series Autoregressive Moving-Average model could be as well used to clothe the covariance dependence of spatial point processes. Nevertheless, the assumptions of causality and invertibility were based on conventional orderings of the vector indexes, as these were introduced by Whittle 1 and generalized by Guyon 2 , and they failed to look natural on the spatial axes. The more general bilateral equations considered by Whittle 1 also related to serious problems during the estimation of the parameters. Besag's 3 automodels as well as Besag's 4 equations for observations on irregular sites moved in a new direction but the weaknesses during the estimation remained.
The inclusion of the time axis has provided more convenience for the spatial analysis that may be summarized in two ways. On the one hand when the observations are collected on irregular sites, the asymptotic properties of the estimators for the spatial dependence parameters can now be established as the number of timings only increases to infinity; for example, Dimitriou-Fakalou 5 has proposed a methodology on how to proceed. On the other hand when the regular recordings are increasing towards all the spatial as well as the time axes, the unilateral ordering of indexes can now be arranged in a meaningful way and the causal and invertible ARMA model may be preferred as a natural spatiotemporal equation.
Special cases of the ARMA model are the pure autoregressive and moving-average equations and they both relate to useful properties as in the spectral density of interest. The autoregressive process uses a spectral density with "finite denominator", which for Gaussian processes 3 translates into a finite conditional expectation of any value based on all other values. A moving-average process uses a "finite numerator" in the spectral density, resulting in the autocovariance function with nonzero values on a finite set of "lags" only.
For a causal and finite autoregressive equation, the method of moments estimators according to the Yule-Walker equations almost coincides with the least squares or maximum Gaussian likelihood estimators and their consistency as well as their asymptotic normality may be established. Nevertheless, if the estimation for a pure moving-average or ARMA multi-indexed process takes place using its infinite autoregressive representation, a problem known as the "edge-effect" 2 resurrects and the asymptotic normality cannot be guaranteed. A modification of the Gaussian likelihood, which has been written from the autoregressive representation of each observation based on its "previous" ones, has been proposed by Yao and Brockwell 6 to confine the edge-effect with a weak condition of a finite second moment, but this works for the case of processes indexed on two dimensions only.
In this paper, for the multidimensional processes that may be clothed via an invertible and finite moving-average equation, a new method of estimation is proposed imitating the Yule-Walker equations. The finite number of nonzero autocovariances of the moving-average process is used most advantageously to eliminate the edge-effect. As a result, the proposed estimators are not only consistent but also asymptotically normal using relaxed conditions of finite second moments.

Theoretical Yule-Walker Equations
For Z {0, ±1, ±2, . . .} and d any positive integer, real-valued processes indexed on Z d will be considered, such as time series or spatial processes on one line transect d 1 , spatial processes in the plane d 2 or in the three-dimensional space, as well as the combination of spatiotemporal processes on Z 2 , Z 3 , or Z 4 , respectively. For any v ∈ Z d , an invertible movingaverage process {Y v } will be defined by the following equation: where {ε v } are uncorrelated, zero-mean random variables with positive and finite variance σ 2 . In 2.1 a unilateral ordering of the fixed lags 0 < i 1 < · · · < i q has been considered. For d 1, this is the standard ordering of integer numbers and, for d 2, the conventional unilateral ordering has been introduced by Whittle 1 . Guyon 2 has generalized the ordering of two vector indexes on Z d , d ≥ 3. As for the invertibility condition, it is only imposed to make sure that it can be written also using the ordering i > 0. Anderson and Jury 7 have provided a sufficient condition on the invertibility of such filters.
that is, the model is invertible. Note that, in either equation, no unnatural ordering is forced over space as the unilateral conventions "run" over the past timings only. Hence the model both uses a space-time unilateral ordering of indexes and is meaningful. As a result, for the special case of a stationary process with autocovariance function identified to take zero values after a reasonably small number of lags on all spatial and temporal dimensions, the invertible moving-average equation should be preferred and the new results presented next will be useful.

Method of Moments Estimators
For S ⊂ Z d being a set of finite cardinality N, there are available {Y v , v ∈ S} from 2.2 with true parameters θ 0 θ i 1 ,0 , . . . , θ i q ,0 τ that are to be estimated. The following condition will be assumed to be true.
Condition C1. The parameter space Θ ⊂ R q is a compact set containing the true value θ 0 as an inner point. Further, for any θ θ i 1 , . . . , θ i q τ ∈ Θ the moving-average model 2.2 is invertible.
Based on S, a set F v will be defined for any v ∈ Z d . More specifically, it will hold that i ∈ F v only if v − i ∈ S. Next, for the corrected set S * with N * elements, it will hold v ∈ S * if v i ∈ S for every i ∈ F. Thus for any v ∈ S * , it holds that F ⊆ F v . The reduction from S to S * is essential as it will guarantee later that the edge-effect will fade away; this is because the set S * ⊂ S includes the locations v that have all their "neighbors" As the source of the edge-effect is the speed of the absolute bias to zero, S * will guarantee that whatever "is missing" from the finite sample will not add at all to the bias of the estimators.
By imitating 2.11 the estimators θ θ i 1 , . . . , θ i q τ are set to be the solutions of equations: will be used additionally, in order to make sure that the estimators are consistent. The proposition guarantees that the use of q equations used as prototypes for estimation can provide a unique way to compute the q parameters of the invertible equation. International Journal of Stochastic Analysis and θ z −1 1 i>0 Θ i z i with i>0 |Θ i | < ∞, as well as γ z θ z θ z −1 ≡ j∈F γ j z j and c z γ z −1 ≡ j∈Z d c j z j , the unique solution that satisfies the q equations j∈F γ j,0 c j−i n 0, n 1, . . . , q, and j∈F γ j,0 c j− i n −i m 0, n, m 1, . . . , q, n > m, is θ i n ≡ θ i n ,0 , n 1, . . . , q.
For mathematical convenience, new variables depending on the sampling set are defined ∈ S, and 3.1 may be rewritten as where the zero subindex relates to θ 0 . In 3.3 and 3.4 , it is J n J n,1 , . . . , J n,q with elements

3.5
The last term in 3.5 will become more obvious when the properties of the estimators will be established but for now it might be justified by a second-order Taylor expansion, where it holds due to the fact that θ will be shown to be consistent estimators in Theorem 4.1. Further, the sums over S * involving the second derivatives of the Taylor expansion, when divided by N * or N N * /N → 1 , will also tend in probability to the relevant finite expectations when {ε v } will be assumed to be independent and identically distributed random variables with a finite variance. Finally by imitating 2.12 , the estimator of the error variance may be defined as Example 3.2. For the simplest one-dimensional case when q * q 1, it is demonstrated in detail how to compute the new estimator; some other cases will briefly be considered later.
Note that the "lags" on which the autocorrelations of Y are not equal to zero are as in the set F {1, 0, −1}; further, for any The autocovariance of the relevant auto-regression at the "lag" j ∈ Z is −θ |j| / 1−θ 2 . As a result, the estimator θ will be the solution of the following equation:

3.8
Furthermore, when {e v } are independent and identically distributed random variables, it will hold in probability that . . and the estimator may be approximated as a solution of the quadratic equation: Later, it will be investigated whether the approximation is worthwhile. The estimator reduces to

3.12
For the actual nonzero auto-correlation ρ θ/ 1 θ 2 , it holds that |ρ| < 1/2 and D 1−4ρρ > 0. As a result, if 0 < θ < 1 and 0 < ρ < 1/2, then 1/2ρ > 1, and the value International Journal of Stochastic Analysis is bigger than 1. If, on the other hand, −1 < θ < 0, and −1/2 < ρ < 0 then 1/2ρ < −1 and the same value is smaller than −1. Thus, the estimator is approximately 3.14 Although the distribution of θ is not known for small sample sizes and | ρ |, | ρ − | risk to be outside the −0.5, 0.5 range with some positive probability, both the consistency and the asymptotic normality will be established next. As opposed to what is known for the maximum Gaussian likelihood estimators, the proposed estimators will be asymptotically normal even for processes on more than one dimension. Thus, the new methods might be more useful for higher number of dimensions as well as large sample sizes.

Properties of Estimators
The asymptotic normality of standard estimators for the parameters of stationary processes indexed on Z d , d ≥ 2, has not been established yet. More specifically as Guyon 2 demonstrated, the bias of the genuine Gaussian likelihood estimators, computed from N regular recordings increasing to infinity on all dimensions and at equal speeds, is of order N −1/d . Nevertheless in order to secure the asymptotic normality, the absolute bias when multiplied by N 1/2 should tend to zero, which is only true for d 1 as in Brockwell and Davis 9 .
Regarding the proposed estimators of this paper, their bias will stem from 3.3 as in the expected value of v∈S * which expresses what is "missing" from the sample. Even when the bias is multiplied by N * 1/2 , the random variable N * −1/2 . , q has a zero mean due to the selection S * and the fact that the "finite" autocovariance function of a moving-average process is being used. Thus, the edge-effect will not be an obstacle to establish the asymptotic normality of the estimators.
A series of ARMA model-based arguments have been used before to deal with the edge-effect by Yao and Brockwell 6 for two-dimensional processes as well as a weak condition of finite second moments. Guyon 2 , on the other hand, used the periodogram as in the Gaussian likelihood of observations proposed by Whittle 1 and required that the fourth moments would be finite. Dahlhaus and Künsch 10 improved certain weaknesses of Guyon's 2 modification on the likelihood but paid the price of having to reduce the dimensionality to secure their results. The proposed estimators of this paper will also use the weak condition of finite second moments, which is a great advantage, since they will refer to moving-average processes on any number of dimensions.

then under Conditions C1 and C2(i) and as N → ∞, it holds that
Also the vector ξ ≡ η −i 1 , . . . , η −i q τ and the variance matrix

Empirical Results
In this section, an empirical comparison of the proposed estimators to maximum Gaussian likelihood estimators for moving-average processes is presented in detail. The theoretical foundations in favour of the new estimators have been provided already when their asymptotic normality on any number of dimensions has been established based on finite second moments only. As a result, the speed of the bias to zero is not expected to cause them to perform worse than maximum likelihood estimators, especially when the dimensionality is large.
On the other hand, Theorem 4.3 attributes to the new estimators the variance matrix W * −1 q when, according to Hannan 11 , efficient estimation results in W −1 q with W q ≡ Var ξ and the same notation as in Theorem 4.3. W q and W * q are defined similarly but they are not, in general, equal. It seems that as the number of moving-average parameters of the process increases, the two types elements get closer. Nevertheless, as the pure moving-average model will be preferred when it is parsimonious, a decision is made here whether the new estimators are efficient then.
The investigation has started with the one-dimensional case by generating {Y 1 , . . . , Y N } from the model with one parameter only. The moments estimator has been approximated as in the example earlier, while the minimizer of v∈S * e 2 * v with respect to the parameter θ has been considered to be the Gaussian likelihood estimator. The true values θ 0.5, −0.8 have been considered and the sample size has been set equal to N 30, 100, 200, 300. Very encouraging results for the efficiency of the proposed estimator are presented in Table 1 as even when the sample size is still small, extreme differences in the variances of the two types estimators cannot be detected. It is the bias of the moments estimator that seems to be the only reason why it might be outperformed by the likelihood estimator in terms of the Mean Squared Error. Nevertheless, the speed with which the bias tends to zero is much faster as one would expect and, eventually, the new estimator for θ 0.5 performs better altogether.
Next for the two-dimensional case, {Y u, v , u, v 1, . . . , n} have been generated from International Journal of Stochastic Analysis 11 with {e u, v , u, v ∈ Z} being independent, Gaussian random variables with zero mean and variance unity. Using similar arguments as in the example, the moments estimators of β and γ might be derived by solving simultaneously as well as the third equation and {e * u, v } are zero mean, uncorrelated random variables all with variance unity. Under the condition | β| | γ| < 1, a special form of "causality" X u, v ≡ ∞ j,k 0 Θ j,k e * u − j, v − k , ∞ j,k 0 | Θ j,k | < ∞ is implied and the standard Yule-Walker equations may be used to easily write the needed autocorrelations c, as functions of β and γ. Again, the maximum likelihood estimators have been approximated by the minimizers of Table 2 verifies once more the conclusions drawn from the one-dimensional case. It is safe to consider that the new estimators are efficient and this is very apparent in the case that the parameters are in absolute value further away from zero β 0.5 and γ 0.45 . The striking case with the small sample size of around 50 points observed in the plane reveals that the variances of the moments estimators may even happen to be smaller. On the other hand, the bias heavily affects the results for the MSE for small sample sizes. Nevertheless as the sample size increases, the absolute bias of the likelihood estimators does not seem to decrease at all versus the bias of the proposed estimators that speedily reaches the zero value. As a result, the new estimators eventually equalize the MSE performance of the standard estimators.

A.1. Proof of Proposition 3.1
Under the invertibility condition for the polynomials θ z , there is a one-to-one correspondence between the q coefficients θ i 1 , . . . , θ i q and the q autocorrelations It is true that the coefficients θ i n ,0 , n 1, . . . , q, generate the numbers c j,0 , j ∈ Z d , and they are a solution to the q equations of interest. If there is another solution, say θ i n ,1 , n 1, . . . , q generating c j,1 , j ∈ Z d , then it must hold that  − ρ j,1 c j−i n ,1 0, n 1, . . . , q, ,1 0, n,m 1, . . . , q, n >  where n,m 1 , is nonsingular and there is a unique solution to the equations above, that is, ρ j,1 ρ j,0 , j ∈ F.

A.2. Proof of Theorem 4.1
For n 1, . . . , q, 3.1 may be rewritten as Under the assumption that {ε v } are independent and identically distributed, it holds as , according to Proposition 6.3.10 or the Weak Law of Large Numbers and Proposition 7.3.5 of Brockwell and Davis 9 , extended to include the cases d ≥ 2. Then for n 1, . . . , q and the first of two terms in A.5 , it may be written as International Journal of Stochastic Analysis For the next term, it may be written as due to the Cauchy-Schwartz inequality and the independence of the random variables Y v , Y v − j , j / ∈ F. Now for any θ ∈ Θ, it holds that c j ≡ c j θ is the corresponding autocovariance function of a causal auto-regression. This guarantees that the autocovariance function decays at an exponential rate and there are constants C θ > 0 and α θ ∈ 0, 1 , Similarly for the estimator θ, it holds that with probability one and E c 2 After combining the results for the two terms in A.5 , it can be written that where the first term has been set equal to zero. Thus, j−i n ∈F c j γ j−i n ,0 P → 0 exactly like the theoretical analogue j−i n ∈F c j,0 γ j−i n ,0 0 implies. Since q instead of q equations have been used, there is a unique solution θ 0 , according to Proposition 3.1 and θ P → θ 0 as N → ∞ and C2 i holds. Finally, the consistency for θ implies, according to 3.7 , that σ 2 P −→ σ 2 j∈F c j,0 γ j,0 σ 2 since j∈F c j,0 γ j,0 1. A.12

A.3. Proof of Proposition 4.2
According to 3.5 , for the n, m th element of J/N, n, m 1, . . . , q, it suffices to look at where the last term tends to zero in probability, thanks to the consistency of the estimators from the use of all the q equations. For the polynomial the second term in A.1 can be written as This comes straight from the fact that For the auto-regression {X v } as it was defined in 2.3 , it can be seen immediately that Y v is uncorrelated to X v i , i > 0, since these are linear functions of "future" error terms only. In general, it can be written according to 2.5 that and back to the general Yule-Walker equations. Thus, International Journal of Stochastic Analysis and Y v is uncorrelated to X v − i , i / 0. As a result, it holds for n, m 1, . . . , q, that and, for n ≥ m, that The proof is completed when it is seen that both Y v and X v are linear functions of members from the sequence {ε v } and, thus, A.9

A.4. Proof of Theorem 4.3
From 3.3 , 4.2 , and 4.3 and for n 1, . . . , q, it can be written that The convergence in probability to zero of the remainder might be justified by the fact that its expected value is equal to zero due to the selection S * , and that its variance is given the sampling set S. First, {ε v } are independent and identically distributed and then without the assumption of a finite third or fourth moment. Under Condition C2 ii , it can be written that and a similar argument may be used for the cross-terms due to the Cauchy-Schwartz inequality. For the case d 2 and recordings on a rectangle, there is a justification for that by Yao and Brockwell 6 . Thus, Var N −1/2 v∈S * u n v → 0 as N → ∞ and Condition C2 holds, which guarantees the convergence in probability to zero. Equation 4.2 may now be rewritten as Then for n, m 1, . . . , q, it can be written for j ≥ 0 that Thus, for X v ≡ X v i 1 , . . . , X v i q τ , C j,0 ≡ E X v X v j τ /σ 2 , it can be written that E U v 0 and Now, for any positive integer K, the set is defined. According to the infinite moving-average representation of X v , for fixed K the new process {X K v } may be defined from International Journal of Stochastic Analysis Similarly, U K v τ ≡ X K v i 1 , . . . , X K v i q Y v , v ∈ Z d , and are defined here. For the same reasons as before, it holds that E U K v 0 and that For any vector λ ∈ R q , it holds that {λ τ U K v } is a strictly stationary and K * -dependent process, for a positive and finite integer K * . The definition of K-dependent processes as well as a theorem for the asymptotic normality for strictly stationary and K-dependent processes on Z d might be given similarly to the one-dimensional case by Brockwell and Davis 9 . Then γ j,0 C j,0 C τ j,0 , A.14 it holds as K → ∞ that λ τ M K λ → λ τ Mλ. Using Chebychev's inequality, it may be verified that as N → ∞ and under Condition C2. The proof will be completed when it is shown that Θ τ 0 Θ 0 W * q . Indeed, for the vector W W −i 1 , . . . , W −i q τ and for ξ Θ τ 0 W R with R being a q × 1 random vector that is independent of W, since it is a linear function of W −i 1 − i , i > 0, i / i 2 − i 1 , . . . , i q − i 1 , the required result might be obtained then, as A.20 Note that the decomposition W * q Θ τ 0 Θ 0 , as it is justified in the end of proof of Theorem 4.3, guarantees that W * q has a nonzero determinant.