See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/266007638

A new procedure for generalized STAR modeling using IAcM approach

Article  in  ITB Journal of Science · July 2012

DOI: 10.5614/itbj.sci.2012.44.2.7

CITATIONS READS

10 159

2 authors:

Utriweni Mukhaiyar Udjianna Sekteria Pasaribu

Bandung Institute of Technology Bandung Institute of Technology
21 PUBLICATIONS   34 CITATIONS    50 PUBLICATIONS   119 CITATIONS   

SEE PROFILE SEE PROFILE

Some of the authors of this publication are also working on these related projects:

Maintenance service and lease contracts View project

Space time analysis & modeling View project

All content following this page was uploaded by Utriweni Mukhaiyar on 09 June 2016.

The user has requested enhancement of the downloaded file.

 ITB J. Sci., Vol. 44 A, No. 2, 2012, 179-192 179

A New Procedure for Generalized STAR Modeling using

IAcM Approach
Utriweni Mukhaiyar & Udjianna S. Pasaribu

Statistics Research Division, Institut Teknologi Bandung,

Jalan Ganesha 10 Bandung, Jawa Barat 40132, Indonesia
Email: [email protected]

Abstract. A new procedure of space-time modeling through the Invers of

Autocovariance Matrix (IAcM) is proposed. By evaluating the IAcM behaviors
on behalf of the Generalized Space-Time Autoregressive (GSTAR) process
stationarity, we may find an appropriate model to space-time data series. This
method can complete the Space-Time ACF and PACF methods for identifying
space-time models. For study case, we apply the GSTAR models to the monthly
tea production of some plantations in West Java, Indonesia.

Keywords: stationarity; space time model; weight matrix; IAcM.

1 Introduction
An observation at a certain spatial location and time may be commonly
influenced by previous-time observations in that location and its neighboring
locations. For example, the oil production of a certain well during this month
can be affected by the production during the previous months at that same
location and at the neighboring wells; the number of today’s criminal actions in
a city may be influenced by the number of criminal actions that happened
during the previous days in that city and the nearest cities; or, then the number
of vehicles at a certain intersection within a certain time range is influenced by
the number of vehicles at that location and related intersections during the
previous time ranges. These examples are the portrait of space-time analysis
problems.

Also, we would like to be able to forecast the oil production of that particular
well during the coming months; how many criminal actions will happen in the
city; and how many vehicles can be found at that particular intersection during
the next range times. Knowing these forecasts, people in charge may take the
necessary actions. We can do forecasting after obtaining the appropriate model.
One way of doing this is space-time modeling.

In this paper we proffer a new procedure for space-time modeling, especially

for Generalized STAR (GSTAR) modeling. For this purpose, we modify the
Pfeifer and Deutsch three stage-iterative procedure for space-time modeling
Received March 6th, 2012, Revised May 10th, 2012, Accepted for publication June 22 nd, 2012.
Copyright © 2012 Published by LPPM ITB, ISSN: 1978-3043, DOI: 10.5614/itbj.sci.2012.44.2.7
180 Utriweni Mukhaiyar and Udjianna S. Pasaribu

(see [1]), which consist of identification, parameter estimation and diagnostic

checking. We improve model selection by checking the stationarity of possible
models using the IAcM. The IACM has been studied for autoregressive models
of time-series by Rahardjo in 2005 [2]. We will explain these processes briefly
in Section 3. Before that, we will describe the GSTAR model and spatial weight
matrices in Section 2. In Section 4, we employ real data of the monthly tea
production in some plantations in West Java. Finally Section 5 contains remarks
and conclusions.

2 Generalized STAR Model

The observation at site i and time t, Yi,t , which is a realization of a stochastic
process, follows the GSTAR model if it can be declared as a linier combination
of past observations for both time and spatial indices. A (N×1)-dimensional
centered process vector Yt  Y1,t ,Y2 ,t , ,YN ,t ' follows the GSTAR(
p; 1 , 2 ,..., p ) model if Yt can be presented as,
p k
Yt    k  W  Yt  k  et

(1)
k 1   0

with k is time autoregressive order (k = 1, 2, ..., p), ℓ is spatial autoregressive

