Spatial autocorrelation analysis of tourist arrivals using municipal data : A Serbian example

Spatial autocorrelation methodologies can be used to reveal patterns and temporal changes of different spatial variables, including tourism arrivals. The research adopts a GIS-based approach to spatially analyse tourist arrivals in Serbia, using Global Moran’s I and Anselin’s Local Moran’s I statistics applied on the level of municipalities. To assess feasibility of this approach the article discusses spatial changes of tourist arrivals in order to identify potentially significant trends of interest for tourism development policy in Serbia. There is a significant spatial inequality in the distribution of tourism arrivals in Serbia that is not adequately addressed in tourism development plans. The results of global autocorrelation suggest the existence of low and decreasing spatial clustering for domestic tourist arrivals and high, relatively stable spatial clustering for international tourists. Local autocorrelation statistics revealed different of domestic and international tourism arrivals. In order to assess feasibility of this approach these results are discussed in their significance to tourism development policy in Serbia.


Introduction
Spatial inequality of tourism development is pronounced in Serbia.To date, no tourism development strategy has systematically dealt with this issue.Tourism, as a sector with significant spillover into the rest of the economy, can be of paramount importance for regional development.Serbia, likewise other developing countries, embrace tourism to jumpstart its socio-economic development and to uphold sustainable regional development.In the more recent studies, spatial statistical techniques were used more extensively to analyze tourism development inequality and tourism flow patterns in various countries (see Klepers, Rozite, 2009;Chhetri, et al., 2013;Yang, Wong, 2013;Yang, et al., 2013;Kang, et al., 2014;Xing-zhu, Qun, 2014).
Tourist flows have been studied at different state levels, from micro-level to the -level of a destination.For example, Jansen-Verbeke (1995) used interregional and intra-regional analysis of tourist flows within Europe to assess particular destinations within Europe.Hence, in order to carry out an adequate research in the context in Serbia, the authors needed a quantitative technique for detecting and analyzing temporal and spatial patterns across whole country's territory.The main spatial trends that are to be noted in tourism context are the clusters of popular and unpopular places.To ensure viability of the approach, the level of clustering should be related to the existing administrative division of the country.As a solution, the authors adopted a GIS approach for spatial analysis using spatial autocorrelation analysis based on the data gathered on the municipality level.The first aim was to investigate trends of spatial clustering of municipalities by the level of tourist arrivals between 2001 and 2013.The second aim was to indicate local deviations from global pattern of spatial association, that is, significant local clusters or outliers and how they have changed over the time.For achieving aims of the research Global Moran's I and Anselin's Local Moran's I statistics are applied using Geographic Information System (GIS) software tools.
This paper extends the literature by using contemporary methodological framework within which spatial changes of tourist arrivals can be studied with respect to statistically significant global and local spatial associations in the variables.Unlike other similar studies, we used municipality level data for the analysis of tourist arrivals.Usually, for this kind of tourism studies the researchers used global and local distribution of cities or point data and examined their spatial aggregation analysis (Yüncü, et al., 2016).However, municipality level data appears to be a compromise between data availability, data aggregation level, analytical capabilities and displaying opportunities of the results for effective decision making.Comparing to studies that employ similar methodology, this paper analyses the topic in the context of a European country as there is no similar research concerning this region.Furthermore, in order not to miss statistically significant results in case of local autocorrelation analysis, temporal statistical tests are presented in the isolation at different level of significance, together with results of false discovery rate (FDR) corrected results.

The case area
The 2005 Tourism Development Strategy divided the country into four main spatial tourism clusters: Vojvodina, Belgrade region, South-West Serbia and South-East Serbia (Official Gazette of the Republic of Serbia, 2006).The suggestions of four tourism clusters were based on the rational strongholds and various kinds of economic activities that could be developed in some parts of the country (Ilić, 2013).Bordering areas can be considered as parts of two bordering clusters (see Figure 1).There are two main basic markets for tourism demand: (1) Belgrade, Novi Sad and other cities dominantly visited by international tourists (primarily motivated by business reasons) and (2) spas and mountains that together make up more than half of the tourist traffic.Spas and mountains are visited dominantly by domestic tourists (Ministry of Trade, Tourism and Telecommunications of the Republic of Serbia, 2014).Domestic tourist arrivals in Serbia showed a significant declining trend in the reported period.After the significant fall in 2002, domestic tourist arrivals suffered another "shock" during global economic crises (Najdić, Sekulović, 2012).Domestic tourist arrivals decreased from 1.88 million in 2001 to 1.27 in 2013.On the other hand, international tourist arrivals have fared better.The significant jump recorded in 2007 was followed by small decline during the global economic crisis (Čerović, et al., 2015).However, the strong increasing trend of international arrivals has been continued from 2010 (see Figure 2).

