ISSN 1440-771X

Department of Econometrics and Business Statistics

http://www.buseco.monash.edu.au/depts/ebs/pubs/wpapers/

The tourism forecasting competition

George Athanasopoulos, Rob J Hyndman, Haiyan Song, Doris C Wu

October 2009

Working Paper 10/08 Revised Oct 2009

The tourism forecasting competition

George Athanasopoulos* Department of Econometrics and Business Statistics and Tourism Research Unit, Monash University, VIC 3800, Australia. Email: George.Athanasopoulos@buseco.monash.edu.au Rob J Hyndman Department of Econometrics and Business Statistics, Monash University, VIC 3800, Australia. Email: Rob.Hyndman@buseco. monash.edu.au Haiyan Song The School of Hotel and Tourism Management, The Hong Kong Polytechnic University, Hong Kong. Email: Haiyan.Song@inet.polyu.edu.hk Doris C Wu The School of Hotel and Tourism Management, The Hong Kong Polytechnic University, Hong Kong. Email: Doris.Wu@inet.polyu.edu.hk

*Corresponding author

7 October 2009

JEL classication: C22,C52,C53

The tourism forecasting competition

Abstract:

We evaluate the performance of various methods for forecasting tourism demand. The data used include 366 monthly series, 427 quarterly series and 518 yearly series, all supplied to us by tourism bodies or by academics who had used them in previous tourism forecasting studies. The forecasting methods implemented in the competition are univariate and multivariate time series approaches, and econometric models. This forecasting competition differs from previous competitions in several ways: (i) we concentrate only on tourism demand data; (ii) we include approaches with explanatory variables; (iii) we evaluate the forecast interval coverage as well as point forecast accuracy; (iv) we observe the effect of temporal aggregation on forecasting accuracy; and (v) we consider the mean absolute scaled error as an alternative forecasting accuracy measure.

Keywords: Tourism forecasting, ARIMA, exponential smoothing, time varying parameter model, dynamic

regression, autoregressive distributed lag model, vector autoregressions.

The tourism forecasting competition

1 Introduction
Over the past four decades, tourism has developed into one of the most rapidly growing global industries. The World Tourism Organization (2008) reports that international tourist arrivals world-wide grew at a rate of 6% in 2007, reaching nearly 900 million, compared to 800 million two years earlier. Both academic interest and the tourism literature have grown, parallel to this growth in the industry, producing many articles that model and forecast tourism ows between various countries. These articles vary immensely in scope, modelling and forecasting techniques, and data types, lengths and frequencies. The three major literature review articles that attempt to summarise these are Witt and Witt (1995), Li et al. (2005) and Song and Li (2008). Despite the authors best efforts, the diversity of the studies has not led to a consensus about the relative forecasting performance of commonly used methods when they are applied to tourism data. In this paper we apply forecasting methods to a very broad collection of series within the eld of tourism. This allows us to draw conclusions for methods within this eld, and also to contribute some general observations and conclusions relevant to the broader eld of forecasting. Since the last of the M series of forecasting competitions was published (the M3 competition, see Makridakis and Hibon, 2000), there have been no more major contributions to this line of research. In this paper we take on the challenge of creating a forecasting competition that overcomes some of the limitations highlighted by the commentators of the M3 competition (see Ord, 2001). For example, in the M3 competition, 3003 series from various business and economic sectors were used. This was larger than the number of series used in any previous study. However, the question of whether bigger means better was raised by several commentators. Suggestions were made calling for mini-competitions along the lines suggested by Fildes and Ord (2004), who proposed the use of a more homogeneous set of series. In an effort to use a homogeneous set of series which are representative of some population, we use 366 monthly, 427 quarterly and 518 yearly series, all from the eld of tourism. We present the data in Section 4. From such a set of series we intend to draw some general inferences on the modelling and forecasting of tourism demand. We will also examine whether the results from previous competitions carry over to this well-dened set of series. According to the M3 competition results, one of the most accurate forecasting methods was Forecast Pro (see Goodrich, 2000). This method was arguably the most accurate for seasonal data and was second only to the Theta method (Assimakopoulos and Nikolopoulos, 2000) for non-seasonal data. The overall disappointing performance of ARIMA-based methods led the authors to conclude that statistically sophisticated methods do not necessarily produce more accurate forecasts than simpler ones. Like some of those who commented on the M3, we challenge this conclusion in the context of forecasting tourism demand. We evaluate forecasts from both Forecast Pro and the Theta method, as well as forecasts from two newly proposed fully automated forecasting algorithms, one which identies and estimates

The tourism forecasting competition

ARIMA models and another which identies and estimates state space models that underly exponential smoothing methods. We present these methods in Section 2 and evaluate their forecasting performance in Section 5.1. One of the issues we address is the effect of temporal aggregation on forecast accuracy. In Section 5.2, the monthly series are aggregated to be quarterly, and the quarterly series are aggregated to be yearly. Therefore, we can directly compare the accuracy of the forecasts made before and after aggregation. The importance of producing statements of uncertainty when forecasting has long been undervalued in the empirical forecasting literature. There is a gap in the literature between the methodological developments in density forecasting (see, for example, Tay and Wallis, 2002; West, 2006) and the applied forecasting papers that do not produce statements of uncertainty to accompany their point forecasts. For instance, a signicant limitation of all forecasting competitions to date has been the lack of assessment of whether the forecasting methods can produce reliable forecast intervals. This is highlighted in the commentaries on the M3 competition (see Armstrong, 2001; Goodrich, 2001; Koehler, 2001; Tashman, 2001). Not observing the uncertainty associated with a point forecast can lead to a false sense of accuracy. Economic planning based on forecast intervals can be very different to that based on mean or median forecasts. In Section 5.3 we attempt to ll this gap in the literature by also evaluating the forecast coverage probabilities of the forecasting methods. Using models with exogenous variables for policy analysis and forecasting is common in both the tourism literature and the tourism industry. These are usually labeled as causal approaches, and take on various functional forms (see Song and Witt, 2000, for a detailed exposition on their use within the tourism literature). A major challenge in forecasting with this type of model is that once a satisfactory relationship has been specied, the user needs to produce forecasts of the exogenous variables to be able to forecast the variable of interest. For the rst time in a large forecasting competition, we consider the forecasting performance of such models. There are 98 quarterly and 128 yearly cases for which we have available the explanatory variables that are typically used in the literature. An alternative approach is to treat all of the variables as endogenous (Sims, 1980). This leads us to also consider multivariate models in our forecasting competition, as was suggested by Granger (2001) and Ord (2001). We present the econometric approaches and the methodologies we implement in Section 3. The forecast performance of these models is evaluated in Section 5.4.

1.1 Literature review, rules and objectives

As we mentioned in the introduction, attempting to draw general conclusions from the existing tourism literature through survey articles is very difcult. Papers vary in scope, data frequencies, modelling

The tourism forecasting competition

frameworks and estimation techniques, forecast evaluation rules (ex-post versus ex-ante for causal models), and forecast horizons, and also in results and conclusions. In this section we establish various rules that will be followed in order to achieve a degree of uniformity and fairness across the application and evaluation of the forecasting methods. We also set some objectives and list the questions we ask in this research. Outline of rules: We aim to apply general modelling frameworks with objective and clearly stated decision rules. We aim to replicate, and hence evaluate, some of the typical modelling procedures both within and outside the tourism literature. All out-of-sample values are used only for evaluating the forecasts generated by the competing methods, never in the modelling stages. No intervention is applied to the pure time series methods in terms of one-off events causing possible structural breaks. However, all such exogenous information is included in the causal approaches. All models are estimated once only, and only one set of forecasts is produced for each series. All forecasts of exogenous variables are ex-ante unless otherwise specied. We should note that all forecasts from the causal approaches were generated by authors Song and Wu, while all forecasts from the time series methods were produced by authors Athanasopoulos and Hyndman. (See the acknowledgement in Section 7 for the Theta method forecasts.) Summary of objectives: The rapid increase in the capacity of computers to store information has generated an abundance of data across all types of industries. For example, Athanasopoulos et al. (2009) generated forecasts for 117 tourism demand series (and that only includes Australian domestic tourism), disaggregated only by selected geographical areas. In total, Tourism Australia generates forecasts for thousands of series every quarter when considering inbound, outbound and domestic travel, as well as numerous levels of disaggregation such as geographical regions, purpose of travel and so on. Hence, accurate automatic forecasting procedures have become a necessity in order to take advantage of such a plethora of information. We evaluate the performance of three fully automated algorithms (with no intervention). ARIMA models have not proven as accurate as other forecasting methods, whether the model identication is automated or implemented manually (refer to Fildes and Ord, 2004, and references therein). In fact, Armstrong (2006) lists these as tested areas with little gain in accuracy. We re-evaluate the forecasting performance of ARIMA models using a recently proposed algorithm that has shown promising performance in smaller scale evaluations.

The tourism forecasting competition