order (ℓ = 0, 1,..., λk), and et is (N×1)-dimensional vector of errors that assumed
have i.i.d normal distribution. The W (ℓ) is ℓ-th order weight matrix which the
main diagonal is zero and the sum of each row is one, and 𝚽𝑘ℓ is (NN)
diagonal matrix which presents autoregressive parameter of k-th time order and
ℓ-th spatial order for each location i = 1, 2, ..., N. We write
(1) (2) (𝑁)
𝚽𝑘ℓ = diag 𝜙𝑘ℓ , 𝜙𝑘ℓ , … , 𝜙𝑘ℓ .

Figure 1 Spatial order/lag in two dimensional radius system. Sites in same

’category of distance’ are in the same spatial order. The category of distance is
classified by a fixed radius value d0.
 A New Procedure of Generalized STAR Modeling 181

The ℓ-th spatial order shows the ordered neighbors of a number of spatial
locations, while spatial weight wij   shows the influence of location j on location
i in spatial lag ℓ-th. The locations which are closer to the reference location (s0)
will have a smaller spatial lag, which is distinguished by fixed distance d0
(Figure 1). Since the configuration in each spatial lag order is distinct, different
spatial lags may give different weight matrices.

2.1 Spatial Weights

There are three types of spatial weight: binary, uniform, and non-uniform.
Binary weight matrix has values 0 and 1 in off-diagonal entries; uniform weight
is determined by the number of sites surrounding a certain site in ℓ-th spatial
order; and non-uniform weight gives unequal weight for different sites. The
element of the uniform weight matrix is formulated as,
 1
 , j is neighbor of i in -th order
wij (  )   ni(  )
0
 , others

ni(  ) is the number of neighbor locations with site-i in  -th order.

The non-uniform weight may become uniform weight when some conditions
are met. One method in building non-uniform weight is based on inverse
distance, which is used by Borovkova in 2008 [3]. The weight matrix of spatial
1
lag k is based on the inverse weights 1+𝑑 for sites i and j whose Euclidean
𝑖𝑗
distance 𝑑𝑖𝑗 lies within a fixed distance range, and otherwise is weight zero.

2.2 Stationarity
GSTAR models are invertible since the vector of observations Yt could be
expressed as a weighted linear combination of past observations with weights
that converge to zero. In order for a GSTAR model to represent a stationary
process, one in which the covariance structure of Yt does not change with time,
certain conditions should be met. We obtain the stationary conditions of the first
and second order of GSTAR models using its IAcM by the following
propositions.

Proposition 1 Suppose that (NN)-dimensional matrix 𝚨 = 𝚽𝟏𝟎 + 𝚽𝟏𝟏 𝐖 and

IAcM, 𝐌𝟏 = 𝐈𝑁 − 𝐀′𝐀, is invers of autocovariance matrix whose elements are
the parameters of GSTAR(1;1) process. If the determinants of all the leading
principal submatrices of IAcM are positive, then GSTAR(1;1) model is
stationary.
182 Utriweni Mukhaiyar and Udjianna S. Pasaribu

Proof. The IAcM is built from the sum square of errors from all past
observations until the first observation. For (N×1)-dimensional vector of
observations, 𝒀1 = 𝑌1,1 , 𝑌2,1 , … , 𝑌𝑁,1 ′ , we can write 𝐘1 ′𝐌1 𝐘1 = 𝟏𝒕=−∞ 𝐞𝑡 ′𝐞𝑡
which has positive and finite values. The 𝐌𝟏 is symmetric and since 𝐘1 ′𝐌1 𝐘1 >
0 then 𝐌𝟏 is positive definite matrix [4]. It implies the determinants of all
leading principal submatrices of 𝐌𝟏 are positive. Furthermore, since 𝐌𝟏 and the
covariance structure is do not change with the time, the GSTAR(1;1) process is
stationary. Q.E.D.

The IAcM is a symmetric positive definite matrix. For GSTAR(1; 1) models
1
with 1  0, we may define A   Φ1 W  , and the IAcM is,


 0

𝐌𝟏 = 𝐈𝑁 − 𝐀′𝐀 (2)
Proposition 2 GSTAR(1; 1) models are stationary if the determinants of all
the leading principal submatrices of M1 are positive.

Proof. Similar to Proposition 1, suppose Y1 is a (N×1)-dimensional vector of

