Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 6s (2025) 340 https://internationalpubls.com Spatial Moving Average Polynomial Model on Irregular Lattice and its Full-Likelihood Based Implementation Lidia Khonglam1 Sanjeeva Kumar Jha2 1Department of Statistics, NEHU, Shillong- 793022, India. Email: lidiaklam123@gmail.com 2Department of Statistics, NEHU, Shillong- 793022, India. Email: skjha@nehu.ac.in Article History: Received: 20-10-2024 Revised: 05-12-2024 Accepted: 12-12-2024 Abstract: The simultaneous incorporation of a contiguity-based neighbourhood matrix of first and second order neighbors in the spatial moving average model on an irregular lattice can be specified in two different ways- as a quadratic polynomial or as a product of two linear polynomials. The present article concerns with the second order spatial moving average polynomial model of the latter form. Expression for the joint probability distribution of the response is derived together with expressions for the first and second order moments of the model using characteristic function. Model parameters are estimated using method of maximum likelihood and expression for their standard errors are obtained. In the absence of an analytical solution of the full likelihood function, optimization of the full likelihood function using Differential Evolution technique has been performed. Confidence interval for drawing inference has been constructed based on the derived standard error expression and bootstrapping method. Using likelihood ratio test statistic and bootstrap results, simultaneous confidence regions for the two spatial dependence parameters are also constructed. Implementation of the model has been executed on simulated data, and relevant comparisons made with linear models and first order spatial moving average model. The advantage of using the proposed model lies in its ability to provide statistically valid inference, as it affects producing spatial dependence of the first and the second lag orders in the residuals, insignificant. Keywords: spatial moving average; spatial lattice; contiguity neigbourhood; maximum likelihood; Differential Evolution; simultaneous confidence region; bootstrap 1. Introduction Economic, biological, health and ecological data on space commonly occur as aggregated data over geographically or administratively defined regions, forming regular or irregular spatial lattices. Spatial regression models, model for regression analysis of such spatially dependent data, incorporate spatial dependence in the concept of a contiguity neighborhood matrix constructed based on the shape of a spatial lattice. The prominent spatial regression models for aggregated data are the Conditional Autoregressive (CAR), Simultaneous Autoregressive (SAR) and the Spatial Moving Average (SMA) model. CAR was developed by Besag (1974) under the specification of a conditional probability on the response variable. The SAR and SMA model are specified by the distributional assumption on the error; hence provide joint probability of the response variable, Cressie (1993). SAR was developed by Whittle (1954) and SMA by Hainning (1978). SAR and CAR consider the observation of the response at the neighboring units of a particular region as an additional covariate and associated with it a spatial dependence parameter. The defined dependence is then incorporated into the mean structure of the linear models, similar to the time series mailto:lidiaklam123@gmail.com Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 6s (2025) 341 https://internationalpubls.com autoregressive model. SAR functional form and CAR model’s property, respectively analogue the autoregressive and Markov property of time series models (Wall, 2004; Cressie, 1993). SMA model on the other hand, resembled the time series moving average model. In this model, the disturbances at the neighboring units of a particular region are considered as an additional covariate with spatial dependence parameter associated with it. Implementation of CAR is more suitable with the Bayesian approach whereas SAR and SMA with likelihood based approaches. All three models have established its significance in enhancing the precision of the parameter estimates and capturing spatial dependence. Aggregated data on spatial lattices exhibit the existence of various kind of spatial interaction which causes spatial dependence. The geographical location at which an event occurs is the main factor of data on space; structuring dependence based on it forms the basis of spatial regression models. It assumed that an observation at one location shows an effects level that are similar to those of its neighbouring units, Anselin et al., 1998. The mentioned spatial regression models usually, established effects level based on the first order neighbours of the observation at each region. In essence, spatial dependence is not confined to the first order neighbours only. Somewise the higher order neighbours may or may not have an effect on the observation in a region. Analysis using the usual spatial regression models may have a substantial impact on the accuracy of the estimates, if higher order neighbours have a significant effect on the data at a region. Therefore, it is important to take into account the effects level based on higher order neighbourhood structure. Moreover, the incorporation of a higher order neighbourhood structure in the model allows the examination of complicated interactions that may exist and, to account for its presence, which may not be possible with just the first order neighbourhood framework. This manuscript study the specification of spatial dependence based on the first and second order neighbours in the SMA model and provides its implementation. The incorporation of a contiguity matrix of second order neighbors ( )2W into the SMA model with first order contiguity neighbors ( )1W can be specified analogous to the moving average time series models. It follows extending the linear polynomial ( )11WI + associated with the error structure of SMA to the form ( )2211 WWI  ++ , a quadratic polynomial as named by Lesage & Pace, 2009. Likelihood based implementation of the second order SMA specified hereby involves the evaluation of 2211 WWI  ++ which consist of both the first and second order spatial dependence parameters. The computation of 2211 WWI  ++ becomes complicated with higher order neighbourhood structure, as the contiguity matrix becomes less sparse with increasing order of the neigbourhood structure, Lesage and Pace, 2009. The author suggested factorizing ( )2211 WWI  ++ into a product of two linear polynomials i.e., ( )( )2211 WIWI  ++ so as to ease the computational complexity. Thus, the second order SMA polynomial model studied in this article is of the form where ( )( )2211 WIWI  ++ is associated with the error structure, as given in section 2. Theoretically, the specification provides plausible ranges of the two spatial dependence parameters. The second order SMA polynomial model proposed is implemented under the likelihood paradigm. The maximum likelihood (ML) estimate of the parameters is derived as shown in Section 3. It may Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 6s (2025) 342 https://internationalpubls.com be noticed that the analytical expression of the parameters ML estimates are functions of the other estimates. In the case of 2 , it is a function of  ˆ,ˆ,ˆ 21 and , and for  , it is a function of 21 ˆˆ  and . Expressing the estimates of 2 and  accordingly as a function of 21  and in the likelihood function, ML estimates of 21  and may be initially obtained. ML estimates of  and 2 may be correspondingly attained. Thus, estimation is based on modified likelihoods. In this paper estimation is not based on modified likelihoods. Differential evolution (DE) algorithm, proposed by Price & Storn (1997) available in R package; an optimization method for any non-differentiable or non-linear function is used as the method of optimization of the full likelihood function of the study model. The DE optimization method consists of the initialization, mutation, recombination and selection steps. The optimization of the full likelihood function provides desirable properties of its estimates relative to consistency, asymptotic efficiency and asymptotic normality, as compared to the other likelihood based estimation. It also made the construction of the simultaneous confidence regions for the two spatial dependence parameters possible. The optimization method has been successfully used for the implementation of skew Normal SAR model, Jha S.K. et al. (2021) and SMA model of first order, Jha S.K. et al. (2023). The main purpose of the paper is to study the implementation of the second order SMA polynomial model and to acknowledge its importance in modeling relationships in spatially dependent data. The remaining of the paper is organized as follows: Section 2 define the form of the second order SMA polynomial, Section 3 gives the ML derivation of the model, Section 4 provides the implementation of the model on simulated data and Section 5 presented the conclusion. 2. Models Let nsss .....,,, 21 be the n-sites in the study region S which forms a regular or irregular spatial lattice. Let iz and ( )Tikiii xxxx .....,,, 10= be respectively the response and explanatory variables observed at each of the site niSsi ,.....3,2,1; = . Haining (1978) defined Spatial Moving Average (SMA) model as: ( )εβZ 11WIX ++= (2.1) Where ( )IN 2,0~  ; ( )k ,....,, 10= is a vector of regression parameters, 1 is the unknown scalar spatial dependence parameter associated with the first order neighbors and 1W is a symmetric contiguity matrix of first order neighbors with elements 1 ijw defined as:         = = ji ji ij ssifandotherwise neighborsorderfirstaresandsif w ,0 1 1 (2.2) 2.1.Second order SMA Polynomial Under the same notations as mentioned above, the second order SMA polynomial may be defined as: Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 6s (2025) 343 https://internationalpubls.com ( )( )εβZ 2211 WIWIX  +++= (2.3) Where ( )IN 2,0~  ; 2 is the unknown scalar spatial dependence parameter associated with the second order neighbors and 2W is a symmetric contiguity neighborhood matrix of second order neighbors with elements 2 ijw defined as:         = = ji ji ij ssifandotherwise neighborsorderondaresandsif w ,0 sec1 2 (2.4) The probability density function of ε , takes the form, 𝑓(휀/0, 𝜎2𝐼) = 𝑒 − 𝑇 2𝜎2 𝜎(2𝜋)𝑛/2 = 𝑒 − [(𝐼+𝜃2𝑊2) −1 (𝐼+𝜃1𝑊1) −1 (𝑍−𝑋𝛽)] 𝑇 [(𝐼+𝜃2𝑊2) −1 (𝐼+𝜃1𝑊1) −1 (𝑍−𝑋𝛽)] 2𝜎2 𝜎(2𝜋)𝑛/2 The corresponding kernel is, 𝐾(𝑍) = 𝑒 − [(𝐼+𝜃2𝑊2)−1(𝐼+𝜃1𝑊1)−1(𝑍−𝑋𝛽)] 𝑇 [(𝐼+𝜃2𝑊2)−1(𝐼+𝜃1𝑊1)−1(𝑍−𝑋𝛽)] 2𝜎2 Normalizing constant can be obtained as; ∫ 𝐾(𝑍)𝑑𝑍 = ∫ 𝑒 − [(𝐼+𝜃2𝑊2)−1(𝐼+𝜃1𝑊1)−1(𝑍−𝑋𝛽)] 𝑇 [(𝐼+𝜃2𝑊2)−1(𝐼+𝜃1𝑊1)−1(𝑍−𝑋𝛽)] 2𝜎2 Let, (𝐼 + 𝜃2𝑊2)−1(𝐼 + 𝜃1𝑊1)−1(𝑍 − 𝑋𝛽)𝜎−1 = 𝑈 ∴ 𝑍 = (𝐼 + 𝜃1𝑊1)(𝐼 + 𝜃2𝑊2)𝜎𝑈 + 𝑋𝛽 The differential of Z with respect to U is 𝜎|(𝐼 + 𝜃1𝑊1)(𝐼 + 𝜃2𝑊2)| ∴ ∫ 𝐾(𝑍)𝑑𝑍 = ∫ 𝑒 − 𝑈𝑇𝑈 2 (2𝜋)𝑛/2 × 𝜎|(𝐼 + 𝜃1𝑊1)(𝐼 + 𝜃2𝑊2)|(2𝜋)𝑛/2𝑑𝑈 = 𝜎(2𝜋)𝑛/2|(𝐼 + 𝜃1𝑊1)(𝐼 + 𝜃2𝑊2)| (2.5) Equation (2.5) gives the normalizing constant for the joint distribution of Z as, 𝜎(2𝜋)𝑛/2|(𝐼 + 𝜃1𝑊1) (𝐼 + 𝜃2𝑊2)|. If ( )11WI + and ( )22WI + is a full rank matrix then the induced joint density function of Z takes the form; 𝑓(𝑍) = 1 𝜎(2𝜋)𝑛/2|(𝐼 + 𝜃1𝑊1)(𝐼 + 𝜃2𝑊2)| × exp {− [(𝐼+𝜃2𝑊2)−1(𝐼+𝜃1𝑊1)−1(𝑍−𝑋𝛽)] 𝑇 [(𝐼+𝜃2𝑊2)−1(𝐼+𝜃1𝑊1)−1(𝑍−𝑋𝛽)] 2𝜎2 } (2.6) Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 6s (2025) 344 https://internationalpubls.com Equation (2.6) gives the required joint probability density function under the proposed second order SMA polynomial. Let the distribution be denoted henceforth as; 𝑍~2 𝑛𝑑 − 𝑜𝑟𝑑𝑒𝑟 − 𝑆𝑀𝐴 − 𝑝𝑜𝑙𝑦𝑛𝑜𝑚𝑖𝑎𝑙 (𝛽, 𝜎2, 𝜃1, 𝜃2) 2.2.Characteristic function and moments: The derivation of the characteristic function and the corresponding moments of the response vector ‘Z’ is given in Appendix A. The obtained characteristic function is: Ψ𝑍(𝑡) = 𝑒𝑖𝑡𝑇𝑋𝛽− 𝑡𝑇𝜔𝑡 2 (2.7) Where 𝜔 = 𝜎2(𝐼 + 𝜃1𝑊1)(𝐼 + 𝜃2𝑊2)(𝐼 + 𝜃2𝑊2)𝑇(𝐼 + 𝜃1𝑊1)𝑇 And the derived mean vector and variance-covariance matrix of ‘Z’ is: 𝐸(𝑍) = 𝑋𝛽 and, 𝐶𝑜𝑣 (𝑍) = 𝜎2(𝐼 + 𝜃1𝑊1)(𝐼 + 𝜃2𝑊2)(𝐼 + 𝜃2𝑊2)𝑇(𝐼 + 𝜃1𝑊1)𝑇 3. Likelihood function and Maximum Likelihood (ML) estimate The log- likelihood function of second order SMA polynomial is: ( )( ) ( ) ( ) ( )( ) ( ) ( ) ( )  ( ) ( ) ( )  )1.3( 2 1 loglog 2 2log 2 /,,,log 1 11 1 22 1 11 1 222 2211 2 21 2    XZWIWIXZWIWI WIWI nn zL T −++−++− ++−−−= −−−− Let ( )( ) 212211 GGWIWIG =++=  where ( ) 222111 WIGandWIG  +=+= Then, the first partial derivative of (3.1) w.r.t the unknown parameters 21 2 ,,,  respectively is, ( ) ( ) ( )  XZGGX L TT −=   −− 11 2 1log (3.2) ( ) ( )  ( )   XZGXZG nL T −−+−=   −− 11 422 2 1 2 log (3.3) ( ) ( ) ( ) ( ) ( ) ( ) ( ) )4.3( 2 1log 1 21 111 21 11 221 1 1   XZGGWGGGGWGGXZGWGtr L TTTT −+−+−=   −−−−−−− And, ( ) ( ) ( ) ( ) ( ) ( ) ( ) )5.3( 2 1log 1 12 111 12 11 212 1 2   XZGGWGGGGWGGXZGWGtr L TTTT −+−+−=   −−−−−−− At the maximum of 21 2 ,,,  respectively, the expected values of partial derivatives (3.2) to (3.4) must be zero. The proofs for the first two are familiar; therefore we will show here only for the two spatial dependence parameters. Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 6s (2025) 345 https://internationalpubls.com ( ) ( ) ( )  ( ) ( ) ( )( )( ) ( ) ( ) 0 2 2 1log 21 1 21 1 2 21 12 21 12 21 1 21 1 21 1 221 1 1 = +−= + +−= ++−=        −− −− − −−− GWGtrGWGtr GWGtrGWGtr GWGtr GWGGWGEGWGtr L E T TTT     Similarly, ( ) ( ) ( )  ( ) ( ) 0 2 1log 12 1 12 1 12 1 12 1 212 1 2 = +−= ++−=        −− −−− GWGtrGWGtr GWGGWGEGWGtr L E TTT   The elements of the information matrix are:    XGXG L E T T 11 2 2 1log −−=        −  (3.6) 0 log 2 2 =        −  L E (3.7) 0 log 1 2 =        −  L E (3.8) 0 log 2 2 =        −  L E (3.9) 422 2 2 log  nL E =        − (3.10) ( ) 2 1 21 1 2 2 log  − =        − GGWtrL E (3.11) ( ) 2 1 22 2 2 2 log  − =        − GGWtrL E (3.12) ( ) ( )( ) ( )( ) 1 21 1 21 21 212 1 2 log −−− +=        − GGWGGWtrGGWtr L E T  (3.13) Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 6s (2025) 346 https://internationalpubls.com ( ) ( )( ) ( )( ) 1 12 1 12 21 122 2 2 log −−− +=        − GGWGGWtrGGWtr L E T  (3.14) ( )21 21 2 2 log WWtr L E =        −  (3.15) The information matrix ( )21 2 ,,, I may then be expressed as:     ( ) ( ) ( ) ( ) ( )( ) ( )( )  ( ) ( ) ( ) ( ) ( )( ) ( )( )  )16.3( 20 20 2 0 000 1 1 12 1 12 21 12212 1 12 21 1 21 1 21 21 212 1 21 2 1 12 2 1 21 4 11 2                       + + −−− − −−− − −− −− GGWGGWtrGGWtrWWtr SGWtr WWtrGGWGGWtrGGWtr GGWtr GGWtrGGWtrn XGXG T T T     The information matrix (3.16) has block-diagonality between group of parameters, similar to the first order SMA model as shown in Hepple (2003); therefore, the variance for ̂ may be calculated using, 𝑉𝑎𝑟 (�̂�) = 𝜎2[((𝐼 + 𝜃2𝑊2)−1(𝐼 + 𝜃1𝑊1)−1𝑋)𝑇((𝐼 + 𝜃2𝑊2)−1(𝐼 + 𝜃1𝑊1)−1𝑋)]−1 (3.17) With the ML estimate of 1 and 2 . Further, the sub-matrix involving only 21 2 ,,  can be used to obtain the variance of 21  and . At the maximum, the partial derivatives of equation (3.2) and (3.3) equals zero. Therefore, �̂� = [𝑋𝑇𝐺−1𝑇 𝐺−1𝑋] −1 𝑋𝑇𝐺−1𝑇 𝐺−1𝑍 (3.18) And, �̂�2 = 𝑛−1𝐷∗𝑇 𝐺−1𝑇 𝐺−1𝐷∗ where, 𝐷∗ = 𝑍 − 𝑋�̂� (3.19) Under the specified values of 21  and , the estimates of  can be obtained as a linear model, using the Generalised Least Squares (GLS) equation. The log likelihood function in equation (3.1) after replacing 𝛽 by �̂� of (3.18) and 2 by 2̂ of (3.19) becomes, log(𝐿(𝛽, 𝜃1, 𝜃2/𝑍)) = − 𝑛 2 log(2𝜋) − 𝑛 2 log (𝑛−1𝐷∗𝑇 𝐺−1𝑇 𝐺−1𝐷∗) − log|𝐺| − 1 2 = − 𝑛 2 [log ((2𝜋)𝑛−1) + 𝑛−1] − 𝑛 2 log [|𝐺|2/𝑛𝐷∗𝑇 𝐺−1𝑇 𝐺−1𝐷∗] (3.20) Thus, numerical optimization of (3.20) gives estimates of 21  and , after which 2ˆˆ  and could be obtained using (3.18) and (3.19) respectively. The present paper however, uses the technique of DE algorithm for optimization of the full likelihood in (3.1). The method does not necessarily require the step wise derivation shown for obtaining the ML estimates of the parameters. It simultaneously estimates all the parameters; hence allows one to construct confidence intervals based on bootstrap method and for construction of simultaneous confidence region for the two or more spatial dependence parameters. Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 6s (2025) 347 https://internationalpubls.com 3.1. Parameter spaces for 𝜽𝟏 and 𝜽𝟐: Hepple, 2003 mentioned that the structure of 1W and its eigen values determine the range of 1 in the first order SMA. The condition arises from 11WI + that appears in the likelihood function of the model, as it should be positive (Anselin, 1988; Mur et al, 2006). Its range should be between the negative of the reciprocal of the largest and the smallest eigen values of 1W . Proceeded in a similar manner, the feasible parameter spaces for the two spatial dependence parameters of the second order SMA polynomial may be obtained theoretically. The jacobian of the model under study is: |𝐺| = |(𝐼 + 𝜃1𝑊1)(𝐼 + 𝜃2𝑊2)| = |𝐼 + 𝜃1𝑊1||𝐼 + 𝜃2𝑊2| = ∏ (1 + 𝜃1𝜆𝑖) 𝑅 𝑖=1 ∏ (1 + 𝜃2𝛿𝑖) 𝑅 𝑖=1 > 0 Thus, ∏ (1 + 𝜃1𝜆𝑖)𝑅 𝑖=1 > 0 and, ∏ (1 + 𝜃2𝛿𝑗)𝑅 𝑗=1 > 0 ⇒ 1 𝜆𝑚𝑎𝑥 < 𝜃1 < 1 𝜆𝑚𝑖𝑛 and, 1 𝛿𝑚𝑎𝑥 < 𝜃2 < 1 𝛿𝑚𝑖𝑛 Where max and min correspondingly represents the maximum and minimum eigen values associated with the weighted contiguity neighbors of first orders; max and min respectively are the maximum and minimum eigen values of the contiguity neighbors of second orders. In the present analysis, the range of 1 is (-1, 1.6108) and for 2 is (-1, 1.9195). 4. Implementation Simulated data is used for implementation of the study second order SMA polynomial. The New York shape file, which has 281 sub-regions available in ‘spData’ R-package, is used for constructing the spatial contiguity matrix of first and second order neighborhoods required for implementation. A vector of random error ( )25,0~ N and a vector of explanatory variable, ( ))30,20(~2 uniformX each having dimensions 281 are generated. 1X is a unitary vector with 281 dimension, such that  21 XXX = . 21 WandW are 281281 spatial contiguity matrix of first and second order neighborhoods respectively. Then, with 10,20 21 ==  , 1.0,3.0 21 ==  a response variable Z with 281 dimensions is obtained using (2.3). With the response variable, Z and the explanatory variable 2X simulated as above, analysis under the linear model assumption signifies the presence of a significant spatial dependence of the first two lag orders in the residuals. The Moran’s I value for the first lag order is 0.55045 and that of second lag order is 0.17189 with both p-value less than 0.05. Therefore, this simulated data is considered fit for implementation of second order SMA polynomial. Optimization of the full-likelihood function is carried out using the ‘DEoptim’ technique for obtaining the required parameter ML estimates. In running the program, convergence of the optimized log likelihood function and its corresponding parameters with respect to the number of iterations is firstly assessed, and the results are presented in Fig. 1 and Fig. 2 accordingly. Both Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 6s (2025) 348 https://internationalpubls.com figures show that the overall convergence is at about 70 iterations. The current analysis reports result of the parameter estimated in running the program for 200 iterations. Results from second order SMA polynomial: The obtained ML estimate of the parameters, its standard error and 95 per cent exact confidence interval is presented in Table 1. The variance for constructing the 95 per cent exact confidence interval for the covariates is obtained using (3.16) and, for the two spatial dependence parameters using the inverse of the information sub-matrix involving 21 2 ,  and . The table also, presents bootstrap results based on 1000 resampled values, and its respective density plots are presented in Fig 3. From table 1, it is observed that both the regression parameters 10  and are significant at 0.05 levels of significance. The estimated spatial dependence is 1.07464 and 0.76755 for the first and second lag orders respectively; both being significant at 5 per cent significance level. Bootstrap confidence interval also gives similar results for the regression parameters and the two spatial dependences. The estimated global variance under bootstrapping is 21.88125 and under the ML estimate is 25.05253. The observed response and the estimated response values under the SMA polynomial of second order are further presented in a choropleth map as shown in Fig 4. It may be noticed from the figure that second order SMA polynomial properly capture the variability in the dataset as ranges of the response values is reduced by 15 per cent from the observed values. The observed response ranges from 212 to 329 whereas the estimated response ranges between 219 and 320. Fig 1: Convergence plot of the optimized log-likelihood function Table 1: Parameter estimates, its standard error and 95 per cent confidence interval Estimates (se) Exact 95 per cent CI Bootstrap Mean se 95 per cent CI 0 18.708 (2.174) (14.447, 22.969) 18.348 3.695 (11.203, 25.689) 1 10.022 (0.008) (10.007, 10.037) 10.039 0.126 (9.795, 10.290) 1 1.075 (0.009) (1.057, 1.093) 1.143 0.111 (0.913, 1.352) 2 0.768 (0.116) (0.539, 0.996) 0.869 0.222 (0.401, 1.268)  5.005 4.678 AIC 1646.873 0 50 100 150 200 -8 23 -8 22 -8 21 -8 20 -8 19 -8 18 iteration lo g lik el ih oo d End value= -818.4365 Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 6s (2025) 349 https://internationalpubls.com *se stands for standard error and CI for confidence interval Simultaneous Confidence region between 21  and : The simultaneous confidence region for the two spatial dependence parameters based on likelihood ratio (LR) test statistic and bootstrap results is constructed. The simultaneous confidence region based on LR test is obtained by firstly generating 21  and within (-1, 1.611) and (-1, 1.919) respectively. The likelihood values with the different combinations of 21  and is calculated; then only those combinations of 21  and which gives the likelihood value greater than or equal to the maximized log likelihood values at 95 per cent confidence level is taken. Hence, the LR based simultaneous confidence region is as presented in Fig 5 (a). Bootstrap based confidence region is constructed by taking the lower 95 per cent values of the Mahalanobis distance between the bootstrap vector values of 21  and and the region is as presented in Fig 5(b). Fig 2: Convergence plot of the parameter estimate From both figures in Fig 5, it may be noticed that there is a positive correlation between 21  and . The calculated correlation coefficient between them is 0.249 for LR test and 0.376 for bootstrap values with both p-values being less than the 0.05 significance level. The 95 per cent confidence interval based on LR statistic is (0.929, 1.184) for 1 and for 2 is (0.426, 1.049). Based on bootstrap values, the 95 per cent confidence interval is (0.913, 1.352) and (0.401, 1.268) for 21  and respectively. 0 50 100 150 200 15 17 19 iteration 0 0 50 100 150 200 10 .0 2 10 .0 8 iteration 1 0 50 100 150 200 -2 00 40 0 iteration 0 50 100 150 200 1. 05 1. 15 iteration 1 0 50 100 150 200 0. 6 0. 8 iteration 2 Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 6s (2025) 350 https://internationalpubls.com Fig 3: Bootstrap result density plots (a) (b) Fig 4: (a) Observed response, and (b) estimated response under second order SMA polynomial 5 15 25 0.0 0 0.0 4 0.0 8 (a) De ns ity 0 9.6 10.0 10.4 0.0 1.0 2.0 3.0 (b) De ns ity 1 0.8 1.2 1.6 0.0 1.0 2.0 3.0 (c) De ns ity 1 0.0 0.5 1.0 1.5 0.0 0.5 1.0 1.5 (d) De ns ity 2 3.0 4.0 5.0 6.0 0.0 0.2 0.4 0.6 0.8 1.0 (e) De ns ity Observed (212,241] (241,261] (261,282] (282,303] (303,329] Estimated (220,242] (242,260] (260,280] (280,302] (302,319] Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 6s (2025) 351 https://internationalpubls.com (a) (b) Fig 5: Simultaneous confidence region between 21  and based on LR test and Bootstrap Table 2: Parameter estimates with its 95 per cent confidence interval under the linear model, first order SMA and second order SMA polynomial Linear models First order SMA Second order SMA polynomial 0 17.579 (11.139, 24.019) 19.205 (14.891, 23.518) 18.7083 (14.447, 22.969 ) 1 10.074 (9.819, 10.328) 10.00709 (9.992, 10.023) 10.02187 (10.007, 10.037) 1 ____ 1.007568 (0.959, 1.056) 1.075 (1.057, 1.092) 2 _____ _____ 0.76755 (0.539, 0.996) 2 39.163 25.417 25.053 AIC 1832.099 1667.936 1646.873 Table 3: Moran’s I statistic of the residual with its p-values under the Linear model, first order SMA and second order SMA polynomial Moran’s I statistic of the residuals Linear model First order SMA Second order SMA polynomial First lag order (p-value) 0.550 (<0.05) 0.103 (0.004) 0.067 (0.102) Second lag order (p-value) 0.172 (<0.05) 0.107 (<0.05) -0.002 (0.956) 0.0 0.5 1.0 1.5 0. 0 0. 5 1. 0 1. 5 LR based Confidence Region 1 2 0.0 0.5 1.0 1.5 0. 0 0. 5 1. 0 1. 5 Bootstrap based Confidence Region 1 2 Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 6s (2025) 352 https://internationalpubls.com Fig 6: Residuals plot under linear models, first order SMA, and SMA polynomial of second order Comparison of results: The simulated data generated is further analyzed under the linear models and first order SMA, and its results are presented in Table 2. It is observed that all three models give similar results in term of the regression parameters, both being significant at 95 per cent confidence interval. The overall gain in using the second order SMA polynomial rather than the first order SMA model and linear model is reflected in the AIC value. The AIC value for the proposed model is 1646.873 whereas for the linear model is 1832.099, and 1667.936 for SMA model of first order. The estimated variance is reduced by only 1 per cent when using the second order SMA polynomial. (a) (b) (c) Fig 7: Moran’s I statistic with its 95 per cent Confidence interval of the residuals under the (a) linear model, (b) first order SMA, (c) Second order SMA polynomial -20 -10 0 10 20 0. 00 0. 02 0. 04 0. 06 0. 08 , D en si ty 2nd order polynomial SMA 1st order SMA linear model 0. 2 0. 3 0. 4 0. 5 0. 6 Linear model lags M or an 's I 1 2 -0 .0 5 0. 00 0. 05 0. 10 0. 15 0. 20 First order SMA lags M or an 's I 1 2 -0 .0 5 0. 00 0. 05 0. 10 0. 15 0. 20 Second order polynomial SMA lags M or an 's I 1 2 Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 6s (2025) 353 https://internationalpubls.com Residual density plot under the second order SMA polynomial overlaid with those of the linear model and first order SMA is given in Fig 6. The first order SMA and second order SMA polynomial shows nearly same pattern in achieving smoothness as compared to linear model. The Shapiro-wilks normality test statistic provides a p-value of 0.073, 0.489 and 0.805 for the residuals under the linear model, first order SMA and second order SMA polynomial model respectively. Table 3 provides results of the analytical test of the presence of a significant spatial dependence of the first and second lag orders in the residuals. Correspondingly, Fig 7 provides a diagrammatical representation of the result with its 95 per cent confidence interval. From both, the table and figure we may be clearly observed that the residuals from the linear model and first order SMA exhibit a significant spatial dependence of both lag orders, whereas those from the second order SMA polynomial do not. 5. Conclusion The paper studies the aspect of modeling relationships between spatially dependent data through the simultaneous incorporation of the first and second order contiguity neighbourhood matrix in SMA, with the model named as second order SMA polynomial. The model provides a way of modeling, in cases where spatial dependence could not be properly explained in a single neighbourhood structure framework. The main achievement of the paper involves the successful implementation of the second order SMA polynomial through the full likelihood optimization using the DE method; construction of confidence intervals for statistical inferences based on the derived standard error expression and bootstrapping; and, the construction of the simultaneous confidence regions for the two spatial dependence parameters based on LR test and bootstrapping. The 95 percent bootstrap based confidence intervals are wider as compared to the exact 95 percent confidence intervals. In the present analysis, the incorporation of the second order spatial dependence parameters reduces the error variance by a small amount only, as compared to the first order SMA model. Hence, not so much change could be observed in the parameters confidence intervals of both models. Nevertheless, the efficiency of the model is reflected in its AIC values. The potential of the second order SMA polynomial to apprehend spatial dependence of both first and second lag orders, validates the precision of the inferences drawn about the data using the second order SMA polynomial. Acknowledgement: We would like to thank the R Core Team (2022). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. The URL address is https://www.R-project.org. We extend our sincere gratitude to Bivand et al. (2018), Bivand et al. (2021), and Mullen et al. (2011) for developing the useful R packages ‘spdep’, ‘spData’, and ‘DEoptim’, respectively. Conflict of interest:There is no conflict of interest among the authors. Appendix A If 𝑍~2 𝑛𝑑 − 𝑜𝑟𝑑𝑒𝑟 − 𝑆𝑀𝐴 − 𝑝𝑜𝑙𝑦𝑛𝑜𝑚𝑖𝑎𝑙 (𝛽, 𝜎2, 𝜃1, 𝜃2), then its characteristic function is, 𝜓𝑍(𝑡) = 𝐸(𝑒𝑖𝑡𝑇𝑍) https://www.r-project.org/ Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 6s (2025) 354 https://internationalpubls.com = 𝑒𝑖𝑡𝑇𝑍𝑒 − [(𝐼+𝜃2𝑊2)−1(𝐼+𝜃1𝑊1)−1(𝑍−𝑋𝛽)] 𝑇 [(𝐼+𝜃2𝑊2)−1(𝐼+𝜃1𝑊1)−1(𝑍−𝑋𝛽)] 2𝜎2 𝜎(2𝜋)𝑛/2|(𝐼 + 𝜃1𝑊1)(𝐼 + 𝜃2𝑊2)| 𝑑𝑍 Letting, 𝑢 = (𝐼 + 𝜃2𝑊2)−1(𝐼 + 𝜃1𝑊1)−1(𝑍 − 𝑋𝛽)𝜎−1 ⇒ 𝑍 = 𝜎(𝐼 + 𝜃1𝑊1)(𝐼 + 𝜃1𝑊1) + 𝑋𝛽 And, 𝐽 = 𝑑𝑍 𝑑𝑢 = 𝜎|(𝐼 + 𝜃1𝑊1)(𝐼 + 𝜃2𝑊2)| ∴ 𝜓𝑍(𝑡) = ∫ 𝑒𝑖𝑡𝑇𝑋𝛽 (2𝜋)𝑛/2 × 𝑒𝑖𝑡𝑇𝜎(𝐼+𝜃1𝑊1)(𝐼+𝜃1𝑊1)𝑢− 𝑢𝑇𝑢 2 𝑑𝑢 = 𝑒𝑖𝑡𝑇𝑋𝛽 (2𝜋)𝑛/2 ∫ 𝑒 [𝑢𝑇−2𝑖𝑡𝑇𝜎(𝐼+𝜃1𝑊1)(𝐼+𝜃2𝑊2)]𝑢 2 𝑑𝑢 Again letting, 𝑦𝑇 = 𝑢𝑇 − 𝑖𝑡𝑇𝜎(𝐼 + 𝜃1𝑊1)(𝐼 + 𝜃2𝑊2) ∴ 𝑢𝑇 = 𝑦𝑇 + 𝑖𝑡𝑇𝜎(𝐼 + 𝜃1𝑊1)(𝐼 + 𝜃2𝑊2) And, 𝑢 = 𝑦 + 𝑖𝜎(𝐼 + 𝜃2𝑊2)𝑇(𝐼 + 𝜃1𝑊1)𝑇𝑡 The differential of ‘u’ with respect to ‘y’ is 1. 𝜓𝑍(𝑡) = 𝑒𝑖𝑡𝑇𝑋𝛽 (2𝜋)𝑛/2 ∫ 𝑒− 1 2 [𝑦𝑇−𝑖𝑡𝑇𝜎(𝐼+𝜃1𝑊1)(𝐼+𝜃2𝑊2)][𝑦−𝑖𝜎(𝐼+𝜃1𝑊1)(𝐼+𝜃2𝑊2)𝑡] 𝑑𝑦 = 𝑒𝑖𝑡𝑇𝑋𝛽 (2𝜋)𝑛/2 ∫ 𝑒− 𝑦𝑇𝑦 2 𝑒 𝑖2𝜎2𝑡𝑇(𝐼+𝜃1𝑊1)(𝐼+𝜃2𝑊2)(𝐼+𝜃2𝑊2) 𝑇 (𝐼+𝜃1𝑊1) 𝑇 𝑡 2 𝑑𝑦 = 𝑒𝑖𝑡𝑇𝑋𝛽+ 𝑖2𝜎2𝑡𝑇(𝐼+𝜃1𝑊1)(𝐼+𝜃2𝑊2)(𝐼+𝜃2𝑊2) 𝑇 (𝐼+𝜃1𝑊1) 𝑇 𝑡 2 ∫ 𝑒 − 𝑦𝑇𝑦 2 (2𝜋)𝑛/2 𝑑𝑦 = 𝑒𝑖𝑡𝑇𝑋𝛽+ 𝑖2𝜎2𝑡𝑇(𝐼+𝜃1𝑊1)(𝐼+𝜃2𝑊2)(𝐼+𝜃2𝑊2) 𝑇 (𝐼+𝜃1𝑊1) 𝑇 𝑡 2 × 1 = 𝑒𝑖𝑡𝑇𝑋𝛽+ 𝑖2𝜎2𝑡𝑇(𝐼+𝜃1𝑊1)(𝐼+𝜃2𝑊2)(𝐼+𝜃2𝑊2) 𝑇 (𝐼+𝜃1𝑊1) 𝑇 𝑡 2 Let, 𝜔 = 𝑡𝑇(𝐼 + 𝜃1𝑊1)(𝐼 + 𝜃2𝑊2)(𝐼 + 𝜃2𝑊2)𝑇(𝐼 + 𝜃1𝑊1)𝑇𝑡 ∴ 𝜓𝑍(𝑡) = 𝑒𝑖𝑡𝑇𝑋𝛽− 𝑡𝑇𝜔𝑡 2 Now, 𝜓𝑍 ′ (𝑡) = 𝑒𝑖𝑡𝑇𝑋𝛽− 𝑡𝑇𝜔𝑡 2 𝑖𝑋𝛽 − 𝑒𝑖𝑡𝑇𝑋𝛽− 𝑡𝑇𝜔𝑡 2 𝜔𝑡 𝜓𝑍 ′′ (𝑡) = 𝑒𝑖𝑡𝑇𝑋𝛽− 𝑡𝑇𝜔𝑡 2 𝑖𝑋𝛽(𝑖𝑋𝛽 − 𝜔𝑡) − 𝜔𝑒𝑖𝑡𝑇𝑋𝛽− 𝑡𝑇𝜔𝑡 2 − 𝜔𝑡(𝑖𝑋𝛽 − 𝜔𝑡)𝑒𝑖𝑡𝑇𝑋𝛽− 𝑡𝑇𝜔𝑡 2 Therefore, 𝐸(𝑍) = (−𝑖)𝜓𝑍 ′ (𝑡)| 𝑡=0 Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 6s (2025) 355 https://internationalpubls.com = (−𝑖)𝑖𝑋𝛽 = −(−1)𝑋𝛽 =𝑋𝛽 𝐸(𝑍2) = (−𝑖)2𝜓𝑍 ′′ (𝑡)| 𝑡=0 = (−1)(𝑖2𝑋𝛽𝛽𝑇𝑋𝑇 − 𝜔) = 𝑋𝛽𝛽𝑇𝑋𝑇 + 𝜔 Thus, 𝐶𝑜𝑣 (𝑍) = 𝐸(𝑍2) − [𝐸(𝑍)]2 = 𝑋𝛽𝛽𝑇𝑋𝑇 + 𝜔 − 𝑋𝛽𝛽𝑇𝑋𝑇 = 𝜔 = 𝜎2(𝐼 + 𝜃1𝑊1)(𝐼 + 𝜃2𝑊2)(𝐼 + 𝜃2𝑊2)𝑇(𝐼 + 𝜃1𝑊1)𝑇 Refrences [1] Anselin, L. Bera, A.K. (1998).Spatial dependence in linear regression models with an introduction to spatial econometrics. Statistics: Textbooks and Monographs. 155, 237-289. [2] Anselin, L. (1988). Spatial econometrics: methods and models. Springer Science & Business Media. [3] Besag, J.(1974). Spatial interaction and the statistical analysis of lattice systems. Journal of the Royal Statistical Society. Series (Methodological). 36, 192-236. [4] Bivand, R. S. Wong, D. W.S. (2018). Comparing implementations of global and local indicators of spatial association test. TEST . 27( 3) , 716-748. https://doi.org/10.1007/s11749-018-0599-x. [5] Bivand R, Nowosad J, Lovelace R. (2021) spData:Datasets for Spatial Analysis. R package version 0.3.10. 2021. https://CRAN.R-project.org/package=spData. [6] Cressie, N. (1993). Statistics for spatial data. New York: John Wiley and sons. [7] Hainning, R.P. (1978). The moving average model for spatial interaction. Transactions of the institute of British geographers. 3( 2), 202-225. [8] Hepple, L.W. (2003). Bayesian and maximum likelihood of the linear model with spatial moving average disturbances. http://www.ggy.bris.ac.uk/personal/LesHepple/bayesian.pdf. [9] Jha, S.K. Begum, R. (2023). A proposed weighting scheme for spatial moving average model in irregular lattice. Journal of the Indian Society for Probability and Statistics. 24, 357-375. [10] Jha, S.K. Singh, N.V. (2021). A skew-normal spatial simultaneous autoregressive model and its implementation. Sankhya A. 85(1), 306-323. https://doi.org/10.1007/s13171-021-00246-3. [11] Lesage and Pace. (2009). Introduction to spatial econometrics. Statistics: Textbooks and Monographs. https://www.researchgate.net/publication/339938974. [12] Mur, J. Angulo, A.M. (2006). Clues for discriminating between moving average and autoregressive process in spatial processes. Verlag Springer. 9, 273-298. [13] Mullen, K. Ardia, D. Gil, D. Windover, D and Cline, J. (2011). 'DEoptim': An R Package for Global Optimization by Differential Evolution. Journal of Statistical Software. 40(6), 1-26. https://doi.org/10.18637/jss.v040.i06. [14] Storn, R. Price, K. (1997). Differential Evolution- A simple and efficient Heuristic for global optimization over continuous spaces. Journal of Global optimization. 11, 341-359. [15] Whittle P. (1954). On stationary process in the plane. Biometrika. 41, 434-449. [16] Wall, M.M. (2004). A close look at the spatial structure implied by the CAR and SAR models. Journal of Statistical Planning and Inference. 121, 311-324. https://doi.org/10.1007/s11749-018-0599-x https://cran.r-project.org/package=spData http://www.ggy.bris.ac.uk/personal/LesHepple/bayesian.pdf https://doi.org/10.1007/s13171-021-00246-3 https://www.researchgate.net/publication/339938974 https://doi.org/10.18637/jss.v040.i06