Fildes and Ord (2004) and Armstrong (2006) highlight the dominance of the damped trend method in previous forecasting competitions. We evaluate what is gained or lost by considering aggregate model selection procedures instead of implementing specic methods, considering forecasts from both the Theta method and the damped trend method. One of the ndings of Witt and Witt (1995) is that, for annual data, the Nave method seems to produce the most accurate forecasts (especially for one year ahead). We revisit this result and also examine whether we can aid the performance of forecasting methods by using higher frequency data. Previous studies have found the forecast intervals obtained to be too narrow; hence, the actual values fall outside the empirical forecast intervals more often than they should (see for example Makridakis and Winkler, 1989; Chateld, 2001; Hyndman et al., 2002). Makridakis et al. (1987) conclude that the lower the frequency of the data, the more the coverage probability of the forecast intervals is over-estimated. We re-examine this conclusion by evaluating the forecast interval coverages for the three automated forecasting algorithms for monthly, quarterly and yearly data. Allen and Fildes (2001) (who collated the results from Fildes 1985 and Armstrong 1985) found that models with exogenous variables forecast better than extrapolating methods when ex-post forecasts are used for the regressors. A surprising result from their study is that the forecasting performance of causal models seems to improve when using ex-ante, rather than ex-post, forecasts. In the tourism literature, Song et al. (2003a) found that econometric models perform better than the no-change, ARIMA and VAR models, using ex-post forecasts. In contrast, Witt and Witt (1995) concluded that causal models are outperformed by the no-change model, regardless of whether ex-ante or ex-post forecasts are used. We evaluate the performance of models with explanatory variables, as implemented in the tourism literature, in a forecasting competition setting in which all out-of-sample values are ex-ante forecasts. We then revisit this result with ex-post forecasts for the regressors in order to evaluate the best possible result from using causal models in a scenario-based forecasting framework. Despite strong warnings about its limitations (see Hyndman and Koehler, 2006), the MAPE remains the most commonly used forecast error measure among both academics and practitioners (see Fildes and Goodwin, 2007), and the tourism forecasting literature is no exception (see Li et al., 2005; Song and Li, 2008). Here we investigate in detail the forecasting results based on the MAPE and the effect that some of its limitations may have on the results. We also consider the MASE (mean absolute scaled error), which was proposed by Hyndman and Koehler (2006) in order to overcome some of the limitations of the MAPE.

The tourism forecasting competition

2 Pure time series approaches

In this section we present the details of the three fully automated forecasting procedures and the two method-specic procedures that we have implemented in this paper.

2.1 ARIMA forecasting

A non-seasonal ARIMA(p, d, q) process is given by (B)(1 B d ) y t = c + (B) t , where { t } is a white noise process with mean zero and variance 2 , B is the backshift operator, and (z) and (z) are polynomials of orders p and q respectively. To ensure causality and invertibility, it is assumed that (z) and (z) have no roots for |z| < 1 (Brockwell and Davis, 1991). If c = 0, there is an implied polynomial of order d in the forecast function. The seasonal ARIMA(p, d, q)(P, D, Q)m process is given by (B m )(B)(1 B m ) D (1 B)d y t = c + (B m ) (B) t , where (z) and (z) are polynomials of orders P and Q respectively, each containing no roots inside the unit circle. If c = 0, there is an implied polynomial of order d + D in the forecast function (Box et al., 2008, pp. 381382). The main task in automatic ARIMA forecasting is selecting an appropriate model order, that is, the values of p, q, P, Q, D and d. We use the automatic model selection algorithm that was proposed by Hyndman and Khandakar (2008), which is summarised below. Diebold and Kilian (2000) nd strong evidence that unit root pretesting for selecting the level of differencing d improves the forecasting accuracy. For non-seasonal data (we treat yearly data as non-seasonal) we consider ARIMA(p, d, q) models, where d is selected based on successive KPSS unit-root tests (Kwiatkowski et al., 1992). That is, we test the data for a unit root; if the null hypothesis of no unit root is rejected (at the 5% signicance level), we test the differenced data for a unit root; and so on. We stop this procedure the rst time we fail to reject the null hypothesis. For seasonal data we consider ARIMA(p, d, q)(P, D, Q)m models, where m is the seasonal frequency. Unlike Hyndman and Khandakar (2008), we set D = 1 for all seasonal data (we treat all monthly and quarterly data as seasonal), as we nd that their suggested seasonal unit root test does not help in selecting the appropriate order of differencing (Osborn et al., 1999, reached a similar conclusion for seasonal unit root

The tourism forecasting competition

testing). We then choose d by applying successive KPSS unit root tests to the seasonally differenced data. Once the value of d has been selected, we proceed to select the values of p, q, P and Q by minimizing the AIC. We allow c = 0 for models where d + D < 2. Once d and D are known, we select the orders p, q, P and Q via Akaikes Information Criterion: AIC = 2 log(L) + 2(p + q + P + Q + k), where k = 2 if c = 0 and 1 otherwise (the other parameter being 2 ), and L is the maximized likelihood of the model tted to the differenced data (1 B m ) D (1 B)d y t . The likelihood of the full model for y t is not actually dened, and so the values of the AIC for different levels of differencing are not comparable. There are a large number of potential ARIMA modelstoo many to allow us to estimate every possible combination of p, q, P and Q. Instead, we need a way to efciently traverse the space of models in order to arrive at the model with the lowest AIC value. Hyndman and Khandakar (2008) proposed the following step-wise algorithm. Step 1: Try four possible models to start with: ARIMA(2, d, 2) if m = 1 and ARIMA(2, d, 2)(1, D, 1) if m > 1. ARIMA(0, d, 0) if m = 1 and ARIMA(0, d, 0)(0, D, 0) if m > 1. ARIMA(1, d, 0) if m = 1 and ARIMA(1, d, 0)(1, D, 0) if m > 1. ARIMA(0, d, 1) if m = 1 and ARIMA(0, d, 1)(0, D, 1) if m > 1. If d + D 1, these models are tted with c = 0. Otherwise, set c = 0. Of these four models, select the one with the smallest AIC value. This is called the current model, and is denoted by ARIMA(p, d, q) if m = 1 or ARIMA(p, d, q)(P, D, Q)m if m > 1. Step 2: Consider up to thirteen variations on the current model: where one of p, q, P and Q is allowed to vary from the current model by 1; where p and q both vary from the current model by 1; where P and Q both vary from the current model by 1; where the constant c is included if the current model has c = 0 or excluded if the current model has c = 0. Whenever a model with a lower AIC is found, it becomes the new current model and the procedure is repeated. This process nishes when we cannot nd a model close to the current model with a lower AIC.

The tourism forecasting competition

There are several constraints on the tted models in order to avoid problems with convergence or near unit roots (see Hyndman and Khandakar, 2008, for details). The algorithm is guaranteed to return a valid model because the model space is nite and at least one of the starting models will be accepted (the model with no AR or MA parameters). The selected model is then used to produce forecasts and forecast intervals.

2.2 Innovations state space models for exponential smoothing

Ord et al. (1997), Hyndman et al. (2002) and Hyndman et al. (2005) (amongst others) have developed a statistical framework for the exponential smoothing methods presented in Table 1. The statistical framework incorporates stochastic models, likelihood calculations, forecast intervals and procedures for model selection. We employ this framework for building innovations state space models. The aforementioned papers have shown that these models generate optimal forecasts for all exponential smoothing methods (including non-linear methods). The classication of the exponential smoothing methods in Table 1 originated with Pegels (1969) and was expanded by Gardner (1985), Hyndman et al. (2002) and Taylor (2003). Each of the fteen methods listed has a trend and a seasonal component. Hence, cell (N,N) describes the simple exponential smoothing method, cell (A,N) Holts linear method, cell (A,A) Holt-Winters additive method, and so on. Table 1: Classication of exponential smoothing methods. Seasonal component A M (Additive) (Multiplicative) N,A N,M A,A A,M Ad ,A Ad ,M M,A M,M Md ,A Md ,M

Trend component N (None) A (Additive) Ad (Additive damped) M (Multiplicative) Md (Multiplicative damped)

N (None) N,N A,N Ad ,N M,N Md ,N

For each method, there are two possible state space models: one corresponding to a model with additive errors and the other to a model with multiplicative errors. Table 2 presents the fteen models with additive errors and their forecast functions. The multiplicative error models can be obtained by replacing t
t t

with

(for further details see Hyndman et al., 2008). We select models by minimising the AIC amongst all

models (both additive and multiplicative). We then compute forecast intervals from the selected models using analytic formulae (Hyndman et al., 2008, Chapter 6), or by simulation if the analytic formulae are not available.

The tourism forecasting competition

Table 2: State space equations for each additive error model in the classication. Multiplicative error models are obtained by replacing t with t t . In each case, t denotes the level of the series at time t, b t denotes the slope at time t, s t denotes the seasonal component of the series at time t, and m denotes the number of seasons in a year; , , and are constants with 0 < , , < 1 and 0 < < ; t+h|t denotes the h-step-ahead forecast based on all of the data up to time t; y h = + 2 + + h; t+h|t denotes a forecast of y t+h based on all the data up to time t; and y + hm = (h 1) mod m + 1.
Trend component N (none) t
t

N (none) = = = = = = = = = = = = = = = = = = =
t1 t1

Seasonal component A (additive) t

t t

M (multiplicative) t
t

st t+h|t y t A (additive)
t t

t+h|t y t
t t

= = = = = = = = = = = = = = = = = = = = = = = =

+ s tm + t s tm + t t + s tm+h+ m
t1 t1

st t+h|t y t
t

= = = = = = = = = = = = = = = = = = = = = = = =

t1 + t /s tm s tm + t / t1 t s tm+h+ m

t1 s tm

bt t+h|t y t

+ b t1 + b t1 + b t1 + t
t1 t1 t

bt st t+h|t y t
t t

+ hb t

+ b t1 + s tm + b t1 + t b t1 + t s tm + t t + hb t + s tm+h+ m
t1 t1

bt st t+h|t y t
t