observations, 𝒀1 = 𝑌1,1 , 𝑌2,1 , … , 𝑌𝑁,1 ′. The 𝐘1 '𝐌𝟏 𝐘1 is the collection of
squared errors from all past observations until the first observation. This means
that 𝐘1 '𝐌𝟏 𝐘1 always has positive value and is only equal to null if 𝐘1 is a null
vector. We write,
𝐘1 ′ 𝐈 − 𝐀′𝐀 𝐘1 > 0
𝟐 𝟐
𝐘1 − 𝐀𝐘1 >0
𝟐 𝟐
𝐀𝐘1 < 𝐘1

The absolute eigen value of A must be less than one, which means that the
GSTAR(1; 1) process is stationary. Q.E.D.

The stationarity condition in Proposition 1 and 2 meets the stationarity

condition that was formulated by Wei [5] for Vector AR models. Furthermore
for GSTAR(2; 1, 2) models, we define (NN)-dimensional matrix
k
A k   Φk  W  , k = 1, 2 and the (2N2N)-dimensional IAcM is,


 0

𝐈𝑁 − 𝐀′𝟐 𝐀𝟐 −𝐀′𝟏 −𝐀′𝟐 𝐀𝟏

𝐌𝟐 = (3)
−𝐀𝟏 − 𝐀′𝟏 𝐀𝟐 𝐈𝑁 − 𝐀′𝟐 𝐀𝟐
Proposition 3 GSTAR(2; 1, 2) models are stationary if the determinants of
all the leading principal submatrices of its IAcM are positive.
 A New Procedure of Generalized STAR Modeling 183

Proof. Using a similar procedure as with Proposition 2 and by letting 𝐘2 =

𝑌1,1 , … , 𝑌𝑁,1 , 𝑌1,2 , … , 𝑌𝑁,2 ′, we may consider Proposition 3 as proved. Q.E.D

These stationarity conditions can provide a significance contribution to the

process of identifying possible GSTAR models.

3 Identification Process GSTAR Model

A GSTAR model is a generalization of a STAR model with respect to its
parameters. In a STAR model, the parameters are assumed to be the same for
every location, while a GSTAR model assumes that different locations have
different parameter values. Although both have different assumptions
concerning the parameters, they still have the same procedures in doing the
identification process of space-time models.

3.1 Space Time ACF and PACF

The STACF and STPACF had been employed by Pfeifer and Deutsch [1] in
identifying Space-Time Autoregressive Moving Average (STARMA) models.
Since STARMA and GSTAR models have many things in common, STACF
and STPACF may be applied in identifying GSTAR models. The STACF is
obtained by standardizing the autocovariance function of lag time s for
observations with lag spatial k-th and ℓ-th,  k  s  . Autocovariance function in
matrix notation is

 k  s  
1  ()
  W 
Yt  s 
'
(k )
E W Ys (4)
N  
Since equation (4) is a scalar then we write

 k  s  
1
N
 1
  
tr W(  ) W( k ) E  Yt ' Yt  s   tr W(  ) W( k ) Γ  s 
N
(5)

T s
ˆ  s   Yt Yt  s ' and by substituting it to (5) then
The estimator for Γ  s  is Γ 
t 1 T  s

estimated covariance , ˆ lk  s  holds. Furthermore the STACF, which is defined

as,
 k  s 
k  s  
    0   kk  0  
1/ 2

can be estimated by using the sample covariance,

184 Utriweni Mukhaiyar and Udjianna S. Pasaribu

N T s

ˆ k  s   L Y 
i ,t L Yi ,t  s
k

ˆ k  s  
  i 1 i 1 (6)
ˆ   0  ˆ kk  0  
1/ 2 1/ 2


 N T 
   L  
2 N T

  L Yi ,t
2
k 
Yi ,t 
 i 1 t 1 i 1 t 1 

with L as the spatial lag operator of spatial order  -th, so that
 Yi ,t , 0
 N
L Yi ,t  
 wij Y j ,t ,   1, 2,..., p


 j 1
the STPACF for spatial order λ is the last coefficient of 𝜙ℎℓ (ℓ = 0,1,…, λ and h
= 1,2,…) in space-time Yule Walker equation as defined in [1].

  00 1    00  0   01  0    0   0   00  1  01  1    0   1    10 

     
  10 1    10  0   11  0    1  0   10  1  11  1    1  1    11 
           1  h     
     
   0 1     0  0   1  0      0    0  1  1  1      1    1  (7)
  2    20 
 00    
  10  2     21 
  
   1 0  2  h     

  0  2    2  
     
            
  00  h      p0 
    
  10  h      p1 
   
    h  1 h  2   0     

  0  h     h 

