Changes in tourist arrivals in Tuscan destinations between 2000 and 2013: A group based trajectory approach

The aim of this paper is to examine the development of tourism in Tuscany’s municipalities, using a methodology to detect common growth curves from yearly time series (2000-2013) of the number of arrivals. The analysis is based on a partition of Tuscany into 262 areas derived from a partial aggregation of the original 287 municipalities. The work is carried out by finding a criterion to detect common time patterns of arrivals across territories, classifying territories according to such similarities in patterns, and verifying whether areas in the same group are also spatially contiguous. The study confirms the presence of heterogeneity in the growth patterns of different areal units. There is a general stability of the classical tourism destinations which are renowned worldwide. Other destinations have been developing due to exploiting the demand for rural tourism, interest in nature, environment and landscape, and so on. The growth patterns are derived by means of a special approach of finite mixture modelling: the group based trajectory model, a semiparametric method that aims to detect common patterns in longitudinal data.


Introduction
The aim of this paper is to examine the development of tourist municipalities in Tuscany through a methodology to detect common growth curves from yearly time series of the number of arrivals. The analysed period starts in 2000 and ends in 2013. Though a mature destination, Tuscany consists of a large variety of tourist products (cultural tourism, spas, mountain resorts, coastal areas, religious sites etc.) due to its traditional, cultural and environmental resources. In fact, discovering and understanding similarities in tourist arrivals across local territories and over time can shed light on the different developmental phases of the various destinations, even when the time series is relatively short. Furthermore, geographical proximity should also be taken into account as it can play an important role.
Recently, a study of short time series of tourist arrivals was carried out in Gursoy et al. (2013), where the methodology of splines was used to detect common growth patterns. This paper follows the steps suggested in Gursoy et al. (2013) but apply a special approach of finite mixture modelling. More specifically, the work is carried out by: (1) finding a criterion to detect common time patterns of arrivals across territories; (2) classifying territories according to such similarities in patterns; (3) verifying whether areas in the same group are also spatially contiguous.
Points (1) and (2) are accomplished through the application of the Latent Class Growth Analysis (LCGA) or Growth Based Trajectory Modeling (GBTM), a special case of the Growth Mixture Modelling or GMM (Muthén, 1989;Muthén and Muthén, 2000;Nagin, 1999;. Point (3) is accomplished by means of a special test which compares the number of contiguous territories in the same cluster with the number that would occur in the case of randomness (Gursoy et al., 2013). The queen contiguity is chosen for this purpose. Point (3) aims to investigate the presence of any spatial homogeneity in the growth patterns across different areal units. The analysis is based on a partition into 262 areas derived from a partial aggregation of the original 287 municipalities of Tuscany (due to the occurrence of zeros and confidentiality issues).
The importance of studying the evolution of tourism at the administrative territorial level, and for all Tuscan municipalities, resides not only in scientific interest but is also based on two facts: one substantial and the other formal fact.
1) Because of the huge number of resources (artistic, cultural, environmental etc.) distributed throughout the entire territory, any municipality can potentially or effectively be a tourism destination; in fact, as we will see below, even marginal territories exhibit a systematic increase in the tourist accommodation facilities and tourist flows over time. Moreover, tourism is also promoted as a driver to develop and reinforce marginal areas and other economic activities (for example, agritourism is a source of relevant revenues for agriculture and rural areas). 2) For the regional government, any Tuscan municipality may apply a special overnight tax and, in this sense, be considered a tourist destination. The only requirement is that the municipality constitute a special tourist observatory to guarantee the management of tourism activities in a sustainable and competitive perspective, in compliance with the indications of the European Union (Dupeyras and MacCallum, 2013;European Commission, 2013;European Communities, 2006;2006a) and the NECSTouR network (www.necstour.eu).
At the end of 2013, the tourist tax was applied in more than 201 of the 287 municipalities (Federalberghi, 2014). The Tuscany Region has provided special financial funds for the municipalities, in order to build the local Tourist Destination Observatory (TDO). The TDO project has also opened the way for the construction of a repertory of comparable data at a municipality level.
In this respect, the paper aims to add a new prospective view of the Tuscan territory: the different growth trajectories of tourism municipalities compared one to the other and with the macro trend and the extent of their geographical heterogeneity may provide important knowledge for the adoption of appropriate growth strategies. In fact, the municipality level is rarely used in systematic analyses of tourism flow patterns although a territory like this is more adequate in the competitive and sustainable development of tourism. This is demonstrated by a large number of contributions -from both academic and non-academic institutions -concerning the implementation of systems of tourism indicators at the municipality level for the management of tourism (Torres-Delgado, López Palomeque, 2012;.
The paper is structured as follows. Sections 2 and 3 contain a brief review of the methodology related to the issue of analysing tourism flows and identifying common trajectories in time series data, while section 4 describes the methodology in detail. Section 5 presents Tuscany data on tourism. Sections 6 and 7 contain the results of the analysis, and section 8 some concluding remarks.

Literature review
Analysis of tourism flows Since the early 1990s, an impressive number of contributions have gradually shed light on the nature and characteristics of tourism competitiveness by analysing the measures of competitiveness, the determinant attributes, and the level (regional, local, firm) of the analysis.
The number of arrivals is often used as a measure of tourism competitiveness though some authors are not in agreement, mainly because of the complexity of the tourism system and the necessity to also address sustainability issues (Hall and Butler, 1995;Hall et al., 2004;Crouch, 2011). In fact, there are various key determinants of competitiveness, also depending on the evolutionary stage of the destination (Ryan, 1991;Crouch and Ritchie, 1999;Buhalis, 2000;Ritchie and Crouch, 2003;Vanhove, 2005;Wilde and Cox;Sharpley, 2009;European Commission, 2006WTO, 2009;Barros et al., 2011;Romão, Guerreiro, Rodrigues, 2013). Nevertheless, tourist flows are seen as an indirect measure of economic returns from tourism activities and for that reason, the issue of explaining the time evolution of tourist flows has been addressed by a large number of studies. Many works are based on yearly time series in order to identify the evolutionary phases of the lifecycle (Ma and Hassink, 2013), or on monthly data to study seasonality (De Cantis et. al., 2011;Giusti and Grassini, 2009).
More sophisticated econometric models explaining tourist flows are employed in order to compare competitiveness among destinations, as discussed by Zhang and Jensen (2006). Moreover, the time pattern of arrivals or nights spent is analysed over space to account for the heterogeneity in the growth patterns of different regions (Yang & Fik, 2014) as the latter might be found at different stages of their lifecycle (Butler, 1980), and/or different mechanisms may exist across locations in exploiting the resources or providing tourist services.
In Italy, some recent studies of tourism flows have been carried out at the level of administrative (regions-NUTS2, provinces-NUTS3) and functional areas (local tourism systems). In particular, Cracolici et al. (2006) explain nights spent in the Italian provinces within the framework of frontier and TFP (Total Factor Productivity) analyses. Cracolici et al. (2008) study the attractiveness of competing regions in southern Italy by also considering subjective tourist data. Massidda and Etzo (2012) apply a dynamic panel data analysis to estimate the main determinants of domestic tourism in the Italian regions. Morrocu and Paci (2013) carry out an econometric study of tourism flows for the Italian provinces based on origin-destination spatial interaction models. Patuelli et al. (2012; study the effect of World Heritage Sites and the cultural offer on tourist arrivals in the Italian regions through a spatial interaction model. Lorenzini et al. (2014) study the determinants of domestic tourist arrivals in Italian provinces, within the framework of gravity models with supply-and demand-side covariates.
Other studies deal with the analysis of common growth patterns in tourism. In particular, Gursoy et al. (2013) analyse tourist arrivals by using a functional classification of the Italian territory in tourism systems (sets of municipalities sharing a common typology of tourism). Recently, Giusti and Grassini (2014) examined the annual change rate of tourist arrivals in Tuscan municipalities in order to detect any spatial dependence.
In conclusion, we can observe that the analysis at the level of administrative areas like regions and provinces (over-municipality level) is prevailing because more available and accurate information is provided for these territorial levels.

Identifying common patterns in time series data
In the perspective of clustering time series data, several approaches can be followed and they differ in relation to the length of the series, the sampling intervals of the data (regular or irregular), and the field of research (Liao, 2005).Some approaches work directly on raw data or a synthetic representation of these, as time series data are static data, such as the procedure proposed by Genolini and Falissard (2010) which is based on the k-means clustering algorithm.
In some approaches (Liao, 2005), the raw time series is firstly converted into a number of model parameters to which a clustering algorithm is then applied. Recently, in a study of annual tourist arrivals, the application of Bspline curves and the subsequent clustering of spline parameters through the k-means algorithm (Abraham et al., 2003) was experimented in Gursoy et al. (2013). The use of B-splines in the representation of the single trajectory is a flexible tool, however it gives rise a large number of parameters to be clustered.
A more sophisticated approach in clustering longitudinal data is Growth Mixed Modeling or GMM (Muthén, 1989;Muthén, 2004), which is introduced when analyzing panel data, but it has also been used with time series (Cabral et al., 2006).GMM typically refers to statistical methods that allow for estimating interindividual variability and intra individual patterns of change (Curran, Obeidat, Losardo, 2010).
Latent Class Growth Analysis (LCGA) is a special case of GMM, in which the variances of latent slope and intercept are fixed at zero in the class (Jung and Wickrama, 2008;Andruff et al., 2009). It follows that all individuals within the same class have the same growth trajectory, only allowing for individual disturbances. The LCGA framework has been extensively developed by Nagin and colleagues (Nagin, 1999;Nagin & Land, 1993;Weisburd et al., 2004) in their study of criminal behavior. With Nagin's approach, this methodology is known as Group Based Trajectory Modelling (GBTM). It may be considered a semi-parametric approach because it provides polynomial functions of time for approximate growth curves.
Even though GBTM was introduced to describe the criminal behavior, wider fields of application have been emerging: evolution of salaries over time (Guigou et al., 2012), clinical research and causal analysis (Haviland et al., 2008). Moreover, despite being originally developed to model individual trajectories, GBTM has been applied to study phenomena in other units of analysis such as areal units (Griffiths and Chavez, 2004;LaFree, Morris, Dugan, 2010;Dugan and Yang, 2012).
From the perspective of clustering time series, it can be observed how, logically speaking, time series data look like panel data, although panels have a larger number of units n than time points T. Moreover, most panel data models are designed to deal with a limited number of repeated observations. It therefore follows that in this analysis, the use of GBTM is a reasonable compromise and a promising tool even when compared with more flexible approaches such as B-splines, because it provides a more synthetic representation of the time series. The next paragraph describes the GBTM approach in more detail.

Methodology: Group Based Trajectory Modelling (GBTM)
In this paragraph, the GBTM approach is described as a special case of LCGA and GMM. A basic functional form of the growth process for a GMM may be: where t=0, …, T is time (the second column of the design matrix), i indicates the single unit, and k indicates the group. The specification of a quadratic term is recommended to account for any non-linear change. In the GMM model, the parameter betas are the latent growth factors and are described as a function of their means in vector  and residual components in vector : The covariance matrix of the latent factors may differ among the groups, and the disturbances in (1) are homoscedastic. The model is completed by an assignment function, which is generally multinomial logit.
GBTM allows no variations in the growth factors. Therefore, equation (2) is absent. It follows that, under LCGA, a trajectory group is defined as a group of units that follow an identical developmental pathway where only random errors account for each individual's deviation from the respective group average. Population variability is captured essentially by differences across groups in terms of shape and level of trajectories. Conversely, for different groups of units in GMM, growth trajectories are allowed to vary around different means (Muthén et al., 2000).
In the analysis, we followed the GBTM specification which leads to a special case of a LCGA model. Let Yi:{yi0, yi1,…, yiT} denote a time sequence of measurements on the unit ith over T+1 periods. These measurements do not necessarily pertain to an individual (or to individual behaviour) but they may refer to a community, or can measure quantity as a poverty rate, or a crime rate, etc. (Nagin, 2005; p. 24).
Let P(Yi) with i=0,…,T, denote the joint probability (function or density) of Yi; i.e. the unconditional probability of observing the unit's time sequence of T+1 measurements. In the case of population mixture we have: Where P k (Yi) is the probability of Yi given membership in group k, where k=1,…,K, and K is the number of groups (K finite). Expression (3) describes a finite mixture model, because the sum works across a finite number of groups that make up the total population. For a given k, conditional independence is assumed for the sequence of the realizations yit of Yi over the T+1 periods: where p k (yit) is the probability (density) function of yit given membership in group k. The assumption of conditional independence implies that for each individual within a given trajectory group k, the distribution of yit at time t is independent from the realised level of the outcome in past periods t1, t2, etc. This assumption reduces the complexity of the model, but it could be questionable even if, for the sake of simplicity, many applications of finite mixture modeling with longitudinal data assume conditional independence (Nagin, 2005;p. 26). In other words, given membership in group k, the unit-level deviations from the group trajectory are assumed to be serially uncorrelated.
A conditional independence assumption is present in the standard GMM because the sequential realizations yit are independent, conditionally upon the individual's random effect. The difference is that in GBTM, serial uncorrelation is invoked at the group-level; in GMM at the unit-level. It means that this assumption is stronger in GBTM than in GMM (Nagin, 2005;p. 27).
The likelihood to be maximised over the n units is: The estimation procedure assumes multinormal variables and for the application carried out in this paper, the model allows for up to a quadratic relationship between ytik and time t: where tik is a disturbance assumed to be normally distributed with a zero mean and a constant variance  2 and the betas are the parameters determining the shape of the group trajectory.
The final component of the model is the specification of the group probabilities k: The requirements that across all K groups the 's sum to 1 is achieved with the estimation of K1 parameters, because the probability of one group is computed as one minus the sum of the probabilities of the other K1 groups.
The assignment to groups is based on the estimated posterior probabilities of group membership, according to the maximum probability rule. The estimated posterior probability for a unit i in group k is denoted by ) Y | k ( P i and Bayes' theorem provides the analytical basis for calculating this: In order to discriminate between models, we can distinguish two types of diagnostics (Nagin, 2005).
1) The accuracy of the group membership classification is based on the maximum posterior probability rule and the following index is used (Nagin, 2005; p. 88): is the posterior probability of the unit i in group k ,after its assignment to group k , and nk is the number of units assigned to group k. The numerator of this ratio is the odd of a correct classification in group k on the basis of the maximum posterior probability classification rule. The denominator is the odd of correct classification based on random assignment. Simulation analyses showed that values OCCk>5 indicate a high assignment accuracy (Nagin, 2005;p. 89).
2) The values of the Akaike information criterion (AIC) and Bayesian information criterion (BIC) that show whether the addition of a new group is statistically significant.
In addition, we also examined the trajectories identified by the procedure, and the transition in the grouping steps.
One disadvantage of GBTM compared to GMM is that due to not allowing variations about the group means, more groups are usually required to explain total variability. In fact, while under GBTM, population variability is captured by differences across groups, GMM adds an additional source of variability at the individual level. For these reasons, an adequate number of groups should be provided in order to encompass the absence of a model parameter accounting for individual variability. However, a recent survey on GBTM applications concluded that the number of units does not affect the number of trajectories identified when the size is 250 or over, and that the majority of studies identified 3 to 5 trajectories (van Dulmen et al., 2009) Tuscany and tourism Tuscany is the Italian region with the third highest number of arrivals (ISTAT, 2013). In Italy, more than 60% of the arrivals and presences of foreigners are mainly concentrated in four regions (Veneto, Trentino-Alto Adige, Tuscany and Lazio), while the figures for residents show a lower concentration (Emilia-Romagna, Veneto, Tuscany and Trentino-Alto Adige are less than 45%). The average number of overnight stays is close to the average values for Italy.
In a study dated 1996 (Formica and Uysal, 1996), the lifecycle of Italy as a tourist destination was in the growing stages in the years 1950-1973; between 1974-1987 it reached maturity, whereas during the years 1988-1996 the beginning of a decline was observed. However, Formica and Uysal (1996) recognised that between North and South there are remarkable spatial variations in forms and levels of accommodation depending on the tourist products. Moreover, different tourist products in the same destination may be at a different stage of the lifecycle.
In Figure 1, we can see a common growth of tourist arrivals for Italy, and Tuscany, although lower figures can be observed for Tuscany. Specifically, in the 2000-2013 period, an overall change of 30% occurred for Italy, while it was 20% for Tuscany. The slowdown of arrivals in Italy which occurred in the late nineties and early years of the new century, picked up in the following years with a definite revitalization. During the same period, Tuscany partially lost its market share: in 2000 there were12.5% of total arrivals in Italy, with 11.7% in 2013. The  Source: Authors' computation using data from the Tuscany Region Statistical Office figure reveals that Tuscany seems more sensitive to the economic cycle, despite remaining one of the leading tourist destinations, thanks especially to its uniquely attractive resources.
An important feature of tourism in Tuscany is its variety, also leading, among its benefits, to a mitigation of seasonality: a major problem for sustainability. In addition to seaside and mountain tourism, the most important resource is represented by cultural tourism, and thermal and rural tourism are also significant. Figure 2 shows the types of tourism per Tuscan municipality; we have chosen the main type for each municipality, even though various  historical cities are classified in the "Art/business" group, although other tourist offers are also important (for example Pisa, the territory of which also comprises important seaside resorts; coastal areas that can offer rural tourism as well);  several municipalities contiguous to historical cities are classified in the "Other" group, even though they are contiguous to historical cities (for example, the municipalities contiguous to Florence).
Due to the extremely high variety of the tourism offer, the time pattern of arrivals in Figure 1 is determined by the evolution of the different tourism products, as illustrated in Figure  2. Figure 3 shows the time series of arrivals in the years 2001-2013 (expressed as index numbers with base year 2000), according to the types of tourism. We can observe the different time patterns of "Rural tourism" and "Others" which show a major growth. This can be explained by the development of agritourism occurring over the last twenty years, and the rise in interest in tourism among destinations not traditionally devoted to tourism activities. Figure 3 also shows the dramatic stagnation of "Mountain" and "Thermal" products and the positive performance of cultural and business tourism (that comprises the tourist flows of almost all the provincial capitals).

