# Application of Generalized Space-Time Autoregressive Model on GDP Data in West European Countries

Authors: Nunung Nurhayati; Udjianna S. Pasaribu; Oki Neswan
Journal: Journal of Probability and Statistics
 (2012)

---

## Abstract

This paper provides an application of generalized space-time autoregressive (GSTAR) model on GDP data in West European countries. Preliminary model is identified by space-time ACF and space-time PACF of the sample, and model parameters are estimated using the least square method. The forecast performance is evaluated using the mean of squared forecast errors (MSFEs) based on the last ten actual data. It is found that the preliminary model is GSTAR(2;1,1). As a comparison, the estimation and the forecast performance are also applied to the GSTAR(1;1) model which has fewer parameter. The results showed that the ASFE of GSTAR(2;1,1) is smaller than that of the order (1;1). However, thet-test value shows that the performance is significantly indifferent. Thus, due to the parsimony principle, the GSTAR(1;1) model might be considered as a forecasting model.


---

## Body

## 1. Introduction

Space-time data are frequently found in many areas of research, for example, monthly tea production from some plants, yearly housing price at capital cities, and yearly per capita GDP (gross domestic product) of several countries in some region. The generalized space-time autoregressive model of order(p;λ1,…,λp),  shortened by GSTAR(p;λ1,…,λp), is one of space-time models characterized by autoregressive terms lagged in the pth order in time and the order of  (p;λ1,…,λp)  in space [1].The term of generalization is associated with the model parameters. When a parameter matrix is diagonal, the GSTAR model is the same as space time autoregressive (STAR) model given by Martin and Oeppen [2] and Pfeifer and Deutsch [3]. The notion of generalization has also been used by Terzi [4] who also generalizes STAR models but in a different context. He generalized the STAR(1;1) by adding the contemporaneous spatial correlation but still preserved the scalar parameters.Whenp=1  and λ1=1, GSTAR(1;1) is called the first order of GSTAR model. The model has interpretation that the current observation in a certain location only depends on the immediate past observations recorded at the location of interest and at its nearest neighbourhood [5]. The order (1;1) is the simplest natural assumption if one wants to forecast future observations in a certain location. The STAR(1;1) model is another simple space time model which also has the same interpretation as GSTAR(1;1) model. However, contrary to the generalization model, its parameters of each spatial order are assumed to be the same for all location though; there is no a priori justification for this assumption [1]. Parameters of GSTAR model can be estimated by the method of least square. This method has been used to model the monthly oil production [5] and to model the monthly tea production [1]. However, when the model was applied to their data, none of the papers included a description about how to assess the model based on forecasting performance which is an important step when the modelling purpose is to build a forecasting model. In this paper, we attempt to put in the idea to optimize the goodness of fit in model selection.This paper is presented as follows. In Section2, the GSTAR model and the parameters least squares estimation is reviewed with the example given for GSTAR(1;1) and GSTAR(2;1,1). To illustrate the estimator properties for finite sample size, simulation study is discussed in Section 3. In the last section, the model is applied to the per capita GDP ratio data in West European countries for the period 1956–1996. The one step ahead forecasting is performed for each model for the period 1997–2006. As a comparison performance measure, it is used the empirical mean of squared forecast error (MSFE) where forecast error is defined as the difference between the actual value and the forecast value.


## 2. The Model