Space time model selection may be obtained by observing autocorrelation

behavior among time lags and spatial lags which are represented by the cut off
and tail off pattern of the STACF and STPACF.

3.2 IAcM as a Part of GSTAR Modeling

In this paper we propose a new procedure for GSTAR modeling by involving
stationarity investigation through the IAcM characteristics as presented in
Proposition 1 to 3. This stage is executed after identifying the model through
the STACF and STPACF, and estimating the parameters. We see it as a first
stage of model validation in addition to checking the mean square of residuals
as the second stage. The residuals are obtained from the differences between
real and estimated observations. The model with the minimum mean square of
residuals will be recommended as the fitting model.The procedure of GSTAR
modeling is shown in the following flow chart (Figure 2).
 A New Procedure of Generalized STAR Modeling 185

This procedure combines and advances the three-stage iterative procedure for
space-time modeling of Pfeifer & Deutsch [1] and time-series modeling of Box
& Jenkins [6]. The new aspect of this procedure is the process stationarity check
with the IAcM *).

START Yes
Transformation/
differentiation of Model Parameter estimation
observations No Both mean and identification (Least Squares
variance are constant Plot of ST-ACF & Estimation)
ST-PACF

Centered Plot of observations Model Validation 1 :

observations stationarity checking Checking stationarity using
(null mean) visually IAcM
Some possible
models

Model is stationary*)
No (Determinants of principal sub-
IAcM are positive)
No

Yes

short time Model Validation 2 :

FINISH Model Yes Model is valid
forecasting Model residual checking.

Figure 2 The new procedure of space time autoregressive modeling.

4 Case Study
For our case study we used the data of the monthly tea production of six
plantation sites (N = 6) in West Java, Indonesia (Purbasari, Santosa, Sedep,
Papandayan, Cisaruni, and Dayeuh Manggung), from January 1992 to
December 2010 (T = 228). We executed space-time modeling by using Matlab
and started by centering the processes in order to obtain the null mean of the
processes. The coordinates of the observed plantations and the distances
between them are displayed in a symmetric (N×N)-dimensional matrix of which
the main diagonal is null (see Figure 3).

From Figure 3 we can draft the numbers of spatial lag order and its members,
and also the weight matrices for each spatial lag. We can build the range of
spatial lag based on the fixed radius system in Figure 1 with a fixed d0 of 0.4
km. The interval distance (0, 0.4] is categorized as spatial lag 1, (0.4, 0.8] as
spatial lag 2, and so on until (0.4(-1), 0.4] as spatial lag ℓ. For this case we got
four spatial lags for eachweight matrix (uniform, euclidean-uniform and
binary). The appropriate weight matrices for the first and second spatial lag are
shown in Table 1.
 186 Utriweni Mukhaiyar and Udjianna S. Pasaribu

y PUR SAN SED PAP CIS DAM

4.5 DAYEUH
Site /1 /2 /3 /4 /5 /6
MANGGUNG
PURBASARI PUR 0.000 0.320 0.750 1.140 1.484 1.401
4.0 SANTOSA SAN 0.320 0.000 0.447 0.832 1.185 1.237
SEDEP 0.750 0.447
SED 0.000 0.391 0.738 0.943
PAPANDAYAN PAP 1.140 0.832 0.391 0.000 0.364 0.901
3.5
CISARUNI CIS 1.484 1.185 0.738 0.364 0.000 0.863
3.0 DAM 1.401 1.237 0.943 0.901 0.863 0.000
x
11.90 12.40 12.90 13.40 13.90

(a) (b)

Figure 3 (a) The coordinate of six observed plantations, and (b) the distances
(in kilometre, km). PUR is Purbasari, SAN is Santosa, SED is Sedep, PAP is
Papandayan, CIA is Cisaruni, and DAM is Dayeuh Manggung.

Table 1 The weight matrices for spatial lag ℓ-th, ℓ = 1,2,3,4. For ℓ = 0, the
weight matrix 𝐖 (0) = 𝐈𝑁 . The sum of each row must be equal to one with a null
main diagonal.
Uniform Euclidean - Uniform Biner
0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0
     
1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0
0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0
W (1) 
0

0 0.5 0 0.5 0 

0 0 0.495 0 0.505 0 
 
0 0 1 0 0 0

0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0
   0 