+ b t1 )s tm + b t1 + t /s tm b t1 + t /s tm s tm + t /( t1 + b t1 ) ( t + hb t )s tm+h+ m
t1 t1

Ad (additive damped)

bt t+h|t y t

+ b t1 + b t1 + b t1 + t
t1 t1 t

bt st t+h|t y t
t

+ h b t

+ b t1 + s tm + b t1 + t b t1 + t s tm + t t + h b t + s tm+h+ m
t1 t1

bt st t+h|t y t
t

+ b t1 )s tm + b t1 + t /s tm b t1 + t /s tm s tm + t /( t1 + b t1 ) ( t + h b t )s tm+h+ m
t1 t1 t1 b t1 + t /s tm b t1 + t /(s tm t1 ) s tm + t /( t1 b t1 ) h t b t s tm+h+ m t1 b t1 s tm t1 b t1 + t /s tm b t1 + t /(s tm t1 ) t1 b t1 ) h t b t s tm+h+ m t1 b t1 s tm

M (multiplicative)

bt t+h|t y t

+ t b t1 + t / t1
t1 b t1 h t bt t1 b t1 t1 b t1 + t b t1 + t / t1 h t bt

t1 b t1

bt st t+h|t y t
t

+ s tm + t b t1 + t / t1 s tm + t h t b t + s tm+h+ m
t1 b t1 t1 b t1 t1 b t1 + s tm t1 b t1 + t b t1 + t / t1 t h t bt

bt st t+h|t y t
t

Md (multiplicative damped)

bt t+h|t y

bt st t+h|t y

s tm +

bt st t+h|t y