Data collection and processing
For the purpose of this research, already aggregated data for domestic and international tourist arrivals were collected for every municipality in Serbia for the period of 13 years (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013).The data is obtained from Statistical Office of the Republic of Serbia, the only official and hence trustful provider of statistical data.
The authors used 160 polygons that represent territories of municipalities and cities -units of administrative division of Serbia without territory of Kosovo and Metohija.GIS data were analyzed using Esri's ArcGIS suite version 10.2.

Global and local autocorrelation analysis
Spatial autocorrelation can be defined "as the correlation of a variable with itself over space" (Burt, et al., 2009).Spatial autocorrelation is concerned with es-

N
tablishing whether the presence of a variable in one region makes the presence of that variable in neighbouring regions more, or less, likely (Thomas, Huggett, 1980).Spatial autocorrelation can be positive and negative.If similar values of a variable tend to cluster in space, the geographical distribution of that variable can be described as positively spatially auto-correlated.Negative spatial autocorrelation exists when nearby variables are dissimilar (Burt, et al., 2009).Furthermore, measures of spatial autocorrelation may be either global or local.Global measures characterize the nature of spatial autocorrelation for the entire study area using one value that summarizes average trends (Nelson, Boots, 2008).Global autocorrelation assumes homogeneity of trends explored.If the assumption is not held, smaller space units can be defined and local autocorrelation can be calculated for each individual unit.
In this study a global autocorrelation analysis was carried out for each year using Spatial Autocorrelation tool in ArcGIS Spatial Statistics toolbox which implement Global Moran's I index.Global Moran's I shows a formal indication of the degree of linear association between the spatial units and their neighbour (Yu, Wei, 2008) and has been used in some recent tourism research (Zhang, et al., 2011;Yang, Wong, 2013;Luo, Yang, 2013;Kang, et al., 2014).To test null hypothesis which claims that feature values are randomly distributed across the study area Moran's score index is calculated.Moran's I value ranges from -1, indicating perfect dispersion to 1, indicating perfect autocorrelation.Value of Moran's close to 0 suggests a random spatial pattern.To evaluate statistical significance of Moran's I z-score and p-value were used.The critical value of z-score is dependable on desired level of confidence.If standardize percentile of normal distribution α is set to 0.1 then z 0.05 =1.64 for 90 percent confidence intervals.For 95 percent confidence intervals when α= 0.05 then percentile of normal distribution value z 0.025 =1.96.On 99 percent of confidence intervals critical value of z for α=0.01 is z 0.005 = 2.58.Further a high positive z-score for a feature indicates that the surrounding features have similar values (spatial cluster).A low negative z-score indicates that the feature is surrounded by dissimilar values (spatial dissimilarity) (Lwin, Murayama, 2009;Yu, Wei, 2006).
For the purposes of this study identification and mapping of local cluster was carried out using Cluster and Outlier Analysis tool in ArcGIS Spatial Statistics toolbox which implements Anselin's Local Moran's I. Anselin's Local Moran's I statistics was described in 1995 (Anselin, 1995) as a technique to identify clusters of features with values similar in magnitude, but also to identify outliers by comparison to neighbouring features and the mean of the entire population (Sugumaran, et al., 2009).This ArcGIS's tool calculates z-scores and p-values that indicate whether the apparent similarity or dissimilarity is more pronounced than one would expect in a random distribution.A high positive z-score for a feature indicates that the surrounding features have similar values (HH cluster -a municipality with a high value surrounded by municipalities of high values or LL cluster -a municipality with a low value surrounded by municipalities of low values).A low negative z-score for a feature indicates a spatial outlier (HL cluster or LH cluster -a municipality with a high value surrounded by municipalities of low values and vice versa for LH).
To overcome issues with both multiple testing and dependency on large number of data sets we used FDR correction (de Castro and Singer, 2006).In order not to miss meaningful findings, corrected data are presented and discussed together with significant clusters and outliers found in isolated tests on 0.01, 0.05 and 0.10 significance levels.
As tourist arrivals represent spatial process that is associated with each polygon this study uses simple polygon contiguity conceptualizations.In other words, municipalities that share a boundary or a node with a target municipality influence target municipality and therefore their influences are included in computations for that target municipality.Row standardization is used to overcome bias from an imposed aggregation scheme of municipality data.Ultimately, the complete data for tourist arrivals was available only at a municipality level.