LetZ(t)=(Z1(t),…,ZNt)′  be an N-dimensional vector process with zero mean with N as is a fixed positive integer. GSTAR(p;λ1,…,λp)  process is a space-time process Z(t) which satisfies(2.1)Z(t)=∑k=1p∑l=0λkΦklW(l)Z(t-k)+e(t),
where p is the autoregressive order, λk is the spatial order of the kth autoregressive term,  W(ℓ)=(wij(ℓ))  is an  N×N  matrix of spatial weight for the spatial order  ℓ  which has a zero diagonal, sum of each row is equal to one, and matrix  W(0)  is defined as the identity matrix I.  An  N×N  matrix Φkℓ is a diagonal parameter matrix of temporal lag k and spatial lag ℓ with the diagonal element(ϕkℓ(1),…,ϕkℓ(N)). Finally, e(t) is an error vector at time t which is assumed to be independent normal with zero mean and constant variance.Model parametersϕkℓ(1),…,ϕkℓ(N)  for  k=1,2,… and  ℓ=1,2,…  can be estimated by the least squares (LS) estimation. The procedure estimation and the asymptotic properties of LS estimators have been discussed extensively in Borovkova et al. [1]. The following examples are an illustration on how to find the LS estimator for GSTAR(1;1) and GSTAR(2;1,1) models, respectively.Example 2.1 (LS estimation for GSTAR(1;1) model).
From (2.1), GSTAR(1;1) model can be expressed as
(2.2)Z(t)=Φ10Z(t-1)+Φ11W(1)Z(t-1)+e(t).
Suppose observation Z(t)=(Z1(t),…,ZN(t))′ is given for t=1,2,…,T. By rearranging the component of  Z(t)for  t=1,2,…,T, by location then by time, model (2.2) can be expressed as a linear model
(2.3)Z=XΦ+e,
where Z=(Z1(1),…,Z1(T),…,ZN(1),…,ZN(T)),  X=diag(X1,…,XN)  with
(2.4)Xi=(Zi(0)∑j≠iwijZj(0)Zi(1)∑j≠iwijZj(1)⋮⋮Zi(T-1)∑j≠iwijZj(T-1)),Φ=(ϕ10(1),ϕ11(1),…,ϕ10(N),ϕ11(N))′, and  e(e1(1),…,e1(T),…,eN(1),…,eN(T))′.
Then, the LS estimator for parameter matrixΦ is a solution of
(2.5)(X′X)Φ=(X′Z).Example 2.2 (LS estimation for GSTAR(2;1,1) model).
For orderp=2, λ1 = 1 and λ2 = 1, model (2.1) is called GSTAR(2;1,1) model which can be written as
(2.6)Z(t)=Φ10Z(t-1)+Φ̂11W(1)Z(t-1)+Φ20Z(t-2)+Φ21W(2)Z(t-2).
By rearranging the component of Z(t)for t=2,3,…,T, by location then by time, model (2.6) can also be expressed as a linear model
(2.7)Z=XΦ+e,
where Z and e are defined as in Example 2.1, X=diag(X1,…,XN) with
(2.8)Xi=(Zi(1)∑j≠iwij(1)Zj(1)Zi(0)∑j≠iwij(2)Zj(0)Zi(2)∑j≠iwij(1)Zj(2)Zi(1)∑j≠iwij(2)Zj(1)⋮⋮⋮⋮Zi(T-1)∑j≠iwij(1)Zj(T-1)Zi(T-2)∑j≠iwij(2)Zj(T-2))
and Φ=(ϕ10(1),ϕ11(1),ϕ20(1),ϕ21(1),…,ϕ10(N),ϕ11(N),ϕ20(N),ϕ21(N))′.
Then, the LS estimator for parameter matrixΦ is a solution of
(2.9)(X′X)Φ=(X′Z).


## 3. Simulation Study