s tm + t /(

+ s tm+h+ m

We label this method ETS in the tables that follow.

The three letters are an abbreviation of

ExponenTial Smoothing, and also specify the three components of the stochastic model: Error, Trend and Seasonality. For example, an ETS(A,A,A) is a Holt-Winters additive method with an additive error component.

2.3 Forecast Pro

Forecasts from several commercial software packages were considered in the M3 competition. Forecast Pro was arguably the most accurate commercial package, as well as having the most consistent performance

10

The tourism forecasting competition

across all data. In this paper we evaluate forecasts from the Forecast Pro Extended Edition, Version 4. The forecasting method choice is set to expert selection. The software evaluates the forecasting performance of several methods and selects amongst them. Given the nature of the data we consider in this paper, the methods considered by the Forecast Pro algorithm were exponential smoothing, ARIMA models and simple moving averages. Although the ner details of the model identication, estimation and selection are not revealed, Goodrich (2000) presents some details. Exponential smoothing methods are tted by minimising the in-sample sum of squared errors. The nal method is selected by minimising the BIC, supplemented by some logical rules. With ARIMA models, a general-to-specic approach is followed. First a non-parsimonious state space model is estimated, and is used in turn to obtain approximate parameter estimates for a large number of potential ARIMA models. The nal model is selected by the BIC (again supplemented by some logical rules), and then re-estimated using unconditional least squares. This method is labeled ForePro in the tables that follow. For further details, refer to Goodrich (2000) or to
www.forecastpro.com.

2.4 Theta method

One method that performed extremely well in the M3 competition (Makridakis and Hibon, 2000) was the Theta method (Assimakopoulos and Nikolopoulos, 2000), which was further analysed and described by Hyndman and Billah (2003). For a given value of , a time series y t is transformed to x t, (dubbed a theta line) through x t, = a + b (t 1) + y t , Estimates of a and b are obtained by minimising
n t=1 [ y t

t = 1, . . . , n. x t, ]2 . As in the M3 competition, forecasts

are obtained by averaging two theta lines using = 0 (which gives a regression time trend) and = 2. The theta line for = 2 has been extrapolated using simple exponential smoothing for which the smoothing parameter has been chosen by minimising the in-sample one-step-ahead mean squared error, with the starting value for the initial level set equal to y1 . Hyndman and Billah (2003) show that in this case the forecasts obtained by the Theta method are equivalent to those generated by simple exponential smoothing with an added trend and a constant, where the slope of the trend is half that of a tted trend line through the original time series y t . All monthly and quarterly data are rst seasonally adjusted by extracting a seasonal component using classical multiplicative decomposition. The seasonal component is then added to the forecasts generated

11

The tourism forecasting competition

by the Theta method. The method is implemented using Delphi 7.0 for Windows XP The forecasting . software is TIFIS CM3, which is a non-commercial Forecasting Support System.

2.5 Damped trend

As was highlighted in the introduction, the damped trend method has been singled out from previous forecasting competitions as performing very well. In this paper we estimate the additive damped trend model ETS(A,Ad ,A) for monthly and quarterly data and ETS(A,Ad ,N) for yearly data, as presented in Table 2.

2.6 Nave approaches

We produce forecasts from two nave approaches which form natural benchmarks. For yearly data, y we use t+h|t = y t . Hence, all forecasts are equal to the most recent observation. This method is labeled Nave in the tables that follow. For monthly and quarterly data, we use t+h|t = y tm+hm , where y hm = [(h 1) mod m] + 1, with m = 4 for quarterly data and m = 12 for monthly data. Hence, all forecasts for seasonal data are equal to the most recent observation of the corresponding season. This method is labeled SNave, standing for Seasonal Nave, in the tables that follow.

3 Models with explanatory variables

The general tourism demand function in the tourism modelling and forecasting literature (e.g., Song and Witt, 2000) takes the form:
i y ti = f (g t , p i , p is , dummy variables, t t t ),

(1)

where y ti is the demand variable measured by tourist arrivals from origin country i to (or expenditure in)
i the destination country; g t is the income level of origin country i in real terms; p i represents the relative t

cost of living in the destination country for tourists from origin country i, measured as the relative CPI of the destination country to that of the origin country in constant prices, adjusted by the relevant exchange rates; p is represents tourism prices in substitute destinations, and is measured by a weighted average t price index of a set of alternative destinations to the destination country. For a detailed exposition on the price variables, refer to Wong et al. (2007) for a case study with Hong Kong as the destination country. The dummy variables include seasonal dummies and one-off events such as terrorist attacks, epidemics, or other events that relate to particular regions, and
t

is a random error term.

12

The tourism forecasting competition

The models we consider are special cases of

k pj p

y t = 0 +
j=1 i=0

j,i x j,ti +
i=1

i y ti +

t,

(2)

where y t is the tourism demand variable (now dropping the superscripts i and s for brevity); x j,t are the exogenous variables included in our model for j = 1, . . . , k; p is the maximum number of lags of the regressand and p j the maximum number of lags of the regressors; and assumed to be white noise.
t

is the error process, which is

3.1 Autoregressive distributed lag model (ADLM)

In the initial specication of the general ADLM given by equation (2), all possible variables are included. For quarterly data, p j = p = 4, and for annual data, p j = p = 1. This unrestricted specication of the ADLM is an error correction model which has been considered in many previous tourism studies (see Song and Witt, 2003, for further details). In the tables that follow, this specication is labeled as ADLM. The model is rened by implementing the dynamic econometric modelling technique known as the general-to-specic approach (as advocated by Hendry, 1986). The least signicant regressor (i.e., the one with the largest p-value) is deleted from the model, then the simplied model is re-estimated. This process is repeated until the coefcients of all of the remaining regressors are statistically signicant at the 5% signicance level (one-tailed). The nal model should be simple in structure and display no autocorrelation or heteroscedasticity, and preferably no non-normality either. In the tables that follow, this specication is labeled as ADLMR . Setting p = p j = 0, equation (2) gives a static regression specication which has been considered in many previous tourism forecasting studies (refer to the literature reviews by Li et al., 2005; Song and Li, 2008). We label this specication as SR in the tables that follow. As alternatives to this static specication, we also consider a static regression tted to the differenced data, which we label SR. In addition, for the quarterly data we consider a static regression tted to the seasonally differenced data, which we label as m SR. In our attempt to be as objective and general as possible within a static regression framework, one might argue that we are missing features of the data that could be captured and exploited through a more rigorous modelling framework, resulting in more accurate forecasts. In an effort to achieve this, we have also implemented the following regression framework, which we label DR (Dynamic Regression) in the tables that follow.

13

The tourism forecasting competition

Starting from a static regression in each case, we examine whether any unexploited dynamics remain by observing the estimated autocorrelation and partial autocorrelation functions of the residuals. We rst account for unit roots and seasonal unit roots by taking rst and/or seasonal differences of all of the variables, guided by the magnitude of the rst order and seasonal autocorrelations. Then, any dynamics left over in the residuals are modelled via ARMA specications, guided by the autocorrelation and partial autocorrelation functions. If the choice between AR and MA components becomes arbitrary (as either one or a combination of them completely captures all dynamics), then this choice is made by minimising the AIC.

3.2 Time varying parameter (TVP) model

The TVP model is used in order to allow the coefcient of the explanatory variables to change over time. This method is more adaptable when the assumption of constant coefcients is not valid, and structural changes in econometric models need tackling. The TVP approach uses a recursive estimation process in which the most recent information is weighted more heavily than the information obtained in the distant past. With the restriction p = 0 being imposed on the coefcients in equation (2), the TVP model can be expressed in state space form as: yt = t xt + (3) (4)

t = t1 + t ,

where x t is a (k + 1)-dimensional vector of the explanatory variables (including a column of ones for the intercept); t is a (k + 1)-dimensional vector of parameters and is known as the state vector;
t

NID(0, 2 ) refers to the temporary disturbance; and t NID(0, ) is the permanent disturbance.

The coefcient matrix is initially assumed to be known. Equation (3) is called the measurement or system equation, while equation (4) is known as the transition or state equation, which is used to simulate the way in which the parameters in the system equation evolve over time. If is an identity matrix, the transition equation becomes a random walk process:

t = t1 + t .

(5)

In most cases, the random walk process is adequate for capturing the parameter changes in various economic models (see, for example, Bohara and Sauer, 1992; Kim, 1993; Greenslade and Hall, 1996; Song and Witt, 2000). We adopt this approach for both quarterly and annual data. For quarterly data, seasonal dummies are also included to capture the seasonal component in the data. Equations (3) and (4)

14

The tourism forecasting competition

are estimated using the Kalman lter (for details of the estimation procedure, see Harvey, 1989). In the tables that follow, we label this specication as TVP . In order to forecast tourism demand using frameworks that incorporate exogenous variables, we must rst produce forecasts for each of the exogenous variables. As is typical in the tourism literature (see Song et al., 2003b; Song and Wong, 2003), we forecast these using exponential smoothing methods. For the quarterly data, we apply Holt-Winters additive method. For the yearly data, we have either aggregated the quarterly forecasts (where possible) or we apply Holts linear method.

3.3 Vector autoregressive (VAR) model

In contrast to the previous modelling frameworks, we now treat all of the variables in equation (1) as endogenous: yt y t1 y tp = 0 + 1 + + p + t , xt x t1 x tp (6)

where t NID(0, ). This framework was rst suggested by Sims (1980), and has been the main workhorse in multivariate economic modelling and forecasting ever since. It is also popular in the tourism literature (see Witt et al., 2003; Song and Witt, 2006). A great advantage of this approach in terms of the effort required for forecasting is that the system generates forecasts for all of the variables. Hence, we do not need to generate forecasts separately for the x variables. We consider three sets of forecasts from the following general approaches. 1. VAR models for all variables in levels. We choose the lag length of the model by minimising the AIC. We label these as VAR(AIC) in the tables that follow. The maximum lag lengths we consider are 4 for quarterly data and 2 for annual data. We also consider the largest of these models: VAR(4) for quarterly data and VAR(2) for annual data. This ensures that at least for quarterly data, the VAR models contain up to and including the seasonal lag. 2. Reduced VAR models for the growth rates of all variables. We rst calculate the growth rates of all variables by considering the rst differences of the natural logarithms. The growth rate of variable z t is calculated as 100 ln(z t /z t1 ). For quarterly data, we select the lag length by minimising the AIC with a maximum lag order of 4. For yearly data we set the lag length to 1, due to the short samples. All insignicant coefcients (at the 5% signicance level) are then restricted to zero. We impose the restrictions one at a time by eliminating the parameter with the lowest t-statistic at each step. The restricted VARs are estimated using the seemingly unrelated regression estimation method (Zellner, 1963), as not all equations include the same regressors. Athanasopoulos and Vahid (2008)

15

The tourism forecasting competition

found this type of VAR model to be successful in forecasting macroeconomic variables. We label them as VARR (AIC) in the tables that follow.

4 Data and forecast error measures

The data we use include 366 monthly series, 427 quarterly series and 518 yearly series. They were supplied by both tourism bodies (such as Tourism Australia, the Hong Kong Tourism Board and Tourism New Zealand) and various academics, who had used them in previous tourism forecasting studies (please refer to Section 7 for acknowledgements and details of the data sources and availability). Descriptive statistics for the data are shown in Table 3. Table 3: Descriptive statistics of the data Monthly 366 298 330 91 333 74
( 200)

Total no. of series mean length median length min length max length No. of series
(No. of observations)

Quarterly 427 99 110 30 130 125

( 100)

Yearly 518 24 27 11 47 112

( 20)

18
(201300)

302
(> 100)

375
(2130)

264
(>300)

31
(>30)

A subset of these series was used for evaluating the forecasting performance of methods that use explanatory variables. There were 93 quarterly series and 128 yearly series for which we had explanatory variables available. For each series we split the data into a test sample and a hold-out sample which was hidden from all of the co-authors. For each monthly series, the hold-out sample consisted of the 24 most recent observations; for quarterly data, it was the last 8 observations; and for yearly data it consisted of the nal 4 observations. Each method was implemented (or trained) on the test sample, and forecasts were produced for the whole of the hold-out sample for each series. The forecasts were then compared to the actual withheld observations. For each forecast horizon h, we rst consider the percentage better (PB) measure (as in Makridakis and Hibon, 2000). The PB shows the percentage of times that each method produces more accurate forecasts than SNave for monthly and quarterly data, and than Nave for yearly data.

16

The tourism forecasting competition

We also consider three alternative forecast error measures: the mean absolute percentage error measure (MAPE) MAPEh = 1 S
S s=1 s yh h ys s yh

and two scaled error measures suggested by Hyndman and Koehler (2006): the mean absolute scaled error (MASE) MASEh = and the median absolute scaled error (MdASE) MdASEh = median(ASEs ), h where S is the number of series and ASEs = h
s yh h ys 1 n1 n s i=2 | yi s yi1 |

1 S

ASEs , h
s=1

5 Results
5.1 Time series forecasting results

The tables that follow present the PB, MAPE, MASE and MdASE results.

The columns labeled

Average 1 h show the average forecast error measure over the forecast horizons 1 to h. The last column of each table, labeled Average rank, shows the average ranking for each forecasting method over all forecast horizons in the hold-out sample. To assist in evaluating the forecasting performance of the methods, we consider three aspects of the results: (i) the average rankings of the methods over all forecast horizons, i.e., the last column in each table; (ii) the rankings of the methods for the average error measures over the subsets 1 h considered in each case; and (iii) the performance of the methods for h = 1 step ahead. These aspects are considered in no particular order of importance, and lead us to some general conclusions.

Monthly data The results for monthly data are presented in Table 4, and are summarized as follows: Forecast Pro, ARIMA and ETS consistently forecast more accurately than SNave.

17

The tourism forecasting competition

When considering the MAPE, Forecast Pro produces the most accurate forecasts, but when considering the MASE, the ARIMA methodology is more accurate. The Theta method and the Damped trend method seem to be inferior to the other methods for forecasting monthly data. However, the Damped trend method is much more accurate than the Theta method and SNave for one-step-ahead forecasting.

Quarterly data A summary of the results for quarterly data presented in Table 5: Forecast Pro and ARIMA consistently forecast better than SNave. The Damped trend method forecasts quarterly data extremely well. It is consistently in the top two methods, regardless of the forecast error measure used. When it is applied to quarterly data, the Theta method is still generally outperformed by SNave. SNave seems to produce forecasts which are more accurate than those of any of the other methods for the seasonal horizons (i.e., h = 4 and h = 8). Even when considering the MdASE, SNave forecasts more accurately than ETS and ARIMA (though only marginally).

Yearly data A summary of the results for yearly data which are presented in Table 6: The Theta method is the only method that is competitive to Nave. When considering the MASE, it forecasts more accurately than Nave.

18

The tourism forecasting competition

Table 4: Forecast accuracy measures for monthly data

Method Forecast horizon (h) 1 2 3 PB relative to SNave ARIMA ForePro ETS Theta Damped 61.48 57.92 53.55 49.18 58.20 MAPE ForePro ETS ARIMA Theta SNave Damped 16.75 17.86 17.38 19.29 19.89 17.90 MASE ARIMA ForePro ETS SNave Theta Damped 1.00 1.02 1.19 1.23 1.35 1.08 MdASE ARIMA ForePro ETS Theta Damped SNave 0.78 0.82 0.89 0.95 0.78 1.01 0.89 0.86 0.87 0.94 0.93 1.18 1.02 0.83 0.91 1.01 0.98 1.05 1.01 0.95 0.96 0.93 1.19 1.05 0.77 0.85 0.82 0.85 0.81 0.85 1.13 1.14 1.23 1.12 1.51 1.21 1.09 1.12 1.19 1.17 1.08 1.13 0.90 0.84 0.89 0.97 0.90 1.08 0.90 0.89 0.91 0.96 0.98 1.02 1.02 1.01 1.03 1.06 1.11 1.14 2.71 2.88 3.13 3.50 4.00 4.79 1.16 1.17 1.26 1.43 1.50 1.36 1.28 1.25 1.32 1.40 1.71 1.67 1.29 1.30 1.40 1.47 1.43 1.71 1.07 1.12 1.14 1.09 1.15 1.08 1.68 1.69 1.88 1.78 1.73 2.19 1.45 1.54 1.61 1.48 1.48 1.47 1.15 1.14 1.26 1.35 1.52 1.37 1.21 1.22 1.30 1.37 1.42 1.47 1.38 1.40 1.49 1.54 1.55 1.66 1.67 2.25 3.75 4.33 4.33 4.67 16.22 17.30 17.65 20.11 21.56 19.03 17.17 18.30 18.45 20.30 20.64 22.25 17.32 20.89 19.13 20.20 20.94 26.53 20.54 20.44 21.09 21.02 21.09 20.70 17.11 19.74 18.02 18.50 19.97 24.19 23.27 23.65 24.29 22.51 22.30 22.35 16.71 17.82 17.83 19.90 20.70 19.73 18.38 19.67 19.37 21.02 21.38 22.30 19.91 21.15 21.13 22.11 22.56 23.47 1.46 2.92 3.17 4.21 4.54 4.71 63.11 62.30 59.02 53.83 61.20 58.20 59.02 58.74 49.45 55.46 55.46 53.28 56.28 52.19 48.09 53.01 52.73 53.28 51.64 50.27 53.01 52.19 50.27 54.92 43.99 58.20 55.46 51.09 55.46 54.10 60.93 59.74 57.10 50.82 58.29 57.10 56.17 54.94 51.68 52.80 56.96 56.01 54.63 53.19 53.20 2.00 2.33 3.00 3.67 3.79 6 12 18 24 Average 13 112 124 Average rank

19

The tourism forecasting competition

Table 5: Forecast accuracy measures for quarterly data

Method Forecast horizon (h) 1 2 3 PB relative to SNave Damped ForePro ARIMA Theta ETS 62.06 62.76 62.30 53.40 60.89 MAPE ForePro Damped ETS SNave Theta ARIMA 11.78 11.91 11.73 13.95 13.89 12.81 MASE ARIMA Damped ForePro SNave Theta ETS 1.10 1.11 1.13 1.34 1.47 1.19 MdASE Damped Theta ARIMA ForePro ETS SNave 0.92 1.04 0.87 0.83 0.90 1.15 0.93 0.99 1.01 1.00 1.01 1.08 0.90 0.88 0.86 0.85 0.81 0.90 0.89 0.92 0.93 0.91 0.94 0.92 1.12 1.26 1.31 1.39 1.34 1.57 1.45 1.34 1.45 1.52 1.50 1.39 0.91 0.96 0.92 0.90 0.92 1.01 1.08 1.11 1.11 1.14 1.14 1.21 2.75 2.88 3.00 3.63 4.25 4.50 1.30 1.18 1.29 1.45 1.43 1.36 1.18 1.21 1.16 1.22 1.22 1.17 1.24 1.23 1.22 1.18 1.28 1.30 1.80 1.60 1.79 2.08 1.74 1.96 1.80 1.81 1.86 1.79 1.79 1.99 1.21 1.18 1.20 1.30 1.35 1.26 1.47 1.43 1.48 1.59 1.56 1.58 2.63 2.63 2.88 4.00 4.13 4.75 12.38 11.68 12.59 14.79 13.90 12.72 13.99 14.85 13.70 14.41 14.47 14.67 14.21 14.21 14.78 13.61 15.07 14.79 15.05 13.83 15.88 18.02 15.24 16.13 22.90 22.28 24.01 21.15 21.71 22.21 13.09 13.16 13.20 14.19 14.33 13.75 15.72 15.56 16.05 16.46 16.15 16.23 2.63 3.13 3.25 3.75 3.88 4.38 60.66 55.74 58.08 51.52 52.46 54.80 52.69 52.22 49.65 53.40 52.69 52.22 49.18 48.48 47.78 67.68 60.42 57.38 59.72 57.38 54.57 52.69 52.46 54.33 50.82 57.55 55.85 55.44 50.76 53.63 58.02 55.91 56.15 53.45 53.81 1.50 2.75 2.75 3.88 4.00 4 6 8 Average 14 18 Average rank

20

The tourism forecasting competition

Table 6: Forecast accuracy measures for yearly data Method Forecast horizon (h) 1 2 3 PB relative to Nave ForePro ARIMA Theta ETS Damped 63.71 53.86 50.58 47.30 44.40 MAPE Nave Theta ForePro ETS ARIMA Damped 21.47 23.06 23.71 23.57 25.06 24.71 MASE Theta Nave ForePro ARIMA ETS Damped 1.32 1.32 1.49 1.56 1.50 1.55 MdASE ForePro Theta ETS Damped ARIMA Nave 1.06 1.10 1.09 1.13 1.14 1.10 1.49 1.56 1.62 1.61 1.70 1.62 2.28 2.21 2.29 2.34 2.35 2.43 2.88 2.69 3.01 2.91 2.92 3.16 1.28 1.33 1.35 1.37 1.42 1.36 1.93 1.89 2.00 2.00 2.02 2.08 1.50 1.75 3.75 3.75 5.25 5.00 1.96 2.08 2.18 2.20 2.22 2.29 2.63 2.95 3.10 3.05 3.13 3.23 3.20 3.64 3.85 3.70 4.01 3.92 1.64 1.70 1.83 1.88 1.86 1.92 2.28 2.50 2.65 2.63 2.71 2.75 1.00 2.00 3.50 4.00 5.00 5.50 20.80 21.17 22.49 23.26 25.32 24.41 24.12 22.94 27.28 28.56 28.06 29.43 28.05 26.61 31.96 35.35 33.69 34.05 21.14 22.12 23.10 23.41 25.19 24.56 23.61 23.45 26.36 27.68 28.03 28.15 1.50 1.50 3.25 4.50 5.00 5.25 70.85 62.55 60.62 55.02 54.83 72.20 65.06 70.85 58.11 59.46 71.24 64.48 71.04 57.53 60.42 51.16 58.20 67.28 55.60 49.61 54.49 61.49 69.50 63.27 54.78 1.00 2.50 2.50 4.50 4.50 Average 12 Average rank

14

21

The tourism forecasting competition

5.2 Does temporal aggregation improve the forecast accuracy?

The analysis so far has shown that Forecast Pro, ETS and ARIMA produce the most accurate forecasts for seasonal data. On average these methods produce more accurate forecasts than SNave. However, for yearly data none of the three methods forecast more accurately than Nave in our evaluatoin. In this section we investigate whether temporally aggregating the forecasts generated from these methods for higher frequency data, can produce more accurate forecasts for yearly data. We use all 366 monthly series and forecast h = 1 and h = 2 years ahead. The results from this investigation are presented in Table 7. Table 7: Comparing forecast errors from forecasting yearly data directly and temporally aggregating the forecasts produced for monthly and quarterly data ETS h=1 MAPE Yearly Quarterly to Yearly Monthly to Yearly Yearly from Nave 11.79 10.32 10.29 10.70 MASE Yearly Quarterly to Yearly Monthly to Yearly Yearly from Nave 1.50 1.37 1.36 1.43 MdASE Yearly Quarterly to Yearly Monthly to Yearly Yearly from Nave 1.21 1.09 1.08 1.27 1.85 1.72 1.71 1.93 1.16 1.06 1.05 1.27 1.72 1.65 1.63 1.93 1.30 1.11 1.11 1.27 1.87 1.78 1.78 1.93 2.25 2.09 2.08 2.17 1.43 1.28 1.28 1.43 2.01 1.89 1.89 2.17 1.49 1.29 1.29 1.43 2.15 2.05 2.04 2.17 16.49 14.32 14.29 15.01 10.99 9.94 9.93 10.70 14.59 13.98 13.96 15.01 11.44 9.95 9.92 10.70 15.36 14.48 14.46 15.01 h=2 ARIMA h=1 h=2 ForePro h=1 h=2

The rows labeled Yearly in Table 7 show the forecast error measures for h = 1 and h = 2 when forecasting the yearly data directly using each of the three methods. The rows labeled Monthly to Yearly and Quarterly to Yearly show the forecast error measures when forecasting the yearly data; however, the forecasts now come from aggregating the forecasts generated for the monthly and quarterly series, respectively. These results show that, in all cases, the aggregated forecasts (whether they are produced from the monthly data or the quarterly data) are more accurate than the forecasts produced from the yearly data directly. The rows labeled Yearly from Nave show the Nave forecast errors from forecasting the yearly data directly. In each case, the aggregated forecast error measures are smaller than the Nave ones.

22

The tourism forecasting competition

5.3 Forecast interval coverage

Producing estimates of uncertainty is an important aspect of forecasting which is often ignored in academic empirical applications, and is even more neglected in business practice. In this section we evaluate the performance of forecasting methods in producing forecast intervals that provide coverages which are close to the nominal rates. Tables 8, 9 and 10 show the percentage of times that the nominal 95% and 80% forecast intervals contain the true observations for monthly, quarterly and yearly data respectively. As was discussed in Section 1.1, forecasting methods often tend to overestimate the coverage probabilities of the forecast intervals they generate. This is the case with the forecast intervals produced by the ARIMA methodology for all frequencies. However, a signicant new nding here is that Forecast Pro and ETS produce coverage probabilities that are very close to the nominal rates for monthly and quarterly data. In fact, these methods slightly underestimate the coverage probabilities for the nominal 80% forecast intervals. As Makridakis et al. (1987) found, as we move to lower frequency data there is an increase in the tendency of methods to overestimate coverage probabilities. This is the case with all methods here. Table 8: Forecast interval coverage for monthly data Nominal coverage Forecast horizon (h) 1 2 3 4 5 Forecast Pro 95 95 93 85 83 82 ETS 95 93 83 82 ARIMA 89 85 77 70 92 84 83 66 Average 112 124

12

15

18

24

95% 80% 95% 80% 95% 80%

91 84 93 86 87 72

90 80 92 83 78 62

94 82 93 84 85 69

94 84 94 87 87 73

94 84 95 86 89 74

93 84 95 84 83 67

92 84 94 85 83 69

92 85 95 84 84 67

93 83 94 85 86 71

93 83 94 85 85 70

Table 9: Forecast interval coverage for quarterly data Nominal coverage Forecast horizon (h) 1 2 3 4 5 Forecast Pro 94 95 94 84 85 87 ETS 94 97 86 87 ARIMA 82 82 66 66 96 87 86 74 Average 14 18

95% 80% 95% 80% 95% 80%

94 84 95 85 83 67

93 83 95 84 82 67

91 83 93 87 78 63

92 78 92 84 82 67

90 77 92 81 78 63

94 85 95 86 83 68

93 83 94 85 82 67

23

The tourism forecasting competition

Table 10: Forecast interval coverage for yearly data Nominal coverage Forecast horizon (h) 1 2 3 4 ForePro 85 82 71 68 ETS 85 69 80 68 Average 14

95% 80% 95% 80% 95% 80% 95% 80%

76 62 76 63 67 49 68 51

74 57 75 62 65 47 64 50

79 65 79 66 69 51 71 55

ARIMA 72 71 55 53 Theta 78 73 61 57

Koehler (2001) stated that it would be interesting to see whether it is possible to nd statistically based forecast intervals for the Theta method. These were subsequently derived by Hyndman and Billah (2003). We have applied that result and produced forecast intervals from the Theta method for the annual data (as there is no seasonal component in the model), which are presented in the last two rows of Table 10. As with the other methods, the forecast intervals produced by the Theta method are badly undersized. We should note that there are alternative methods for constructing prediction intervals for each of the forecasting methods. We leave such comparisons for a separate study.

5.4 Forecast evaluation results: cases with explanatory variables

Quarterly data Table 11 shows the PB results, and Tables 1214 show the MAPE, MASE and MdASE results respectively. Some general observations from the analysis of these results are as follows. The pure time series approaches are consistently the most accurate. The best ranked of these approaches are Damped, Forecast Pro, ARIMA and ETS. Of the frameworks that use explanatory variables, the TVP DR and VARR (AIC) perform best. , For all forecast horizons h = 1 to 8, and for all forecast error measures, one of the pure time series approaches is always the most accurate.

24

The tourism forecasting competition

For h = 4, not many methods can forecast more accurately than SNave. The most inaccurate forecasts are those generated by the frameworks that do not perform any differencing. This sends a warning to forecasters who use quarterly variables in levels.

Yearly data Table 15 shows the PB results, and Tables 1618 show the MAPE, MASE and MdASE results respectively. Some general observations from the analysis of these results are as follows. No method can forecast more accurately than Nave when considering the MAPE. The pure time series approaches are consistently more accurate than the methods with explanatory variables. For the models with explanatory variables, TVP generates the most accurate forecasts. The three methods that compete with Nave are Theta, TVP and Forecast Pro. As was the case with the quarterly data, the results send a clear message to practitioners who model these types of variables in levels, as they are consistently the most inaccurate. Table 11: PB values for the quarterly cases with explanatory variables Method Damped ForePro ARIMA ETS Theta VARR (AIC) TVP m SR DR SR VAR(AIC) VAR(4) ADLM ADLMR SR Forecast horizon (h) 1 2 4 64.52 67.74 53.76 67.74 51.61 48.39 64.52 51.61 46.24 69.89 53.76 39.78 55.91 52.69 40.86 65.59 45.16 41.94 53.76 46.24 44.09 49.46 50.54 54.84 62.37 46.24 43.01 54.84 40.86 50.54 56.99 45.16 35.48 56.99 41.94 30.11 56.99 38.71 25.81 47.31 40.86 34.41 40.86 33.33 25.81 Average 14 60.22 55.91 55.65 53.76 47.85 51.88 46.51 52.15 47.04 46.77 44.35 41.40 39.52 39.52 30.38 Average rank 1.88 3.13 4.00 4.88 5.75 6.00 6.63 6.75 7.63 9.25 10.38 11.38 12.25 12.25 14.63

6 73.12 52.69 50.54 51.61 56.99 51.61 55.91 48.39 45.16 49.46 35.48 34.41 37.63 39.78 32.26

8 58.06 47.31 52.69 48.39 48.39 43.01 47.31 51.61 49.46 45.16 40.86 37.63 32.26 35.48 35.48

18 61.29 55.91 54.84 53.36 50.81 51.75 50.00 51.08 47.72 46.51 43.15 40.73 39.11 40.19 32.80

25

The tourism forecasting competition

Table 12: MAPE values for the quarterly cases with explanatory variables Method Damped ETS ForePro Theta SNave TVP DR VARR (AIC) ARIMA SR m SR VAR(AIC) ADLMR VAR(4) ADLM SR Forecast horizon (h) 1 2 4 11.10 8.33 9.83 9.58 10.00 10.71 11.42 10.62 9.85 12.03 10.75 10.78 13.83 11.78 8.96 10.35 12.12 9.89 10.23 10.55 12.24 10.94 11.14 11.25 14.25 11.26 11.21 10.94 12.46 10.63 14.90 12.95 9.19 10.79 12.25 17.11 12.54 14.29 15.83 11.43 13.21 17.62 11.92 14.06 17.87 25.18 27.03 29.67 Average 14 9.67 10.13 10.32 11.28 11.10 11.38 11.46 10.91 11.77 12.20 12.04 14.21 14.56 14.59 15.15 28.13 Average rank 3.63 4.38 5.13 5.25 5.38 5.63 6.50 7.63 8.13 8.38 10.00 11.25 11.75 13.25 13.75 16.00

6 11.79 15.31 14.99 14.15 16.14 15.37 17.96 17.08 16.62 17.31 20.76 22.03 21.53 23.68 22.92 35.11

8 29.21 32.65 30.00 26.52 27.47 27.65 26.31 38.12 27.08 32.73 38.06 35.40 33.40 37.00 38.70 48.56

18 11.90 11.81 12.72 12.14 12.62 12.58 13.25 12.47 14.13 13.62 14.08 17.41 16.65 17.90 18.09 30.56

Table 13: MASE values for the quarterly cases with explanatory variables Method Damped ForePro ARIMA ETS Theta SNave TVP VARR (AIC) DR m SR SR VAR(AIC) VAR(4) ADLMR ADLM SR Forecast horizon (h) 1 2 4 0.85 0.99 1.01 0.82 1.27 0.94 0.86 1.35 1.01 0.86 1.22 1.08 1.27 1.33 1.03 1.22 1.39 0.92 1.18 1.48 1.00 0.94 1.34 1.11 0.99 1.31 1.19 1.33 1.60 0.97 1.23 1.48 1.09 1.10 1.52 1.61 1.04 1.64 1.64 1.31 1.75 1.74 1.11 1.74 1.74 3.30 3.39 3.27 Average 14 18 0.99 1.28 1.03 1.34 1.09 1.40 1.08 1.45 1.28 1.48 1.18 1.51 1.30 1.53 1.13 1.56 1.24 1.69 1.29 1.76 1.36 1.75 1.52 2.12 1.52 2.19 1.65 2.20 1.61 2.23 3.29 3.94 Average rank 2.25 2.38 4.50 5.25 5.50 6.38 6.75 7.25 8.38 9.50 10.00 11.88 12.75 13.25 14.00 16.00

6 1.48 1.84 1.97 1.92 1.76 2.06 1.90 2.08 2.31 2.59 2.13 2.70 2.93 2.83 2.97 4.66

8 1.65 1.68 1.66 1.94 1.62 1.69 1.75 2.21 2.13 2.02 2.26 2.89 2.98 2.97 3.05 4.64

26

The tourism forecasting competition

Table 14: MdASE values for the quarterly cases with explanatory variables Method Damped ARIMA ForePro ETS TVP VARR (AIC) Theta DR SNave m SR SR VAR(AIC) ADLMR VAR(4) ADLM SR Forecast horizon (h) 1 2 4 0.61 0.74 0.69 0.61 1.02 0.66 0.62 1.02 0.64 0.60 0.98 0.74 0.73 0.94 0.70 0.65 1.04 0.76 0.67 1.14 0.84 0.69 0.90 0.79 0.92 0.97 0.75 1.01 1.24 0.68 0.81 0.99 0.78 0.72 1.14 1.22 0.83 1.24 1.03 0.78 1.36 1.23 0.76 1.16 1.28 1.42 2.32 1.57 Average 14 18 0.71 0.88 0.77 1.04 0.76 1.03 0.83 1.08 0.83 1.07 0.78 1.07 0.95 1.09 0.86 1.17 0.87 1.16 0.95 1.25 0.97 1.25 1.06 1.45 1.14 1.41 1.13 1.57 1.19 1.60 1.89 2.09 Average rank 2.50 4.25 4.75 5.38 5.63 6.25 7.00 7.13 7.75 8.75 9.38 11.63 12.38 13.25 13.75 16.00

6 0.96 1.49 1.63 1.54 1.31 1.62 1.41 1.54 1.66 1.86 1.52 1.88 1.93 2.12 2.17 2.20

8 1.22 1.30 1.40 1.33 1.54 1.58 1.34 1.30 1.45 1.14 1.61 1.88 1.76 1.88 1.95 2.51

Table 15: PB values for the yearly cases with explanatory variables Method ForePro ARIMA Theta VARR (AIC) Damped ETS TVP DR ADLM SR ADLMR SR VAR(AIC) VAR(2) Forecast horizon (h) 1 2 3 75.97 77.52 81.40 58.91 59.69 65.89 45.74 55.81 69.77 51.94 51.94 60.47 42.64 45.74 53.49 42.64 46.51 52.71 43.41 42.64 53.49 44.19 43.41 45.74 42.64 42.64 46.51 49.61 41.86 42.64 41.09 40.31 50.39 27.91 34.11 45.74 37.21 41.09 37.98 32.56 37.21 34.88 Average 12 76.74 59.30 50.78 51.94 44.19 44.57 43.02 43.80 42.64 45.74 40.70 31.01 39.15 34.88 Average 14 78.88 61.82 60.27 56.40 49.81 49.61 48.45 45.35 44.57 44.38 45.35 39.53 40.31 35.85 1.00 2.50 3.00 3.75 6.00 6.50 6.75 8.25 9.00 9.75 10.00 11.50 12.00 13.50

4 80.62 62.79 69.77 61.24 57.36 56.59 54.26 48.06 46.51 43.41 49.61 50.39 44.96 38.76

27

The tourism forecasting competition

Table 16: MAPE values for the yearly cases with explanatory variables Method Nave ForePro Theta TVP ETS Damped SR DR VARR (AIC) ARIMA ADLMR ADLM VAR(AIC) VAR(2) SR Forecast horizon (h) 1 2 3 49.13 32.65 27.93 51.94 33.23 28.61 54.22 35.16 27.82 51.27 36.16 30.86 51.60 34.45 35.01 54.59 36.52 32.17 49.11 35.80 35.67 53.11 35.39 35.02 55.97 35.31 36.20 55.45 45.03 33.88 63.47 42.70 36.71 64.72 40.18 39.18 65.23 43.48 47.29 72.59 47.73 50.73 99.77 76.70 49.41 Average 12 40.89 42.58 44.69 43.71 43.03 45.55 42.45 44.25 45.64 50.24 53.09 52.45 54.35 60.16 88.23 Average 14 35.67 36.88 37.56 39.16 40.59 40.16 40.86 41.31 42.05 44.71 46.52 47.43 52.62 57.33 72.49 1.50 3.25 3.50 5.00 5.25 6.50 6.50 7.00 7.75 9.75 10.75 11.50 12.75 14.25 14.75

4 32.97 33.74 33.05 38.34 41.31 37.39 42.84 41.70 40.72 44.49 43.22 45.64 54.47 58.29 64.07

Table 17: MASE values for the yearly cases with explanatory variables Method Theta Nave TVP ForePro Damped ARIMA ETS VARR (AIC) ADLMR SR DR ADLM VAR(AIC) SR VAR(2) Forecast horizon (h) 1 2 3 0.63 0.80 0.90 0.60 0.82 1.01 0.59 0.84 1.01 0.64 0.83 1.01 0.68 0.88 1.06 0.69 0.92 1.06 0.70 0.94 1.14 0.73 0.98 1.13 0.79 0.98 1.12 0.66 1.01 1.22 0.74 0.99 1.19 0.87 1.01 1.25 0.92 1.15 1.44 1.10 1.23 1.39 1.01 1.29 1.57 Average 12 0.72 0.71 0.72 0.74 0.78 0.80 0.82 0.85 0.88 0.84 0.86 0.94 1.03 1.17 1.15 Average 14 0.88 0.94 0.95 0.95 0.99 1.00 1.06 1.07 1.09 1.13 1.12 1.20 1.33 1.38 1.47 1.50 2.25 3.00 3.50 5.50 5.75 7.75 8.25 9.00 9.75 10.00 11.75 13.50 13.75 14.75

4 1.17 1.32 1.35 1.32 1.35 1.34 1.46 1.46 1.46 1.62 1.56 1.68 1.83 1.79 2.01

28

The tourism forecasting competition

Table 18: MdASE values for the yearly cases with explanatory variables Method Theta ForePro TVP Nave ARIMA Damped ADLMR ADLM ETS VARR (AIC) SR DR VAR(AIC) SR VAR(2) Forecast horizon (h) 1 2 3 0.46 0.41 0.72 0.42 0.49 0.76 0.42 0.62 0.72 0.41 0.53 0.75 0.48 0.60 0.78 0.48 0.52 0.84 0.63 0.68 0.73 0.67 0.66 0.79 0.50 0.69 0.98 0.53 0.75 0.85 0.49 0.78 1.00 0.54 0.77 1.02 0.63 0.75 1.04 0.98 0.93 1.03 0.71 0.90 1.13 Average 12 0.44 0.46 0.52 0.47 0.54 0.50 0.65 0.67 0.60 0.64 0.63 0.66 0.69 0.95 0.80 Average 14 0.63 0.65 0.68 0.68 0.71 0.71 0.74 0.77 0.84 0.81 0.92 0.93 0.92 1.08 1.08 2.25 2.75 3.50 4.25 5.50 6.00 6.25 8.00 9.25 9.25 11.25 11.50 11.75 14.00 14.50

4 0.94 0.91 0.94 1.02 0.98 1.00 0.92 0.98 1.20 1.12 1.41 1.37 1.24 1.38 1.58

5.5 Ex-ante versus ex-post forecasting for models with exogenous variables

The results from the models that include exogenous variables indicate that these models cannot forecast as accurately as pure time series approaches. In this section we perform ex-post forecasting for these models, i.e., we use the observed out-of-sample values for the exogenous variables. This eliminates any uncertainty related to the forecasting of the exogenous variables. Comparing the ex-post forecasts with the ex-ante forecasts allows us to evaluate whether the forecasting performance of these models improves enough for it to be deemed reasonable to use models with exogenous variables for scenario-based forecasting. In Table 19 we present the percentage improvements in forecast accuracy generated by the use of expost rather than ex-ante forecasting for models that use explanatory variables. Here, we only present the results for the three approaches we consider to be the most accurate based on the ex-ante results. In general, we nd the same puzzling result as Allen and Fildes (2001), namely that these methods become even less accurate when actual out-of-sample values are used for the exogenous variables. There are only two exceptions: for quarterly data, when h = 1 for TVP and SR for both MAPE and MASE. However, the improvement in these models for h = 1 does not change their rankings against the pure time series approaches. We should note that the methods that we have not presented here do not show any improvement when performing ex-post forecasting.

29

The tourism forecasting competition

Table 19: Percentage improvements of ex-post over ex-ante forecasting

Method Forecast horizon (h) 1 2 3 MAPE Quarterly data TVP DR SR 5.24 4.24 3.76 2.89 19.51 17.75 9.07 10.42 22.54 15.54 13.94 2.63 10.79 19.38 7.39 5.08 5.23 3.60 3.18 8.77 4.36 19.76 5.12 19.58 9.35 12.77 5.50 6.44 11.44 5.30 4 Average 1 2 3 MASE 4 Average

Yearly data TVP DR SR 15.73 46.56 63.45 77.44 47.36 20.66 51.58 59.09 74.85 50.03 19.07 47.83 50.33 57.71 41.61 18.20 33.95 43.26 46.74 38.50 23.10 36.95 44.82 48.51 41.20 23.00 38.35 44.54 45.11 39.81

5.6 Further warnings about the MAPE

The MAPE is clearly the most commonly used forecast accuracy measure in the forecasting literature. Hyndman and Koehler (2006) provide some warnings regarding this measure and highlight the conditions under which it is unsuitable and should not be used. The MAPE is dened only when all of the actual values being forecast are non-zero. In the M3-competition (Makridakis and Hibon, 2000), excessively large or undened MAPEs were avoided by using positive data only. Also, the MAPE (and all other measures based on percentage errors) assumes a meaningful zero. Both of these conditions hold with the tourism data we use here, and there is no reason for us to think that the MAPE is unsuitable. However, the MAPEs for the pure time series methods applied to the yearly data (Table 6) are less than half the size of the MAPEs obtained for those cases where explanatory variables are used (Table 16). This is due to nine series that contain low values and are amongst the cases with explanatory variables. These cases produced very large percentage errors, which caused the distribution of the MAPEs to be highly positively skewed. For example, the one-step-ahead forecast errors from the ETS models have a MAPE for all yearly series of 23.6%, whereas it is 51.6% for the cases with explanatory variables. When these nine series are excluded, the MAPEs are 16.3% and 23.5% respectively. Furthermore, when recalculating the MAPEs for the yearly series, excluding these nine cases, we nd that the rankings of the methods change. In particular, the Theta method forecasts more accurately than Nave, which makes the rankings of the methods identical to the rankings given by the MASE. Consequently, even though the MAPE can formally be applied here, we caution against its use due to the numerical instability that results whenever some series contain small values.

30

The tourism forecasting competition

6 Conclusion
We have designed a forecasting competition that uses data from the eld of tourism only. The forecast methods we consider are three fully automated time series algorithms (Forecast Pro, ARIMA and exponential smoothing based algorithms), two method-specic approaches (the Theta method and the damped trend), and ve general frameworks that incorporate explanatory variables (static and dynamic regression, autoregressive distributed lag models, time varying parameter models and vector autoregressions). We conclude that pure time series approaches forecast tourism demand data more accurately than methods that use explanatory variables. The forecasting performance of the models that use explanatory variables deteriorates when ex-post forecasting is performed instead of ex-ante forecasting. This has immediate practical consequences, as models with explanatory variables are commonly used in both the tourism literature and the tourism industry (especially for scenario-based forecasting). At the very least, the forecasting performance of these models should be evaluated against that of pure time series alternatives. The most consistent performance of any of the methods that use explanatory variables came from the TVP model. Of the pure time series forecasting approaches, we nd that Forecast Pro, ARIMA and ETS forecast more accurately than the seasonal Nave approach for seasonal data (both monthly and quarterly). It is interesting that this is the rst time in the empirical forecasting literature that an ARIMA based algorithm has produced forecasts as accurate as, or more accurate than, those of its competitors. For both of these seasonal frequencies, Forecast Pro and ETS produce forecast coverage probabilities which are satisfactorily close to the nominal rates. For yearly data, the Theta method is the only method that is competitive to Nave. Aggregating monthly or quarterly forecasts from one of Forecast Pro, ARIMA or ETS to give yearly forecasts also produced more accurate forecasts than Nave. Finally, we nd that the mean absolute percentage error distribution becomes highly positively skewed, due to a handful of small scaled series. Hence, we verify and endorse statements from previous researchers that the mean absolute scaled error should replace the mean absolute percentage error as the standard measure of forecast accuracy.

7 Acknowledgments
We would like to extend our gratitude to everyone who responded to our Request for data letter and supplied us with data. All of the time series data will be available from http://www.forecasters.org/
ijf. In order to adhere to all condentiality agreements with all parties, the data are presented under

coded titles.

31

The tourism forecasting competition

We acknowledge valuable comments from Michael Clements (editor), two anonymous referees and Andrew Maurer. We acknowledge Konstantinos Nikolopoulos and Nikolaos Bougioukos for supplying the forecasts from the Theta method. We thank Andrew Maurer and Claude Calero at Tourism Research Australia and Mike Chan and Bruce Bassett at New Zealands Ministry of Tourism for providing data and explanations. Athanasopoulos and Hyndman acknowledge nancial support from Tourism Research Australia and the Sustainable Tourism Cooperative Research Centre. Song and Wu would like to thank the Research Grant Committee of Hong Kong for nancial support (Grant No. BQ-04H).

References
Allen, P G. and R. Fildes (2001) Econometric forecasting, in J. Armstrong (ed.) Principles of Forecasting: A . Handbook for Researchers and Practitioners, pp. 271327, Kluwer Academic Publishing, Boston. Armstrong, J. (2001) Should we redesign forecasting competitions?, International Journal of Forecasting, 17, 542545. Armstrong, J. S. (1985) Long range forecasting, John Wiley, New York, 2nd ed. Armstrong, J. S. (2006) Findings from evidence-based forecasting: methods for reducing forecast error, International Journal of Forecasting, 22, 583598. Assimakopoulos, V and K. Nikolopoulos (2000) The theta model: a decomposition approach to forecasting, . International Journal of Forecasting, 16, 521530. Athanasopoulos, G., R. A. Ahmed and R. J. Hyndman (2009) Hierarchical forecasts for Australian domestic tourism, International Journal of Forecasting, 25, 146166. Athanasopoulos, G. and F. Vahid (2008) VARMA versus VAR for macroeconomic forecasting, Journal of Business and Economic Statistics, 26, 237252. Bohara, A. K. and C. Sauer (1992) Competing macro-hypotheses in the United States: a Kalman ltering approach, Applied Economics, 24, 389399. Box, G., G. Jenkins and G. Reinsel (2008) Time series analysis, Wiley, New Jersey, 4th ed. Brockwell, P J. and R. A. Davis (1991) Time series: theory and methods, Springer-Verlag, New York, 2nd ed. . Chateld, C. (2001) Prediction intervals for time series forecasting, in J. Armstrong (ed.) Principles of Forecasting: A Handbook for Researchers and Practitioners, pp. 475494, Kluwer Academic Publishing, Boston.

32

The tourism forecasting competition

Diebold, F. and L. Kilian (2000) Unit-root tests are useful for selecting forecasting models, Journal of Business and Economic Statistics, 18, 265273. Fildes, R. (1985) Quantitative forecasting the state of the art: Econometric models, The Journal of the Operational Research Society, 36, 549580. Fildes, R. and P Goodwin (2007) Against your better judgement? How organisations can improve their . use of management judgment in forecasting, Interfaces, 37, 570576. Fildes, R. and J. K. Ord (2004) Forecasting competitions their role in improving forecasting practice and research, in M. P Clements and D. F. Hendry (eds.) A companion to economic forecasting, pp. 322353, . Blackwell, Oxford. Gardner, Jr, E. S. (1985) Exponential smoothing: the state of the art, Journal of Forecasting, 4, 128. Goodrich, R. (2000) The Forecast Pro methodology, International Journal of Forecasting, 16, 533535. Goodrich, R. (2001) Commercial software in the M3-competition, International Journal of Forecasting, 17, 560565. Granger, C. W. J. (2001) Comments on the M3 forecast evaluation and a comparison with a study by Stock and Watson, International Journal of Forecasting, 17, 565567. Greenslade, J. V and S. G. Hall (1996) Modelling economies subject to structural change: the case of . Germany, Economic Modelling, 13, 545559. Harvey, A. C. (1989) Forecasting, structural time series models and the Kalman lter, Cambridge University Press, Cambridge. Hendry, D. F. (1986) Empirical modelling in dynamic econometrics, Applied Mathematics and Computation, 20, 201236. Hyndman, R. J. and B. Billah (2003) Unmasking the Theta method, International Journal of Forecasting, 19, 287290. Hyndman, R. J. and Y. Khandakar (2008) Automatic time series forecasting: The forecast package for R, Journal of Statistical Software, 26, 122. Hyndman, R. J. and A. B. Koehler (2006) Another look at measures of forecast accuracy, International Journal of Forecasting, 22, 679688. Hyndman, R. J., A. B. Koehler, J. K. Ord and R. D. Snyder (2005) Prediction intervals for exponential smoothing using two new classes of state space models, Journal of Forecasting, 24, 1737.

33

The tourism forecasting competition

Hyndman, R. J., A. B. Koehler, J. K. Ord and R. D. Snyder (2008) Forecasting with exponential smoothing: the state space approach, Springer-Verlag, Berlin-Heidelberg. Hyndman, R. J., A. B. Koehler, R. D. Snyder and S. Grose (2002) A state space framework for automatic forecasting using exponential smoothing methods, International Journal of Forecasting, 18, 439454. Kim, C. J. (1993) Sources of monetary growth uncertainty and economic activity: the time-varyingparameter model with heteroskedastic disturbances, Review of Economics and Statistics, 75, 483492. Koehler, A. B. (2001) The asymmetry of the sAPE measure and other comments on the M3-Competition, International Journal of Forecasting, 17, 570574. Kwiatkowski, D., P C. Phillips, P Schmidt and Y. Shin (1992) Testing the null hypothesis of stationarity . . against the alternative of a unit root, Journal of Econometrics, 54, 159178. Li, G., H. Song and S. F. Witt (2005) Recent developments in econometric modelling and forecasting, Journal of Travel Research, 44, 8299. Makridakis, S. and M. Hibon (2000) The M3-competition: Results, conclusions and implications, International Journal of Forecasting, 16, 451476. Makridakis, S., M. Hibon, E. Lusk and M. Belhadjali (1987) Condence intervals: an empirical investigation of the series in the M-Competition, International Journal of Forecasting, 3, 489508. Makridakis, S. and R. Winkler (1989) Sampling distributions of post-sample forecasting errors, Applied Statistics, 38, 331342. Ord, J. K., A. B. Koehler and R. D. Snyder (1997) Estimation and prediction for a class of dynamic nonlinear statistical models, Journal of the American Statistical Association, 92, 16211629. Ord, K. (2001) Commentaries on the M3-Competition: an introduction, some comments and a scorecard, International Journal of Forecasting, 17, 537584. Osborn, D. R., S. Heravi and C. R. Birchenhall (1999) Seasonal unit roots and forecasts of two-digit European industrial production, International Journal of Forecasting, 15, 2747. Pegels, C. C. (1969) Exponential smoothing: some new variations, Management Science, 12, 311315. Sims, C. A. (1980) Macroeconomics and reality, Econometrica, 48, 148. Song, H. and G. Li (2008) Tourism demand modelling and forecasting A review of recent literature, Tourism Management, 29, 203220.

34

The tourism forecasting competition

Song, H., S. Witt and T. Jensen (2003a) Tourism forecasting: accuracy of alternative econometric models, International Journal of Forecasting, 19, 123141. Song, H., S. Witt and G. Li (2003b) Modelling and forecasting the demand for Thai tourism, Tourism Economics, 9, 363387. Song, H. and S. F. Witt (2000) Tourism demand modelling and forecasting: modern econometric approaches, Oxford, Pergamon. Song, H. and S. F. Witt (2003) Tourism forecasting: The general-to-specic approach, Journal of Travel Research, 42, 6574. Song, H. and S. F. Witt (2006) Forecasting international tourist ows to Macau, Tourism Management, 27, 214224. Song, H. and K. F. Wong (2003) Tourism demand modeling: A time-varying parameter approach, Journal of Travel Research, 42, 5764. Tashman, L. (2001) The M3-Competition and forecasting software, International Journal of Forecasting, 17, 578580. Tay, A. S. and K. F. Wallis (2002) Density forecasting: A survey, in M. P Clements and D. F. Hendry (eds.) . A companion to economic forecasting, Wiley-Blackwell. Taylor, J. W. (2003) Exponential smoothing with a damped multiplicative trend, International Journal of Forecasting, 19, 715725. West, K. D. (2006) Forecast evaluation, Handbook of Economic Forecasting, 1, 99134. Witt, S. F., H. Song and P Louvieris (2003) Statistical testing in forecasting model selection, Journal of . Travel Research, 42, 151158. Witt, S. F. and C. A. Witt (1995) Forecasting tourism demand: A review of empirical research, International Journal of Forecasting, 11, 447475. Wong, K., H. Song, S. Witt and D. Wu (2007) Tourism forecasting: to combine or not to combine?, Tourism Management, 28, 10681078. World Tourism Organization (2008) UNWTO World Tourism Barometer, 6, 144. Zellner, A. (1963) Estimators for seemingly unrelated regression equations: Some exact nite sample results, Journal of the American Statistical Association, 58, 977992.

35