0   0 
0 0 0 0 1 0   0 0 0 1  0 0 0 1 0 
 0 0 1 0 0 0  0 0 1 0 0 0 0 0 1 0 0 0
     
 0 0 1 0 0 0
 0 0 1 0 0 0 0 0 1 0 0 0
1 / 3 1 / 3 0 0 1/3 0  0.311 0.376 0 0 0.313 0  0 0 0 0 1 0
W (2) 
 0 1 0 0 0 0
 
 0 1 0 0 0 0
 
0 1 0 0 0 0


 0 0  0 0 0 1 0 0 0
0 1 0 0 0 1 0 0 0  
 0 
0 
  0 
 0 0 1 0
 0 0 0 1 0 0  0 0 0 1 0

 0 0 0 1 0 0   0 0 0 1 0 0  0 0 0 1 0 0
     
 0 0 0 0.5 0.5 0   0 0 0 0.544 0.456 0  0 0 0 0 1 0
 1 0 0 0 0 0   1 0 0 0 0 0  1 0 0 0 0 0
W (3) 
1 / 3 1 / 3 0 0 0 1 / 3
 
 0.304 0.355 0 0 0

0.342 

0 0 0 0 0 1

 0 0.5 0 0 0 0.5   0 0.460 0 0 0 0.540  0 0 0 0 0 1
     
 0 0 1 / 3 1 / 3 1 / 3 0   0 0 0.326 0.334 0.340 0  0 0 1 0 0 0 
 0 0 0 0 0.5 0.5   0 0 0 0 0.491 0.509  0 0 0 0 0 1
     
 0 0 0 0 0 1   0 0 0 0 0 1  0 0 0 0 0 1
 0 0 0 0 0 1   0 0 0 0 0 1  0 0 0 0 0 1
W (4) 
 1 0 0 0 0 0 
 
 1 0 0 0 0 0 
 
1 0 0 0 0 0

 1 0 0 0 0 0   1 0 0 0 0 0  1 0 0 0 0 0
   0.482 0.518 0 
  
 0.5 0.5 0 0 0 0   0 0 0 0 1 0 0 0 0 

The appropriate weight of the matrices in Table 1 has similar values. This is
acceptable because they have the same criterion for determining the spatial lag
interval. This happens because there are only six observed sites, which should
be categorized into four levels of spatial authority (p = 4). Since there is
subjective judgment involved in building the weight matrices, we should pay
more attention to the rules, especially when the ratio between the number of
locations and the number of spatial lags is close to 1:1. In these cases, it will be
more difficult to obtain appropriate and obedient weight matrices.
 A New Procedure of Generalized STAR Modeling 187

The first stage in GSTAR modeling is to identify possible models by computing

the STACF and STPACF values as in equation (6) and (7). These values
arepresented in Table 2.

Table 2 The STACF and STPACF values for several time and spatial lag.
STACF STPACF
spatial lag 0 1 2 3 4 spatial lag 0 1 2 3 4
time lag time lag
1 0.4697 0.3102 0.3315 0.1989 0.3329 1 0.3578 -0.0086 0.1975 -0.0984 0.0634
2 0.1554 0.006 0.0105 -0.0339 0.0401 2 0.0191 -0.0385 -0.2131 0.1136 -0.0886
3 0.0699 -0.0779 -0.0912 -0.1194 -0.0532 3 0.1688 0.076 -0.088 -0.0345 -0.0845
4 -0.0102 -0.1232 -0.0859 -0.0821 -0.0374 4 -0.1734 -0.0166 -0.0854 0.0221 0.1226
5 0.0651 -0.0633 0.0083 0.0136 0.0473 5 0.0486 -0.0135 0.1478 -0.0569 -0.0015
6 0.0408 -0.0868 -0.004 0.0331 0.0372 6 -0.0401 -0.0529 -0.0472 0.0905 -0.045
7 0.0032 -0.1322 -0.0478 0.0236 -0.0063 7 0.0488 -0.0038 -0.0297 0.1858 -0.1577
8 -0.0343 -0.1714 -0.1305 -0.0862 -0.0808 8 0.036 -0.0311 -0.0982 0.0292 -0.0511