Under the stationary assumption, the LS estimator for the GSTAR parameters is a consistent estimator [1]. To get insight into the LS properties for finite sample we performed a Monte Carlo simulation with 1000 replications. Artificial data were generated from GSTAR(1;1) model where the error e(t) is normally distributed with mean 0 and covariance matrix I4. Spatial weight matrix and model parameters that used in the simulation, respectively, were (3.1)W(1)=(00.50.500.5000.50.5000.500.50.50),ϕ10(1)=0.2,ϕ10(2)=0.5,ϕ10(3)=0.3,ϕ10(4)=0.2,ϕ11(1)=0.4,ϕ11(2)=0.3,ϕ11(3)=0.5,ϕ11(4)=0.7
such that the parameter matrix for the generated model was
(3.2)Φ10=diag⁡(0.2,0.5,0.3,0.2),Φ11=diag⁡(0.4,0.3,0.5,0.7)The LS estimator vectorϕ̂=(ϕ̂10(1),ϕ̂10(2),ϕ̂10(3),ϕ̂10(4),ϕ̂11(1),ϕ̂11(2),ϕ̂11(3),ϕ̂11(4)) was calculated using R functions in Algorithm 1. Simulation was carried out for sample size T= 40, 50, 100, 500, 1000, 10000, and for each T, simulations were repeated 1000 times to obtain the estimates average(3.3)ϕ̂̅T=11000∑r=11000ϕ̂T(r)Algorithm 1:R codes for least square estimation of GSTAR model.
# =============================================================# [FUNCTION]: OLS estimation for GSTAR(p;L1,…,Lp) models# =============================================================# 3 dimension zeros matrix# –––––––––––––––-zeros <-function(m,n,p){W<-rep(0, m*n*p)dim(W)<-c(m, n, p)W}# “vec” operatorvec<-function(X){a<-dim(X)Y<- t(X[1, ])for (i in 2:a[1]){Y<- cbind(Y,t(X[i, ]))}t(Y)}#––––––––––# Inverse of matrix#––––––––––inv<-function(X){if(dim(X)[1]! = dim(X)[2]) stop(“THE MATRIX MUST BE SYMMETRIC!!!”)else{if (det(X)==0) stop("THE MATRIX IS SINGULAR!!!")else { n<-dim(X)[1]solve(X)}}}# –––––––––––––––––––––––––––––––# Construction of vector Zi, for each i=1,…,N# –––––––––––––––––––––––––––# Construction of vector Zi, for each i=1,…,N# –––––––––––––––––––––––––––# suppose x = c(p,L1,..., Lp) represent the model orderZi<-function(Zt, x){N<-dim(Zt)[1] #number of sitesT<-dim(Zt)[2]-1 #number of time periodsp<-x [1]Zi<-matrix(0, T-p+1, N)for (i in 1: N)Zi[,i]<-Zt[i,(p+1):(T+1)]Zi}# –––––––––––––––––––––––––––––––-# Construction of matrix Xi, for each i = 1,…,N# ––––––––––––––––––––––-Xi<-function(Zt,x){N<-dim(Zt)[1] #number of sitesT<-dim(Zt)[2]-1 #number of time periodsp<-x[1]La<-x[2: length(x)]r<-lmd+1 # where lmd = the greatest order for weight matricesWZ<-zeros(N, T, r)for (k in 1: r)WZ[,, k]<-W [,, k]%*%Zt[, 1: T]Xi<-zeros((T-p+1),sum(La+1), N)if (p==1){ for (i in 1:N){TR<-WZ[i,p:T,1:(La[1]+1)]Xi[,, i]<-TR}}if (p>=2){for (i in 1:N){TR<-WZ[i,p:T,1:(La[1]+1)]for (s in 2:p)TR<-cbind(TR,WZ[i,(p-s+1):(T-s+1),1:(La[s]+1)])Xi[,, i]<-TR}}Xi}# –––––––––––––––––––––––––––––––-# OLS parameter of GSTAR model# ––––––––––––––gstar<-function(Zt,x){p<-x[1]La<-x[2:length(x)]r<-lmd+1 # where lmd = the greatest order for weight matricesN<-dim(Zt)[1] #number of sitesT<-dim(Zt)[2]-1 #number of time periodsXi<-Xi(Zt,x)Zi<-Zi(Zt,x)coef.OLS<-matrix(0,sum(La+1),N)col.name<-array(0,N)for (i in 1:N){coef.OLS[, i]<-inv(t(Xi[,, i])%*%Xi[,,i])%*%t(Xi[,, i])%*%Zi[,i]col.name[i]<-paste("site",i)}colnames(coef.OLS)<-col.nameround(coef.OLS,4)}# –––––––––––––––––––––––––––––––-# Residuals of GSTAR model# ––––––––––––––# (1). To find the LS estimates only, for example GSTAR(2;1,1), use#    the command:#     > gstar(Zt,c(2,1,1))#    where Zt is data matrix.# (2). To find the estimates, prediction values, and residuals#    vector respectively, call the function by the following#    commands:#       > as.2<-res(Zt,c(2,1,1))#   > as.2$coef#   > as.2$pred#   > as.2$res# –––––––––––––––––––––––––––––––-res<-function(Zt,x){coef<-gstar(Zt,x)Xi<-Xi(Zt,x)p<-x[1]La<-x[2:length(x)]N<-dim(Zt)[1] #number of sitesT<-dim(Zt)[2]-1 #number of time periodsZ.OLS<-matrix(0,T-p+1,N)res.OLS<-matrix(0,N,T-p+1)if (p==1){for (i in 1:N){if (La[1]!=0)Z.OLS[,i]<-Xi[,,i]%*%coef[,i]else Z.OLS[,i]<-Xi[,,i]*coef[,i]res.OLS[i,]<-t(Z.OLS[,i]-Zt[i,(p+1):(T+1)])}}if (p!=1){for (i in 1:N){Z.OLS[,i]<-Xi[,,i]%*%coef[,i]res.OLS[i,]<-t(Z.OLS[,i]-Zt[i,(p+1):(T+1)])}}az<-new.env()az$Xi<-Xi # matrix Xiaz$coef<-coefaz$pred<-t(Z.OLS)az$res<-res.OLSax<-as.list(az)}# –––––––––––––––––––––––––––––––-and empirical mean squared error (MSE)(3.4)MSET=11000∑r=11000‖ϕ̂T(r)-ϕ‖2,where∥ϕ̂T(r)-ϕ∥2=∑i=14(ϕ̂10(i)-ϕ10(i))2+(ϕ̂11(i)-ϕ11(i))2.The result is presented in Table1. It can be seen that the parameters estimates (in average) approaches the true parameters as T increases while the empirical MSE is getting smaller and approaching 0 as  T  increasing. It means that behavior of the LS estimator in the simulation exhibits the consistent property. In general, we can notice that the LS estimation could give fair estimates even for moderate sample size such as  T=40 and T=50.Table 1
The LS estimated values (in average) for the data generated from GSTAR(1;1) model with 1000 replications for various sample sizesT compared to the theoretical parameter.
T= 40T = 50T = 100T = 500T = 1000T = 10000ϕ10(1)=0.20.17990.18220.18910.19960.20040.2000ϕ10(2)=0.50.46390.46840.48190.49710.49870.5002ϕ10(3)=0.30.27860.28150.29420.29660.29830.2998ϕ10(4)=0.20.18020.18280.19290.19830.19860.1998ϕ11(1)=0.40.40000.40250.40280.39990.40020.4005ϕ11(2)=0.30.29610.29440.29700.0050.29980.2999ϕ11(3)=0.50.48800.49060.49310.49880.50010.5002ϕ11(4)=0.70.69890.69800.69610.70080.70020.7002MSE0.02790.02190.01050.00020.00010.0001