Results
In order to reveal temporal changes of global patterns of spatial autocorrelation in domestic and international tourism arrivals Global Moran's I statistics is calculated (Table 1).
From 2001 to 2003 the results of global Moran's I tests for domestic tourist arrivals indicate the clustered patterns.Later, the intense of clustering gradually decline.From 2007 the global pattern does not appear to be significantly different whereas null hypothesis could not be rejected.On the other hand, for international tourist arrivals the results of the global Moran's tests indicate highly clustered pattern for entire period.
Linear trend of Moran's I tests for international tourist arrivals during investigated period remained unchanged compared to decline of linear trend for domestic tourist arrivals (Figure 3).
The results of isolated Cluster and Outlier tests are summarized in Table 2 and 3. Figures 4 and 5 show locations of significant clusters and outliers of municipalities up to three levels of significance (for the first and the last year of investigated period).
Municipalities from Belgrade and South-West regions were found to be the most relevant HH clusters for domestic tourist arrivals based on isolated tests on the significance level of sig<0.01,as presented in Table 2. FDR correction partially confirmed findings of isolated tests.As shown in isolated tests, between 2001 and 2006 only municipalities from Belgrade region created a significant HH cluster.Later, HH cluster of municipalities from South-West region appeared.However, FDR correction confirmed only Užice as significant HH cluster.
Mostly lower levels of significance on isolated tests revealed Novi Sad as HL outlier from Vojvodina region, Niš and Sokobanja as HL outliers from South-East region and Rakovica as LH outlier from Belgrade region, until 2006.
Due to high intense of global spatial autocorrelation for international tourist arrivals the most relevant clusters in isolated tests are shown at sig<0.01 significance level (Table 3).Similar results, except in case of HL outlier from South-East region, are evident with reduced sig.values in the FDR correction.In the period between 2001 and 2013 six out of eight municipalities

Discussion of results
In order to assess feasibility of this approach the spatial changes of tourist arrivals will be discussed to identify potentially significant trends of interest for tourism development policy in Serbia.The starting point for this approach to discussion is a viewpoint of Kang and colleagues (2014) suggesting that implementation of national tourism policies should be manifested in the patterns of tourism activities and the patterns of destination development.
Although established 10 years ago, the above-mentioned four tourism clusters still only exist as a conceptual category.The underlying tourist arrivals have own clustered patterns that do not correspond to official tourism cluster propositions.Despite that fact, the draft of the new Strategy from 2016 to 2025 concludes that former four tourism clusters are still relevant for the next planning period.Klepers and Rozite (2009) emphasize that sometimes it is not rational to organize tourism efforts in a single local administrative territory in isolation from contiguous territories.In line with the above discussion, it is important to acknowledge global as well as local spatial patterns of important tourism flow indicators when evaluating the adequacy of tourism development plans.The approach used in this research is based on municipalities that are units of state division small enough to give a detailed image of this spatial process (Dolinaj, et al., 2009).These insights can be a valuable source of information for Serbian tourism administration within state sector (Cimbaljević, et al., 2014).The findings are especially important for tourism policies development of 116 local tourism organizations operated by municipality or regional governments (OECD, 2012).
Generally, due to lack of centralized organization, financial support and expert knowledge most of them are often left to their own devices.
From the findings, the authors raise some important issues regarding current policymaking implications for Serbia: • Because of separate spatial patterns identified for domestic and international arrivals, different spatial tourism planning approaches could be employed for these niches.• Unfavourable spatial distribution of tourist arrivals is clearly identified.These results are supported by the evidences form other tourism studies (Stankov, et al., 2016;Stankov, et al., 2011).For example, South-East Serbia does not have strong enough cluster of municipalities that can act as a catalyst for stronger tourism development.• Identified clusters and outliers with higher number of tourists (mostly cities) can be used as relevant hotspots for focusing efforts of tourism promotion to spill over positive effects of tourism on neighbouring regions.Some other studies also show that tourism hot-spots tent to be associated with cities (Yang, Wong, 2013;Xing-zhu, Qun, 2014).Therefore, Belgrade, as a capital city and most important population growth hub (Nikitović, et al., 2016) and Novi Sad as constant hotspots of international business travellers can be marked as priorities for addressing this issue.
Just with the results of this study in hand, even more specific questions could be further asked regarding spatial dimension of tourist arrival data.Therefore, the use of spatial autocorrelation method- ology for creating better insight into tourist data for planning activities have to be considered in order to achieve a balance in the tourism development.