From Table 2, the STACF values for all spatial lags already have cut off after
first time lag. This alert us to consider Moving Average (MA) factors in the
models. But if we plot those values for larger time lags, it seems like sinusoidal
graphs are coming out. This means that we can temporarily neglect the MA
factors and more focus on the STPACF values that represent autoregressive
terms. This is also because this paper discusses about the GSTAR modeling.
The STPACF values are cut off after the first and second lag time for zero,
second, and third spatial lag. From these lags (both spatial and time) we obtain
some possible models, namely: GSTAR(1;0), GSTAR(1;1), GSTAR(1;3),
GSTAR(2;1,1), GSTAR(2;1,3), etc.

The next step is parameter estimation with the least squares method. To execute
this step we should rearrange the GSTAR model to be a linear structure,
Y=X+e, with 𝐘 = 𝐘1′ , 𝐘2′ , … , 𝐘𝑁′ ′, 𝐗 = diag 𝐗1 , 𝐗 2 , … , 𝐗 𝑁 , 𝚽 =
𝚽1′ , 𝚽2′ , … , 𝚽𝑁′ , and 𝐞 = 𝐞1′ , 𝐞′2 , … , 𝐞′𝑁 ′. For each location we have linear
structure 𝐘𝑖 = 𝐗 𝑖 𝚽𝑖 + 𝐞𝑖 , 𝑖 = 1,2, … , 𝑁, with 𝐘𝑖 = 𝑌𝑖,𝑝 , 𝑌𝑖,𝑝+1 , … , 𝑌𝑖,𝑇 ′,
𝐞𝑖 = 𝑒𝑖,𝑝 , 𝑒𝑖,𝑝+1 , … , 𝑒𝑖,𝑇 ′,
(𝑖) (𝑖) (𝑖) (𝑖) (𝑖) (𝑖)
𝚽𝑖 = 𝜙10 , … , 𝜙1𝜆 1 , 𝜙20 , … , 𝜙2𝜆 2 , … , 𝜙𝑝0 , … , 𝜙𝑝 𝜆 𝑝 ′, and

 Vi (0) V (1)  Vi (,p1)1  Vi (0) Vi (1)  Vi ,0 p 

( )

 ,p 1 i ,p 1 ,0 ,0
( )

 V (0) Vi (1)  Vi (,p1 )  Vi (0) Vi (1)  Vi ,1 p 
Xi   i ,p ,p ,1 ,1

         
 V (0) V (1)  Vi (,T1)1  Vi (0) Vi (1)  Vi ,T p p 
( )
 i ,T 1 i ,T 1 ,T  p ,T  p

(𝑘) 𝑁 (𝑘) (0)

with 𝑉𝑖,𝑡 = 𝑗 =1 𝑤𝑖𝑗 𝑌𝑗 ,𝑡 for 𝑘 ≥ 1 and 𝑉𝑖,𝑡 = 𝑌𝑖,𝑡 (see [3]).
188 Utriweni Mukhaiyar and Udjianna S. Pasaribu

For 𝐗 ′ 𝐗 nonsingular, the least square (LS) estimator may be obtained by

𝚽 = 𝐗′ 𝐗 −1
𝐗′𝐘 (8)
The LS estimator is unbiased with minimum variance.
After obtaining the LS estimators, we can construct the matrix of parameters
(A) and the IAcM for checking the process stationarity. The IAcM determinant
values of possible models with 1st, 2nd time lags and 0th, 1st, 2nd, 3rd spatial lag
are presented in Figure 4.

Based on Proposition 1 to 3, the model that produce at least one negative value
of the leading principal submatrix determinant value in the IAcM, will be
eliminated from appropriate model candidates. Looking at Figure 4, we notice
that models of which at least one of its leading principal matrix has negative
determinant are: GSTAR(2;0,2), GSTAR(2;0,3), GSTAR(2;1,3),
GSTAR(2;3,1), GSTAR(2;3,2), and GSTAR(2;3,3). We may leave these
models and move forward to the next stage, diagnostic checking (Model
Validation 2).

In diagnostic checking stage, the residuals of rest possible models are

compared. Models with least mean square of residuals (MSR) is recommended
as the appropriate models. MSR is formulated as,
1
MSR = 𝐘𝑡 − 𝐘𝑡 ′ 𝐘𝑡 − 𝐘𝑡 (9)
𝑁𝑇−𝑝

𝐘𝑡 is a (NT×1)-dimensional vector of estimated observations that follows