## 4. Application of GSTAR Model to the Ratio of Per Capita GDP Data

In this section, we apply the GSTAR model to the ratio of per capita GDP data in 16 West European countries. The data was kindly given by Maddison [6] from Faculty of Economics, University of Groningen, the Netherland, who passed away on  April 24 2010. The per capita GDP of a country is the country GDP value divided by its population, and the per capita GDP of total West Europe is the sum of each West European country GDP divided by the total population in West Europe. The ratio of the per capita GDP of a country is the country per capita GDP value divided by the per capita GDP of total West Europe, multiplied by 100. Hence, the unit of the per capita GDP ratio is the percentage. For the data analysis in the following subsection, we will use the ratio of the per capita GDP data and for simplicity it will be called the GDP ratio data.


### 4.1. Dataset and Preliminary Model Building

The dataset is the GDP ratio data for periods 1955–2006. It consists of 52 observations of 16 dimensional vectors. For the purpose of forecasting the data was grouped into the training data set and test data set. The training data is the first 42 observations that will be used for model building and the test data is the last ten data that will be used in forecasting performance comparison.Clearly, the 42 observations in the training data, depicted in Figure1(a), are not stationary though they tend to converge to the value between 50% and 150%. Therefore, to achieve the zero mean stationary data, the first difference transformation and data centralization must be applied.(a) The ratio of per capita GDP in 16 West European countries for the period 1955–2006 and (b) plot of the centralized data of the difference data.
(a)(b)SupposeYi(t) represents the per capita GDP ratio for country  i,i=1,2,…,16  at time  t,t=-1,0,1,…,50. The first difference transformation of Yi(t) denoted by Di(t) is(4.1)Di(t)=Yi(t)-Yi(t-1)and the centralization of the differenced data is(4.2)Zi(t)=Di(t)-D̅i,whereD̅i is the average of  Zi(t)  at location  i,fori=1,2,…,16. Plot of the 16 series of this transformation is displayed in Figure 1(b) and their behaviour has seemed to represent stationary series.As a preliminary model building, we set some notations used in model (2.1). The length of time period is T= 40, and the number of sites is  N=16. Time period  t=0 will correspond to 1956,  t=1  to 1957,  t=2  to 1958, and so on. Since the time dimension is one, the time lag can be ordered naturally by the sequence of  k=1,2,…. On the other hand, the spatial order  ℓ  may be defined in a different ways because in a two-dimensional space there is no specific order just as in one-dimensional space. For the GDP data, there are 16 countries. The first order and the second order neighbours of the countries are given in Table 2. These are defined based on the geographical location of the countries. The first order neighbours of a country are those which have a common border with the country or within a close distance along a sea route. A second order neighbour of a country is the union of all first order neighbourhood countries of its first order neighbours, excluding itself.Table 2
Geographical neighbourhood of order 1 and order 2.
No.CountryCountries in the 1st neighborhoodCountries in the 2nd neighborhood1Austria6,9,152,3,5,7,10,142Belgium5,6,10,161,3,8,9,13,14,153Denmark6,10,11,141,2,4,5,15,164Finland11,143,65France2,6,9,13,15,161,3,7,8,10,12,146Germany1,2,3,5,10,14,154,9,11,13,167Greece91,5,158Ireland162,5,109Italy1,5,7,152,6,13,1610Netherland2,3,6,161,5,8,11,14,1511Norway3,4,146,1012Portugal13513Spain5,122,6,9,15,1614Sweden3,4,6,111,2,5,10,1515Switzerland1,5,6,92,3,7,10,13,14,1616United Kingdom2,5,8,103,6,9,13,15Supposeni(ℓ),ℓ=1,2  represents the number of countries which are the  ℓth neighbours of a country i. The spatial weight of order ℓ  between countries  i  and j can be defined as(4.3)wij(l)={1ni(l),ifjisthelthneighboursofi,0,ifjisnotthelthneighboursofi.
For example, from Table 2 it can be seen that Austria has 3 the first order neighbours, Germany, Netherland and Switzerland, and has 6 second order neighbours, Belgium, Denmark, France, Greece, Netherland, and Sweden. Then, the first order of spatial weight between Austria and each nearby country is(4.4)w16(1)=w19(1)=w1,15(1)=13and the second order of spatial weight isw12(2)=w13(2)=w15(2)=w17(2)=w1,10(2)=w1,14(2)=1/6.


### 4.2. GSTAR Model Building for the GDP Data

After transforming the data and constructing spatial weight matrix, the next step is identification of the model order. In STAR model-building, [3] used the sample space time autocorrelation function (STACF) and space time partial autocorrelation function (STPACF) as the primary tools for model identification. The order (p;λ1,…,λp) model can be characterized by the tail-off behaviour of the autocorrelations and the cut-off behaviour of the partial after  p  time lag and  λp  spatial lag. Since STAR model is a special case of GSTAR model, these autocorrelation functions will be adopted to identify the order of GSTAR model.Figure2 presents sample STACF and STPACF for the differenced data  Z(t)  up to time lag 10 and spatial lags 0, 1, and 2. The pattern is not clearly suggested an exact order. However, since the sample STPACF cut off after time lag 2 and spatial lag 1, order (2;1,1) can be considered as the space-time order candidate. In addition, the space-time partials also cut off at time lag 5 and spatial 2, but applying this GSTAR model to the data will result too many estimated parameters because there will be at least 160 parameters which have to be estimated.Figure 2
Sample space-time ACF and PACF of the differenced data.The 64 parameters and the error variance in this model were estimated using the least square method and the result is presented in Table3. The empirical MSE for GSTAR(2;1,1) model is 2.735 counted based on 16×39  or 624 values. The residuals histogram and normal probability plot on Figures 4(a) and 4(b) show that the GSTAR(2;1,1) residuals are approximately normal distributed with zero mean and constant variance. Meanwhile, fitted value versus residuals plot in Figure 4(c) exhibits that the residuals do not show a significant pattern. From Figure 3 we can observe that the STACF of the residuals is significantly almost zero except for time lag 5 and 10, and spatial lag 2. The exception points at time lags 5 and 10, suggesting that the seasonal difference of order 5 might be useful for further model analysis. But the seasonal model analysis is not discussed here because it is out of the research scope.Table 3
Least square estimates for GSTAR(2;1,1) parameter.
Siteiϕ̂10(i)ϕ̂11(i)ϕ̂20(i)ϕ̂21(i)1−0.0830.0580.045−0.0462−0.0920.0280.4040.0523−0.2180.730−0.1130.63840.566−0.534−0.2620.25350.1520.563−0.0300.25560.1760.211−0.0520.04670.106−0.048−0.0180.26180.332−0.126−0.1140.54090.271−0.424−0.2170.46110−0.026−0.2500.097−0.332110.549−0.278−0.1660.239120.4420.198−0.254−0.204130.3760.242−0.0540.037140.390−0.120−0.4570.304150.0420.564−0.250−0.949160.2930.2260.119−0.318σ̂LS2=2.969;MSE=2.735.Figure 3
Sample space-time autocorrelation of the residuals from the GSTAR(2;1,1) model.Figure 4
Histogram and normal probability plot of the GSTAR(2,1,1) residuals (a) and (b) and plot of residuals versus fitted value (c).For forecasting purpose, the estimated parameters in Table3 can be used to predict the h-step-ahead forecast value at time forecast origin T defined by(4.5)ẐT(h)=Φ̂10ẐT(h-1)+Φ̂11WT(1)ẐT(h-1)+Φ̂20ẐT(h-2)+Φ̂21WT(2)ẐT(h-2).


## 4.1. Dataset and Preliminary Model Building

The dataset is the GDP ratio data for periods 1955–2006. It consists of 52 observations of 16 dimensional vectors. For the purpose of forecasting the data was grouped into the training data set and test data set. The training data is the first 42 observations that will be used for model building and the test data is the last ten data that will be used in forecasting performance comparison.Clearly, the 42 observations in the training data, depicted in Figure1(a), are not stationary though they tend to converge to the value between 50% and 150%. Therefore, to achieve the zero mean stationary data, the first difference transformation and data centralization must be applied.(a) The ratio of per capita GDP in 16 West European countries for the period 1955–2006 and (b) plot of the centralized data of the difference data.
(a)(b)SupposeYi(t) represents the per capita GDP ratio for country  i,i=1,2,…,16  at time  t,t=-1,0,1,…,50. The first difference transformation of Yi(t) denoted by Di(t) is(4.1)Di(t)=Yi(t)-Yi(t-1)and the centralization of the differenced data is(4.2)Zi(t)=Di(t)-D̅i,whereD̅i is the average of  Zi(t)  at location  i,fori=1,2,…,16. Plot of the 16 series of this transformation is displayed in Figure 1(b) and their behaviour has seemed to represent stationary series.As a preliminary model building, we set some notations used in model (2.1). The length of time period is T= 40, and the number of sites is  N=16. Time period  t=0 will correspond to 1956,  t=1  to 1957,  t=2  to 1958, and so on. Since the time dimension is one, the time lag can be ordered naturally by the sequence of  k=1,2,…. On the other hand, the spatial order  ℓ  may be defined in a different ways because in a two-dimensional space there is no specific order just as in one-dimensional space. For the GDP data, there are 16 countries. The first order and the second order neighbours of the countries are given in Table 2. These are defined based on the geographical location of the countries. The first order neighbours of a country are those which have a common border with the country or within a close distance along a sea route. A second order neighbour of a country is the union of all first order neighbourhood countries of its first order neighbours, excluding itself.Table 2
Geographical neighbourhood of order 1 and order 2.
No.CountryCountries in the 1st neighborhoodCountries in the 2nd neighborhood1Austria6,9,152,3,5,7,10,142Belgium5,6,10,161,3,8,9,13,14,153Denmark6,10,11,141,2,4,5,15,164Finland11,143,65France2,6,9,13,15,161,3,7,8,10,12,146Germany1,2,3,5,10,14,154,9,11,13,167Greece91,5,158Ireland162,5,109Italy1,5,7,152,6,13,1610Netherland2,3,6,161,5,8,11,14,1511Norway3,4,146,1012Portugal13513Spain5,122,6,9,15,1614Sweden3,4,6,111,2,5,10,1515Switzerland1,5,6,92,3,7,10,13,14,1616United Kingdom2,5,8,103,6,9,13,15Supposeni(ℓ),ℓ=1,2  represents the number of countries which are the  ℓth neighbours of a country i. The spatial weight of order ℓ  between countries  i  and j can be defined as(4.3)wij(l)={1ni(l),ifjisthelthneighboursofi,0,ifjisnotthelthneighboursofi.
For example, from Table 2 it can be seen that Austria has 3 the first order neighbours, Germany, Netherland and Switzerland, and has 6 second order neighbours, Belgium, Denmark, France, Greece, Netherland, and Sweden. Then, the first order of spatial weight between Austria and each nearby country is(4.4)w16(1)=w19(1)=w1,15(1)=13and the second order of spatial weight isw12(2)=w13(2)=w15(2)=w17(2)=w1,10(2)=w1,14(2)=1/6.


## 4.2. GSTAR Model Building for the GDP Data

After transforming the data and constructing spatial weight matrix, the next step is identification of the model order. In STAR model-building, [3] used the sample space time autocorrelation function (STACF) and space time partial autocorrelation function (STPACF) as the primary tools for model identification. The order (p;λ1,…,λp) model can be characterized by the tail-off behaviour of the autocorrelations and the cut-off behaviour of the partial after  p  time lag and  λp  spatial lag. Since STAR model is a special case of GSTAR model, these autocorrelation functions will be adopted to identify the order of GSTAR model.Figure2 presents sample STACF and STPACF for the differenced data  Z(t)  up to time lag 10 and spatial lags 0, 1, and 2. The pattern is not clearly suggested an exact order. However, since the sample STPACF cut off after time lag 2 and spatial lag 1, order (2;1,1) can be considered as the space-time order candidate. In addition, the space-time partials also cut off at time lag 5 and spatial 2, but applying this GSTAR model to the data will result too many estimated parameters because there will be at least 160 parameters which have to be estimated.Figure 2
Sample space-time ACF and PACF of the differenced data.The 64 parameters and the error variance in this model were estimated using the least square method and the result is presented in Table3. The empirical MSE for GSTAR(2;1,1) model is 2.735 counted based on 16×39  or 624 values. The residuals histogram and normal probability plot on Figures 4(a) and 4(b) show that the GSTAR(2;1,1) residuals are approximately normal distributed with zero mean and constant variance. Meanwhile, fitted value versus residuals plot in Figure 4(c) exhibits that the residuals do not show a significant pattern. From Figure 3 we can observe that the STACF of the residuals is significantly almost zero except for time lag 5 and 10, and spatial lag 2. The exception points at time lags 5 and 10, suggesting that the seasonal difference of order 5 might be useful for further model analysis. But the seasonal model analysis is not discussed here because it is out of the research scope.Table 3
Least square estimates for GSTAR(2;1,1) parameter.
Siteiϕ̂10(i)ϕ̂11(i)ϕ̂20(i)ϕ̂21(i)1−0.0830.0580.045−0.0462−0.0920.0280.4040.0523−0.2180.730−0.1130.63840.566−0.534−0.2620.25350.1520.563−0.0300.25560.1760.211−0.0520.04670.106−0.048−0.0180.26180.332−0.126−0.1140.54090.271−0.424−0.2170.46110−0.026−0.2500.097−0.332110.549−0.278−0.1660.239120.4420.198−0.254−0.204130.3760.242−0.0540.037140.390−0.120−0.4570.304150.0420.564−0.250−0.949160.2930.2260.119−0.318σ̂LS2=2.969;MSE=2.735.Figure 3
Sample space-time autocorrelation of the residuals from the GSTAR(2;1,1) model.Figure 4
Histogram and normal probability plot of the GSTAR(2,1,1) residuals (a) and (b) and plot of residuals versus fitted value (c).For forecasting purpose, the estimated parameters in Table3 can be used to predict the h-step-ahead forecast value at time forecast origin T defined by(4.5)ẐT(h)=Φ̂10ẐT(h-1)+Φ̂11WT(1)ẐT(h-1)+Φ̂20ẐT(h-2)+Φ̂21WT(2)ẐT(h-2).


## 5. GSTAR(1;1) Modeling

GSTAR(1;1) model is the simplest model of GSTAR(p;λ1,…,λp) model class defined in (2.1) because it is only characterized by the autoregressive terms lagged in time and spatial of order one. From (2.1) we can write GSTAR(1;1) as (5.1)Zi(t)=ϕ0iZi(t-1)+ϕ1i∑j=1NwijZj(t-1)+ei(t),
or in matrix notation (5.2)Z(t)=Φ10Z(t-1)+Φ11W(1)Z(t-1)+e(t),
where  Z(t)=(Z1(t),…,ZN(t)),  Φ10=diag(ϕ10(1),…,ϕ10(N)), and  Φ11=diag(ϕ11(1),…,ϕ11(N)). To simplify the notation,  Φ10  and Φ11  are frequently written, respectively, by  Φ0 and  Φ1.The model has an interpretation that the current observation in a certain location only depends on the immediate past observations recorded at the location of interest and at the nearest locations. The GSTAR(1;1) model has2N parameters since it is assumed that the parameter for each location is allowed to be unequal. For the per capita GDP data the least square estimator for GSTAR(1;1) model is presented in Table 4. The residuals distribution of the model, presented in Figure 5, give the conclusion that the residuals are approximately normally distributed with zero mean and the variance is nearly constant. In general, the STACF plot in Figure 6 shows that the residuals are significantly uncorrelated. The insignificant value was only found at the last time lag.Table 4
Least square estimates for GSTAR(1;1) parameters.
Siteiϕ̂10(i)ϕ̂11(i)1−0.0730.0862−0.0900.2503−0.2010.87140.422−0.45250.1920.58260.1870.22170.1160.01180.3550.07790.183−0.261100.012−0.123110.465−0.175120.3580.054130.3170.287140.285−0.162150.0060.119160.3400.203σ̂LS2=3.1556;MSE=2.992.Figure 5
: GSTAR(1;1) residuals histogram (a), normal probability (b), and residuals versus fitted values (c).Figure 6
Space-time autocorrelation of the GSTAR(1;1) residuals.For the GDP data case, GSTAR(2;1,1) model has 64 parameters while GSTAR(1;1) has 32 parameters. Hence, it is not a surprise if the empirical MSE of GSTAR(2;1,1) is less than that of GSTAR(1;1). However, though the number of parameters of GSTAR(1,1) is half of the other one, the empirical MSE of GSTAR(1,1) is only decreasing 0.257 compared to the GSTAR(2;1,1).The distribution of the MSE difference for each country is presented in a bubble plot (Figure7). Center of the bubble is the value of MSE difference and the radius is its absolute value. The bubble placed under the zero axis indicates that the GSTAR(2,1,1) MSE for the associated country is smaller than that of the GSTAR(1;1) MSE. From the figure, we can see that the decrease of MSE values is mostly contributed by five countries, Belgium, Denmark, Ireland, Sweden, and Switzerland. Though the MSE value seems different, the paired t-test for the result gives the  P values < 0.678. It means that the MSE values for both models are significantly indifferent.Figure 7
Bubble plot for the difference value of empirical MSE between GSTAR(2;1,1) and GSTAR(1;1). The bubble placed under the zero axis indicates that GSTAR(2,1,1) has a smaller value of MSE difference than that of the GSTAR(1;1) model.


## 6. Comparison of Forecast Performance

For the purpose of forecasting model, it would be useful if we also consider their forecast performance. Therefore, in this section we will examine the one-step-ahead forecasting performance for each model candidate using the last ten actual data points of the per capita GDP ratio data set. Result of this section is expected to become a supplementary reference in finding the most parsimony space-time model for the case of per capita GDP ratio.In (4.5) we have defined the h-step-ahead forecasting for GSTAR(2;1,1) model. For GSTAR(1;1) model, it can be estimated by  (6.1)ẐT(h)=Φ10ẐT(h-1)+Φ11W(1)ẐT(h-1),(6.2)ẐT(h)=0.1674ẐT(h-1)+0.0453W(1)ẐT(h-1),
where Φ10=diag(ϕ10(1),…,ϕ10(N)),   Φ11=diag(ϕ11(1),…,ϕ11(N)),  and for  h≤0, ẐT(h)=Z(T-h).The one-step-ahead forecastẐT+j-1(1)  of each data points can be calculated using (4.5), (6.1), and (6.2) for  h=1. Note that calculation of the forecast value for each time  j=1,2,…,10, is performed without updating the parameter estimates and average data. The one-step-ahead forecast for GSTAR(2;1,1) model is  (6.3)ẐT+j-1(1)=Φ10Z(T+j-1)+Φ11W(1)Z(T+j-1)+Φ20Z(T+j-2)+Φ12W(1)Z(T+j-1),
where  j=1,2,…,10. Meanwhile the one-step-ahead forecast for GSTAR(1;1) model is  (6.4)ẐT+j-1(1)=Φ10Z(T+j-1)+Φ11W(1)Z(T+j-1).To compare the forecast performance, we use mean of square forecast error (MSFE) which is defined by(6.5)MSFE=1N∑i=1NE̅i,
where  (6.6)E̅i=110∑j=110[Yi(T+j)-Ŷi,T+j-1(1)]2,Ŷi,T+j-1(1)=(Ẑi,T+j-1(1)+Di)+Yi(T+j-1),
and  Ẑi,T+j-1(1)  is the element of  Ẑi,T+j-1(1).To measure the performance closeness between two models, we also calculate the MSFE difference between modelM1 and model  M2 for country  i which is defined by  (6.7)Ci=E̅i(M1)-E̅i(M2).The performance ofM1 and  M2  model is said to be close if the value of Ci is near to zero. The negative value of Ci indicates that for location i, the  M1 model is better than M2  model. Since the residual is approximately normally distributed, the paired t-test can also be applied to test the null hypothesis that model  M1  and model  M2  have the same forecast performance.SupposeM1  and  M2  in (6.7) are represented by GSTAR(2;1,1) and GSTAR(1;1) model, respectively. The MSFE difference values are displayed in a bubble plot (Figure 8). Center of the bubble is the value of MSFE difference and the radius is its absolute value. The bubble placed under the zero axis indicates that the MSFE value for model  M1  is smaller than that of the model M2. From the figure, it can be seen that most of the countries have the MSFE value which is scattered around the zero axis. The high value of the P value < 0.9858 for the paired t-test gives the evidence that the forecast performance of each model is not significantly different.Figure 8
Bubble plot for the difference of mean of square forecast error (MSFE) between GSTAR(2;1,1) and GSTAR(1;1). The bubble under the zero axis indicates that the GSTAR(2,1,1) has a smaller value of MSFE difference than that of the GSTAR(1;1) model.The behavior of both differences for each time is displayed in Figure9. From the figure, it can be observed that most points in Ireland (site 8) have a negative value. It means that in Ireland GSTAR(2;1,1) model has a better performance than GSTAR(1;1). Conversely, the most points in Norway (site 11) and Switzerland (site 15) are positive when the difference is counted for the pair of GSTAR(2;1,1) and GSTAR(1;1). On the other hand, Figure 7 also shows that in Norway there is a point with high jump difference value. It means that the worst performance of GSTAR(2;1,1)  occurred at this point. In addition, the ratio of the negative value and the positive value is almost the same. It gives the interpretation that forecast performance of GSTAR(2;1,1) and GSTAR(1;1) model is almost similar.Figure 9
Square forecast error difference (SFED) for GSTAR(2;1,1) and GSTAR(1; 1).


## 7. Conclusions

GSTAR modeling has been built to the ratio of per capita GDP data in West European countries. The model of order (2;1,1) model has been identified as the candidate model. However, when the forecast performance is compared with GSTAR(1;1), it is found that the performance is significantly indifferent. Due to parsimony principle, we recommend that the GSTAR(1;1) might be considered as a forecasting model.



---

*Source: 867056-2012-05-08.xml*