Conclusions
The paper makes two main contributions to knowledge.First, it demonstrates suitability of data gathered on municipality level to characterise the spatial patters of tourists.Thus, the novelty of this study is in that it has undertaken a spatial autocorrelation analysis of tourist arrivals at municipality level in the context of a European country where these topics have not been examined extensively.Second, the paper extends the on-going discussion on spatial changes of tourist arrivals with respect to statistically significant global and local spatial associations in the variables using spatial autocorrelation analysis.Other authors, such as Becken and colleagues (2007) emphasize the necessity for the development of methodologies that will allow viewing tourist data spatially, based on statistical modelling in order to reinforce decision making.
The results of local statistics as obtained from cluster and outlier analysis reveal different spatial grouping for domestic and international tourism arrivals in Serbia.Decrease of domestic tourism arrivals was followed by the process of spatial randomization.On the other hand, the increase of international tourist arrivals has not been followed by more significant changes in the spatial patterns.Using similar methodology, Yang and Wong (2013) found a strong positive and significant spatial autocorrelation both for inbound and domestic tourist flows to China's cities, as well as significant hot-spots areas.This research also reinforce the suggestion of Kang and colleagues (2014) that spatial autocorrelation methodologies can be use to analyse spatial outcomes of national tourism policies of countries that undergone some kind of social, political, economic, technological growth or change.This was the case of Serbia.During investigated period Serbia was in the process of economic and political transition with the aim of keeping pace with modern economies and preparing state system for EU membership (Stankov, Dragićević, 2015).
From a practical perspective, this study shows that results of spatial autocorrelation methodologies, based on municipality data, can be used to spatially present changes of tourist arrival in order to make more informed decision in planning.Furthermore, this paper tries to facilitate the adoption of spatial statistics by other tourism researchers.
However, the proposed concept highlights some important limitations.Most importantly, the key challenge for adopting spatial autocorrelation in tour-ism data analysis rests in the scope of spatial units for data aggregation (for example, too small city municipalities that attract large portion of tourists).Another challenge is the complexity of tourism phenomenon itself.Researchers must have as much as wide insight on various processes that could influence the clustering of tourist arrivals.Having that in mind, the results of the presented methodology can successfully be analyzed only in specific space-time context.

Figure 1 .
Figure 1.Main tourism clusters according to Tourism Development Strategy (2005-2015) and major cities in Serbia.

Figure 3 .
Figure 3. Evolution of Moran's I statistics for domestic and international tourist arrivals from 2001 to 2013.

Figure 4 .Figure 5 .
Figure 4. Moran significance map for domestic tourist arrivals in Serbia (without the territory of Kosovo and Metohija) in 2001 and 2013.

Table 1 .
Global Moran's I, z-scores and p-values from 2001 to 2013.

Table 2 .
Municipalities representing significant clusters and outliers using isolated Anselin's Local Moran's I from 2001 to 2013 for domestic tourist arrivals.

Table 3 .
Municipalities representing significant clusters and outliers using isolated Anselin's Local Moran's I from 2001 to 2013 for international tourist arrivals.