GSTAR(𝑝; 𝜆1 , … , 𝜆𝑝 ) process.
Table 3 MSR of possible GSTAR( 𝑝; 𝜆1 , … , 𝜆𝑝 ) models, p = 1, 2. There are no
patterns in MSR values as the order becomes larger or smaller.
GSTAR(𝑝; 𝜆1 , … , 𝜆𝑝 ) Models
Order MSR Order MSR Order MSR
(1;0) 98.6716 (2;0,1) 96.6197 (2;1,2) 89.9513
(1;1) 89.8531 (2;1,0) 87.3273 (2;2,1) 88.3205
(1;2) 88.0149 (2;2,0) 86.7508 (2;2,2) 88.9287
(1;3) 82.5222 (2;3,0) 81.9143 (2;2,3) 90.6263

The model with the smallest MSR value, GSTAR(2;3,0), can be proposed as the
most appropriate model for this case and be employed to execute the short
forecasting. Since there are some models that have close MSR values, we will
perform an additional step. We let the forecasting stage serve as an advanced
checking diagnostic stage. For efficiency, we choose the five possible models
 A New Procedure of Generalized STAR Modeling 189

that produced the lowest MSR from Table 3. Those models are: GSTAR(2;3,0),
GSTAR(1;3), GSTAR(2;2,0), GSTAR(2;1,0), and GSTAR(1;2).

GSTAR(2;0,1) GSTAR(2;0,2)
GSTAR(1;0) GSTAR(1;1) GSTAR(2;0,3) GSTAR(2;1,0)
GSTAR(2;2,0) GSTAR(2;3,0)
GSTAR(1;2) GSTAR(1;3)
1 1

0.8 0.8

0.6 0.6
Determinant

Determinant
0.4 0.4

0.2 0.2

0 0
0 2 4 6 0 1 2 3 4 5 6 7 8 9 10 11 12

-0.2 -0.2
Leading principal sub matrices -th Leading principal sub matrices -th

(b)
(a)
GSTAR(2;1,1) GSTAR(2;1,2) GSTAR(2;2,1)
GSTAR(2;2,2) GSTAR(2;2,3) GSTAR(2;1,3)
GSTAR(2;3,1) GSTAR(2;3,2) GSTAR(2;3,3)
1

0.8

0.6
Determinant

0.4

0.2

0
0 1 2 3 4 5 6 7 8 9 10 11 12
-0.2
Leading principal sub matrices -th
(c)
Figure 4 The determinant values of all the leading principal matrices of the
IAcM tend to be smaller as the dimensions of the leading principle matrices
become larger, (a) possible models with 1st time lag, (b) possible models with 2nd
time lag and one of spatial lag is 0, and (c) others possible models with 2nd time
lag.
190 Utriweni Mukhaiyar and Udjianna S. Pasaribu

For the forecasting stage, we prepared the ten latest data, collected from January
to October 2011. In the advanced checking diagnostic stage, we estimate these
data using the five selected models in order to obtain what we call forecast-
estimated observations. These forecast-estimated observations will be compared
with real observations using (9), so that we have the MSR of the forecast-
estimated observations. The results of the MSR can be seen in Table 4.

Table 4 MSR of forecast-estimated of five possible appropriate GSTAR

models. The smallest MSR is given by GSTAR(2;3,0) model.
GSTAR MSR
(1;2) 175.3843
(1;3) 139.4374
(2;1,0) 155.788
(2;2,0) 141.6173
(2;3,0) 108.9293

Site 1 - Purbasari Site 2 - Santosa

Quantiles of Input Sample

Quantiles of Input Sample

200 200

100 100

0 0

-100 -100

-200 -200
-3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3
Standard Normal Quantiles Standard Normal Quantiles
Site 3 - Sedep Site 4 - Papandayan
Quantiles of Input Sample

Quantiles of Input Sample

200 200

100 100

0 0

-100 -100

-200 -200
-3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3
Standard Normal Quantiles Standard Normal Quantiles
Site 5 - Cisaruni Site 6 - Dayeuh Manggung
Quantiles of Input Sample

Quantiles of Input Sample

200 200

100 100

0 0

-100 -100

-200 -200
-3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3
Standard Normal Quantiles Standard Normal Quantiles

Figure 5 QQ Plot of Sample Data versus Standard Normal Data for each site i, i
= 1,2,...,6 (from top left to bottom right). Quantile-quantile plot of the sample
quantiles of residuals versus theoretical quantiles from a normal distribution.
 A New Procedure of Generalized STAR Modeling 191