GBTM: data and description of the clustering process
This paragraph describes the results of the GBTM applied to the time series. The data set provided by the Statistics Office of the Tuscany Region contains the greatest territorial details for data on tourist flows. This paper analyses the total number of annual arrivals in the years 2000-2013, with no distinction between nationality and the types of accommodation facilities.
Although bed-nights better characterise economic returns of tourism (see for example, Barros et al., 2011), this variable is influenced by the type of tourism and its effect is far more remarkable at a local level. For this reason, we decided to use tourist arrivals as a more suitable measure to represent the destination performance and allow for a correct comparative analysis among municipalities.
For confidential issues, data are not released for all of the 287 municipalities of Tuscany. Additional aggregations are performed in order to construct a complete time series with positive numbers. A total of 262 areal units are used in the application, resulting from the The comparison between tourist arrivals over time is influenced by scale and magnitude effects which relate to the extension of the territories and development of tourist activities. Because we are interested in identifying the trajectories (i.e. the shape of the time patterns), the time series is transformed as suggested in Gursoy et al. (2011) in the following manner: instead of the annual change rate: where yit is the number of arrivals in the area I and year t, mi is the annual average of arrivals in the area i. Expression (9) is a raw index of performance, which expresses the time pattern with respect to the average of the period (in this case, [2000][2001][2002][2003][2004][2005][2006][2007][2008][2009][2010][2011][2012][2013]. Due to the remarkable heterogeneity of the size of the areas in terms of tourist arrivals, the use of rit tends to generate exceptional values (Giusti and Grassini, 2014).
Discovering and understanding similarities in tourist arrivals across municipalities and over time can shed light on the different developmental phases of the various destinations. Therefore, we can expect to detect at least the following ideal time patterns: (i) a stable trend, denoting mature destinations that maintain the same tourist flows, contrasting the crisis of 2008; (ii) a weakly increasing trend (also typical of mature destinations); (iii) a weakly decreasing trend, for mature destinations in a declining phase and/or unable to contrast the 2008 crisis; (iv) an increasing trend, for destinations that were in the developing phase during the years 2000-2013 (for example the destinations with agritourism activities).
After classifying the areal units into trajectory groups, a test is performed to verify whether spatial contiguity between areas of the same group is statistically significant. As in Gursoy et al. (2013), joint count statistics are used, where a "join" is understood as a connection between areas according to a queen contiguity grid. The method used compares the number of joins (contiguous territories in the same cluster) with the number that would occur in the case of Source: Authors' computation using data from the Tuscany Region Statistical Office randomness. This test was carried out through the spdep package from R (Bivand, 2012).
The analysis presented in this section describes the steps from a 2-group partition up to a 6-group partition which constitutes the final result. For each k-group partition, statistics are reported in Table 1 and each transition phase is shown from k1 to k group (Figures 5-8). Finally, the equations of the six group trajectories and the assignment function are illustrated in Table 2.
The figures 5-8 describe the transition from the k to k+1 groups. The transition from 2 to 3 isolates an increasing pattern (group 3/3), and reveals that the shape of group 1/2 is determined by two different paths: one increasing (1/3) and the other decreasing (2/3). The transition from 3 to 4 helps specify more clearly the nonlinear growing trend, whereas the decreasing one is already identified with the  Finally, in the 6-group partition, group 5/5, apparently with an increasing, concave trend, was broken down into two new groups: (2/6)  Source: Authors' computation using data from the Tuscany Region Statistical Office with a declining concave shape, and (5/6) with a weakly increasing trend. Table 1 shows some diagnostics of the grouping process. While on one hand the 6-group configuration confirms the presence of an almost stationary trend (3/6 and 2/5) and two exceptional trajectories (4/6 and 3/5) already detected in the 5-group partition, on the other, it reveals that trajectory (5/5) is actually a Group 1/6: rapid increase Group 2/6: slowdown (   Source: Authors' computation using data from the Tuscany Region Statistical Office mixture of two different shapes: (2/6) and (5/6). Despite the lowest assignment probabilities, the other statistics in Table 1 show a good performance of the final 6-group configuration.
The average assignment probabilities are in any case high and therefore, also the values of OCCk. Table 3 shows the comparison with the results of the k-means algorithm. The main differences are found in the GBTM groups 3/6 and 4/6, but few units are involved. Finally, Figure 9 displays the six group trajectories with the individual trajectories (grey lines).

Discussion of the results
The interpretation of the six trajectory groups is carried out by considering the types of tourism (Table 4), the time evolution of the number of bed places (Figure 10), and the geographical distribution ( Figure 9). In particular, Table 4 illustrates the cross-tabulation of each tourism type per trajectory group. As described above, the original dataset assigns a single type to each of the 287 municipalities. The assigning of the types to the 21 new areas was carried out on the basis of the component municipalities.
Group 1/6 comprises marginal areas most of which are rural destinations (5 over 7) with a rapid growth and also a 20% increase in the number of beds (4% in hotels, 116% in other accommodation facilities). One municipality (Bagno a Ripoli), contiguous to both Florence and the Chianti area, belongs to this group.
Group 2/6, consists of 35 municipalities mostly located in rural and mountain areas of the south (in the proximity of Mount Amiata and the Colline metallifere). It is historically less industrialised also due to its geomorphological characteristics. Although this group has witnessed an increase in the number of beds ( Figure 10) equal to 117% (55% hotels, 152% non-hotels), its attractiveness has been diminishing. The join-test is significant, confirming a proximity effect.
Group 3/6 is very large (96 areas) and comprises very heterogeneous municipalities: coastal zones, industrialised areas, important tourist destinations. On the whole, the group consists of mature destinations which recorded a stationarity pattern in the number of arrivals and a weaker growth (13%) in the number of beds (4% hotels, 20% non-hotels). The most famous tourist destinations of Tuscany are included in this group: Siena, Arezzo, searesorts like Viareggio, the Island of Elba, and the historical spa resorts of Montecatini Terme and Chianciano Terme. The join-test is not significant throughout the whole group, however we can observe that in actual fact, some local proximity effects seem to be present in the eastern (province of Arezzo) and northern parts of the region (province of Lucca). In fact, although some isolated areas are present, the figure shows larger coloured areas determined by contiguous municipalities.
Group 4/6 comprises marginal units (only 13) mainly located in mountain areas with a decline both in the number of arrivals and number of beds. The decrease in the number of beds is 3% due to a decline in the non-hotel establishments (8%). The join-test is significant, showing a proximity effect for which the far northern part of Tuscany is responsible (Figure 9).
Group 5/6 is quite large (80 areas) and comprises municipalities with a noteworthy tourist industry and a relative major occurrence (with respect to the whole of Tuscany) of destinations devoted to art/business and rural tourism. These areas were able to counteract the crisis and reinforce the tourism offer with a slight increase in the number of beds, especially in the non-hotel establishments (34% hotels, 62% non-hotels).The destinations visited by greatest number of foreign tourists, like Florence, Pisa, Lucca, San Gimignano, and the Chianti area, are included in this group.
Other municipalities with a special interest in tourism are Livorno (with an important harbour), Montepulciano and Montalcino (with renowned cellars and wine production). This cluster includes 63% of coastal destinations: the majority are areas of the central coast of Tuscany (province of Livorno). It is interesting to note that even contiguous municipalities of these coastal areas are also included. The jointest is not significant as the areas are numerous and dispersed over the territory, however some linked areas can be observed (larger coloured areas). In particular, one on the coast, one in the centre of Tuscany (Chianti), and one in the south-east (in the proximity of Montepulciano, in the Province of Siena).
Finally, group 6/6 consists of internal areas devoted to rural or thermal tourism. The increase of arrivals and the number of beds denotes a strategy for (and the possibility of) enlarging the offer (in terms of number of beds: 115% hotels, 216% non-hotels).
The shapes of the trajectories are consistent with the evolution of the size of the tourist offer expressed by the number of beds offered. Group 3/6 and 4/6 show a decline, but also a stability in the number of beds, indicating that a threshold has been reached.
A more in-depth analysis of the types of tourism (Table 4 and Figure 9) reveals that most art/business destinations share a stationary or increasing pattern, confirming that this tourism offer is healthy. Mountain products are characterised by a contrasting situation. The pattern of arrivals for rural tourism is similar to that of art/business. This tourism segment is an important component of the tourism industry and it serves tourists interested not only in rural or natural landscapes, but also in art/business or sea. The traditional areas with ski resorts appear to be less attractive than other areas (southern), where mountain and rural landscapes are combined. Coastal tourism is definitely stable (group 3/6) with the presence of growing areas in the province of Livorno (group 5/6) however.
Finally, a few remarks must be made about thermal tourism. The results of the GBTM approach confirm the weakness of the traditional spas characterised by a concentration of large accommodation structures dating back to the 60's/70's (group 3/6) and more modern facilities with a mixed offer of spas and agritourism, typical of the small thermal municipalities (group 6/6). Ultimately, areas with no specific tourist offer ("other" types) still seem to be able to promote themselves as tourist destinations with an increase in the number of beds and number of arrivals.
In conclusion, the study confirms the presence of heterogeneity in the growth patterns of different areal units, and reveals three main findings. The first concerns the shapes of the trajectories identified. Three types of trends (stationary, increasing, decreasing) are detected, which, despite a short time series, can be linked with the different stages of the destination lifecycles (Butler, 1980). Stationary and regular increasing patterns (groups 3/6 and 5/6) occur for the majority of destinations (176 out of 262). This result was expected for a region like Tuscany, where tourism is at the mature stage and has played a remarkable role in the local economies. The presence of a decreasing trend in the period 2000-2013 can be attributed to a weakness of some destinations to counteract the global financial crisis and/or the stage of decline. The occurrence of increasing trends confirms, on one hand, the attractiveness of typical destinations of Tuscany (group 5/6) that also show an increase in the number of beds, while on the other, it reveals the presence of destinations (although few and marginal) which are at the very early stages of their lifecycle, as can be observed in the remarkable growth in the number of beds (group 1, Figure 10).
A second finding concerns the association between trajectory groups and tourism typology. There is a general stability of the classical tourism destinations (business/art and urban segments), which are renowned worldwide. The analysis also shed light on the development of agritourism activities (taking place during the period under study), most of which are included in group 5/6. In the case of the thermal offer, a contrasting situation emerges. Traditional thermal stations, aimed at attracting mass tourism, were not able to renovate their properties and facilities, or relaunch their offer, and as a result they show a slowdown or even a decrease. Smaller destinations with a mix of thermal and agritourism offers have been developing due to exploiting the demand for rural tourism, interest in nature, environment and landscape (e.g. thermal destinations in the province of Siena or Grosseto).
The third finding concerns the effects of geographical proximity. This is only significant in the negative situation (slowdown and decline patterns: groups 2/6 and 4/6), even though clusters of areas seem to emerge from the maps. In any case, the maps of groups 3/6 and 5/6 seem to show the presence of spatial homogeneity: along the coast for group 3/6 and in the south-eastern area (province of Arezzo) for group 5/6. This spatial homogeneity might be also determined by the joint projects and agreements among municipalities. In fact, there are some organizations which manage the planning and marketing of the local territories (for example: the Chianti area, Versilia, etc.) despite the fact that the power and resources for implementing actions are mainly determined by administrative institutions. This can be a limitation as political barriers may fail to take into consideration consumer preferences or tourism industry functions (Buhalis, 2000).

Conclusions
The analysis developed in this paper addresses the dynamics of spatial tourist flows (arrivals) within and across municipalities of Tuscany over time. Despite the relative shortness of the time series, it has been possible to draw conclusions even about the development stage of some destinations which do not conform with the overall trend of the region. Two original features characterise this work: 1) Unlike most works, which analyse countries, regions or provinces, the municipality is the unit of analysis used here. Although it is an administrative area, it constitutes the basic elementary unit for studying the territory, forming networks of interdependence within a geographic proximity, and designing development strategies. 2) As far as the statistical methodology is concerned, the trajectories are derived from the GBTM approach with a univariate analysis of the time series of arrivals with no covariates outside the time variable. This is due to the interest in an exploratory approach. In this respect, the GBTM approach, compared with other exploratory ones, such as the k-means clustering, has the advantage of providing an algebraic interpretation of the trajectory through the estimated parameters of the functions representing the trajectories. Despite that fact that each trajectory group is characterised by a parametric function of time, we understood trajectory groups as providing approximations for a more complex reality, and tried to move away from their interpretation as distinct units.
The use of GBTM with small areal units called for the definition of a relatively large number of groups in order to account for any exceptional observations, also due to the area size in terms of number of arrivals. The GBTM methodology is combined with a test of spatial contiguity that is significant in the presence of a downward trend or slowdown pattern, though some local proximity effects seem to be present in the case of weak increases. The combination of two graphical products (trajectory and territorial map of the groups) enables an intuitive interpretation of the regional tourism development. Moreover, the fact that destinations with the same tourism offer (for example, thermal) are included in different trajectory groups, can allow for comparisons and benchmarking exercises.
As with other works, this study is not free of limitations. The empirical analysis is only based on the number of arrivals, while the performance of the tourism sector should be assessed with a multifaceted approach. In this respect, further development should be devoted primarily in the construction of a richer dataset at the municipality level, also with information about political variables such as the effective presence of forms of cooperation among territories, in order to include additional covariates in the trajectory formula and go beyond the simple exploratory approach applied in this paper. For the Tuscany region, this issue is not far from being accomplished, because the TDO experience is going to provide a repertory of comparable information at the municipality level.
Endnote: 1 The graphs in Figures 5-8 display the values of dit of expression (6). The minimum and maximum values of the vertical axis are 1 and 2 respectively. Time (years 2000-2013) is the variable on the horizontal axis.