Based on the smallest MSR value, we have GSTAR(2;3,0) as the most

appropriate model. For convenience, we display a quantile-quantile plot of
residuals data versus normal distribution data in Figure 5. Since the plot close to
linear, we may say that the residuals have normal distribution.

In order to do forecasting of (t+1)th observations, the GSTAR(2;3,0) may be

presented as the following equation.
𝐘𝑡+1 = 𝚽10 + 𝚽11 𝐖 (1) + 𝚽12 𝐖 (2) + 𝚽13 𝐖 (3) 𝐘𝑡 + 𝚽20 𝐘𝑡−1 (10)
where 𝐘𝑡+1 is a estimated (N×1)-dimensional vector of observations at time
t+1, 𝐘𝑡+1 = 𝑌1,𝑡+1 , 𝑌2,𝑡+1 , … , 𝑌𝑁,𝑡+1 ′, 𝚽1ℓ is a estimated (N×N)-dimensional
matrix of autoregressive parameters for the first time lag, with spatial lag
ℓ = 0,1,2,3. The ℓ = 0 shows the influence of a location on itself.

Site 1 - Purbasari Site 2 - Santosa

300 300
production

production

200 200

100 100

0 0
200 205 210 215 220 225 230 235 240 200 205 210 215 220 225 230 235 240
time time
Site 3 - Sedep Site 4 - Papandayan
300 300
production

production

200 200

100 100

0 0
200 205 210 215 220 225 230 235 240 200 205 210 215 220 225 230 235 240
time time
Site 5 - Cisaruni Site 6 - Dayeuh Manggung
300 300
production

production

200 200

100 100

0 0
200 205 210 215 220 225 230 235 240 200 205 210 215 220 225 230 235 240
time time
real forecast-real forecast-estimated

Figure 6 Plot of forecast-estimated observations using GSTAR(2;3,0) model.

for each site i, i = 1,2,...,6 (from top left to bottom right).
 192 Utriweni Mukhaiyar and Udjianna S. Pasaribu

The GSTAR(2;3,0) model tells us that a location will be influenced by itself

until the second time lag, while its neighbors can exert influence during the first
time lag within three spatial lags.

The plots of the 28 latest observations and the 10 next forecasted observations,
the total MSR of which is 108.9293, are shown in Figure 6. The forecast-
estimated plots display similar patterns and are close enough to the real
observations.

5 Conclusions
An IAcM is suggested to be used in selecting the approriate space-time models
in completion of STACF and STPACF values. Using an IAcM decreases the
subjectivity of the judgment about which models are appropriate. In the case
study of the monthly tea production of six plantations in West Java, we found
that GSTAR(2;3,0) is the most appropriate model. We also observed that three
types of matrices (uniform, euclidean, and binary) will have the same influence
on the modeling process, since the ratio between the number of spatial lags and
observed locations is close to one. They are supposed to have different results if
the number of observed locationsis much higher than the number of observed
spatial lags.

Acknowledgement
The research of the first author was supported by the Fundamental Research of
DIKTI DIPA ITB program in 2011. The data were supplied by PT. Perkebunan
Nusantara VIII (Persero) in Bandung. The authors are also grateful to Assoc.
Prof Wono Setya Budhi and Khreshna Syuhada, Ph.D for their valuable
discussions.

References
[1] Pfeifer, P.E., & Deutsch, S.J., A Three-Stage Iterative Approach for
Space-Time Modeling, Technometrics, 22(1), pp. 35-47, 1980.
[2] Rahardjo, S. Identifikasi Kestasioneran Proses Autoregresif Secara
Analitis, PhD Dissertation, Mathematics, Inst. Teknol. Bandung, 2005.
[3] Borovkova, S.A., Lopuhaä, H.P., & Nurani, B., Consistency and
Asymptotic Normality of Least Squares Estimators in Generalized Space-
Time Models, Statistica Neerlandica, 62, pp. 482-508, 2008.
[4] Laub, A.J., Matrix Analysis for Scientists & Engineers, SIAM,
Philadelphia, 2005.
[5] Wei, W.W.S., Time Series Analysis: Univariate and Multivariate
Methods 2nd. Ed., Pearson Addison Wesley, Boston, 2006.
[6] Box, G.E.P., Jenkins, G.M. & Reinsel, G., Time Series Analysis,
Forecasting and Control, 3rd ed., Prentice Hall, New Jersey, 1994.

View publication stats