Adaptive tourism forecasting using hybrid artificial intelligence model: a case study of Xi’an international tourist arrivals

Method framework

Annual tourism demand typically exhibits a general trend with fluctuations, which reflect observations of complex dynamic systems. Factors like tourism price, exchange rate, and security can influence annual tourism demand. A visual representation of this framework is shown in Fig. 1. The time sequence is divided into fixed-length segments using a rolling window. The tendency is then extracted from each segment by a first-order gray model. The residual fluctuations are decomposed in the time series. The LSTM mechanism is used to predict the residual fluctuation, which has been learned from historical data. The prediction is made one step ahead, and this process can be iterated for multi-step ahead forecasting using a rolling window. To enable adaptive prediction that captures real-time system characteristics, Fixed-length time series are selected by a rolling window algorithm, on the basis of which analysis and prediction tasks are implemented. Our hybrid framework involves several steps. Firstly, we extract the tendency $\hat{x}$ of the data within the rolling window using a first-order gray model. Further, the residual fluctuations $\hat{\varepsilon }$ are decomposed out in the time series. In addition to this, these fluctuations are fitted by an LSTM neural network. Having learned from the historical data, we can then make one-step-ahead predictions, which can be iterated for multi-step-ahead forecasting using the rolling window approach. The models and related algorithms involved in this framework are described in detail in the following subsections.

Gray model

Gray models are constructed for application to complex systems and short-term tourism forecasting through differential equation methods. In grey theory, $GM\left(n,m\right)$ model is defined, where n is the order of the difference equation and m is the number of variables. The GM(1, 1) model, which utilizes first-order differential equations for univariate forecasting, is the most commonly used. Assume that the original time series is (1)${\displaystyle x^{\left(0\right)}=\left[x^{\left(0\right)}\left(1\right),x^{\left(0\right)}\left(2\right),\dots ,x^{\left(0\right)}\left(j\right),\dots ,x^{\left(0\right)}\left(t\right)\right]}$

where GM(1, 1) procedure forecasts the value $x^{\left(0\right)}\left(t+1\right)$ using the gray model, where $x^{\left(0\right)}\left(j\right)$ represents the j-th observation. And t is the total number of training data. The steps are as follows:

Construction of a sequence of first-order cumulative generation operations (1-AGO): (2)${\displaystyle x^{\left(1\right)}=\left[x^{\left(1\right)}\left(1\right),x^{\left(1\right)}\left(2\right),\dots ,x^{\left(1\right)}\left(j\right),\dots ,x^{\left(1\right)}\left(t\right)\right]}$

where $x^{\left(1\right)}\left(k\right)={\sum }_{j=1}^k{x}^{\left(0\right)}\left(j\right),$ for k = 1, 2, …, t. With this step we can obtain a smoother time series for the purpose of predicting the trend of the data.

Calculate the first-order differential equation: (3)${\displaystyle \frac{d{x}^{\left(1\right)}}{dt}+a{x}^{\left(1\right)}=u}$

where a represents the developing coefficient. u is the gray input.

Use ordinary least squares (OLS) to estimate the parameters a and u, with the solution being: (4)${\displaystyle \left[\begin{array}{c}\hfill a\hfill \\ \hfill u\hfill \end{array}\right]={\left({\mathit{B}}^T\mathit{B}\right)}^{-1}{\mathit{Y}}_N}$

where (5)${\displaystyle \mathit{B}=\left[\begin{array}{cc}\hfill -0.5\left(x^{\left(1\right)}\left(1\right)+x^{\left(1\right)}\left(2\right)\right)\hfill & \hfill 1\hfill \\ \hfill -0.5\left(x^{\left(1\right)}\left(2\right)+x^{\left(1\right)}\left(3\right)\right)\hfill & \hfill 1\hfill \\ \hfill \dots \hfill & \hfill \dots \hfill \\ \hfill -0.5\left(x^{\left(1\right)}\left(\mathrm{t}-1\right)+x^{\left(1\right)}\left(\mathrm{t}\right)\right)\hfill & \hfill 1\hfill \end{array}\right]}$

and (6)${\displaystyle {\mathit{Y}}_N={\left[x^{\left(0\right)}\left(2\right),x^{\left(0\right)}\left(3\right),\dots ,x^{\left(0\right)}\left(\mathrm{t}\right)\right]}^T.}$

Forecast model. (7)${\displaystyle {\hat{x}}^{\left(1\right)}\left(t+1\right)=\left(x^{\left(0\right)}\left(1\right)-\frac{u}{a}\right)e^{-at}+\frac{u}{a}}$

and (8)${\displaystyle {\hat{x}}^{\left(0\right)}\left(k\right)={\hat{x}}^{\left(1\right)}\left(k+1\right)-{\hat{x}}^{\left(1\right)}\left(k\right)}$

where ${\hat{x}}^{\left(0\right)}\left(k\right)$ indicates the predicted value at time k.

Rolling window scheme for gray model

Despite its ability to perform well with small samples, GM(1, 1) is limited in capturing real-time trends as it relies on the entire time-series for training. To overcome this issue and improve forecasting accuracy, a rolling window scheme (Akay & Atak, 2007) can be implemented. The approach involves training the model using only the most recent l observations, where l is the window size. This means that instead of using the entire time-series, $\left[x^{\left(0\right)}\left(t-l\right),x^{\left(0\right)}\left(t-l+1\right),\dots ,x^{\left(0\right)}\left(t\right)\right]$ is the basic work of the prediction ${\hat{x}}^{\left(0\right)}\left(t+1\right)$. The optimal window size 1 is computed by a grid search algorithm. The algorithm measures and evaluates the prediction performance using different window sizes.

Brief review on LSTM neural network

In 1997, Hochreiter and Schmidhuber proposed LSTM as a solution to the convergence issues that traditional neural networks face in time sequence prediction (Hochreiter & Schmidhuber, 1997). LSTM is a type of recurrent neural network (RNN) and its general structure is illustrated in Fig. 2A. RNNs excel in processing sequence data by incorporating the previous ‘memory’ of each state as input. The LSTM model comprises three parts: the input layer, the hidden layer, and the output layer. The inputs consist of x data from the previous $T=\left[x\left(t-l\right),x\left(t-l+1\right),\dots ,x\left(t\right)\right]$ historical moments. The output represents the target at the (t + 1)-th timestep, denoted as y (t + 1). H = (h₁, h₂, …, h_j, …, h_n) represents the neurons in hidden layer, where h_j refers to the j-th neuron in the layer. The weight parameter carried by the hidden layer is represented by W_hh.W_ih denotes the weight parameter between the input layer and the hidden layer. By the same token, the weight between the hidden layer and the input layer is represented by W_ho. The model is iteratively computed as follows:

(9)${\displaystyle h_n=\sigma \left(W_{ih}T+W_{hh}{h}_{n-1}+b_h\right)}$ (10)${\displaystyle y\left(t+1\right)=W_{ho}H+b_y}$

where b_hrefers to the size of the bias vector in the hidden layer. And the bias vector of the output layer is indicated as b_y.

Figure 2A illustrates that each neuron is comprised of three gates: the input gate, forget gate, and output gate, as depicted in Fig. 2B. The input gate selects the input information for the current cell state c_t from $x\left(t\right)$, c_t−1, and h_j(t − 1). The forget gate decides which information should be discarded based on the current input $x\left(t\right)$, the previous cell state c_t−1, and the previous output of the j-th neuron h_j(t − 1). The value of the current cell output gate is determined by the input size of c_t, h_j(t − 1) and $x\left(t\right)$. The input gate i_t and forget gate f_t are computed as follows:

(11)${\displaystyle i_t=\sigma \left(W_{ix}x\left(t-1\right)+W_{ic}{c}_{t-1}+W_{ih}{h}_j\left(t-1\right)+b_i\right)}$ (12)${\displaystyle f_t=\sigma \left(W_{fx}x\left(t\right)+W_{fc}{c}_{t-1}+W_{fh}{h}_j\left(t-1\right)+b_f\right).}$

The calculation of the neuron’s updated state is performed as follows: (13)${\displaystyle c_t=f_t\cdot c_{t-1}+I\cdot g\left(W_{ci}x\left(t\right)+W_{ch}{h}_j\left(t-1\right)+W_{cc}{c}_{t-1}+b_c\right).}$

The computation of the output gate is: (14)${\displaystyle o_t=\sigma \left(W_{ox}x\left(t\right)+W_{hh}{h}_j\left(t-1\right)+W_{oc}{c}_{t-1}+b_o\right)}$

where W_ab represents the value of the weight matrix from layer a to layer b. In addition to that, σ(.) denotes the sigmoid activation function, and each bias term with subscript is the bias vector. The weights of LSTM model are updated using back-propagation algorithms with training samples of the time series. LSTM networks utilize a gating mechanism to implement a selective memory function, which enables them to better handle time sequence prediction tasks compared to conventional neural networks.

GM-LTSM hybrid model

The primary objective is to capture the general trend using the gray model and represent the non-linear residual fluctuations with the LSTM model. For the training phase, we utilize the complete time series data $\left[x\left(1\right),x\left(2\right),\dots ,x\left(t\right)\right]$. In line with ‘Gray Model’, we start by extracting training samples S_k using a rolling window of length l, for one-step prediction. (15)${\displaystyle S_k=\left[x\left(k\right),x\left(k+1\right),\dots ,x\left(k+l\right)\right],k=1,2,\dots ,t-l-1.}$

where S_k is the k-th training sample.

To forecast the next data point $x\left(k+l+1\right)$ for each sample, a hybrid model combining GM(1, 1) and LSTM is employed. The time series data within the rolling window is first fitted using the GM(1, 1) model described in ‘Method Framework’, which produces a prediction $\hat{x}\left(k+l+1\right)$. However, as the gray model fails to capture fluctuations, a LSTM model is used to refine the prediction by incorporating the fluctuation. To calculate the residual error of the one-step-ahead prediction by the gray model for the k-th rolling window, we employ the following formula: (16)${\displaystyle \varepsilon \left(k+l+1\right)=x\left(k+l\right)-\hat{x}\left(k+l+1\right),k=1,2,\dots ,t-l-1.}$

During the prediction stage, we first collect the time series S_t−l $=\left[x\left(t-l\right),x\left(t-l+1\right),\dots ,x\left(t\right)\right]$ for the latest rolling window to predict the next data point $\hat{x}\left(t+1\right)$ using gray model. Additionally, we use LSTM to predict the corresponding residual fluctuation $\hat{\varepsilon }\left(t+1\right)$ based on the latest sample series. Thus, the predictions of these two machine learning models are properly combined and we obtain reasonable predictions. (17)${\displaystyle \tilde{x}\left(t+1\right)=\hat{x}\left(t+1\right)+\hat{\varepsilon }\left(t+1\right).}$

The same process used for prediction can be iteratively applied to make predictions for subsequent time series at multiple time steps ahead while reducing their weight.

ARIMA model

In order to evaluate the performance of our hybrid AI model, we have included this linear model for comparison purposes. When it comes to forecasting annual tourism demand, a non-seasonal ARIMA model with the form ARIMA(p, d, q) is typically employed. Among them, p represents the auto-regressive term, d is the integrated term and q is the moving-average term. The general formula can be expressed as follows: (18)${\displaystyle \left(1-{\varphi }_{1}B-{\varphi }_2{B}^2\dots -{\varphi }_p{B}^p\right){\left(1-\mathrm{B}\right)}^d{x}_t=\left(1-{\theta }_{1}B-{\theta }_2{B}^2-\dots -{\theta }_q{B}^q\right){\varepsilon }_t.}$

The data and random error terms at time t are represented as x_t and ɛ_t, respectively. The backward shift operator is denoted as B, and it is used to shift a time series backward by one period. Specifically, Bx_t = x_t−1. The difference levels are denoted as d. The autoregressive coefficients are represented as ϕ₁, ϕ₂, …, ϕ_p, and they capture the effect of previous values of the time series on its current value. The moving-average coefficients are represented as θ₁, θ₂, …, θ_q, and they capture the effect of past forecast errors on the current value of the time series.

Dataset description

According to the latest data, Xi’an is one of the most famous historical and cultural cities in China and one of the world’s most famous tourist destinations. Xi’an has a very large number of visitors every year, reaching over 360 million in 2019. For the dataset, we collected the number of international visitors to Xi’an from 1980 to 2018 as raw data. This information, which has been collected by the government of Xi’an and is presented in Table 1, will be used to assess the accuracy of our forecasting model. The data from 1980 to 2013 are used as the training set and the data from 2014 to 2018 are test set. Test data are {Year: Arrivals}: {2014: 124.23; 2015: 110.72; 2016:133.88; 2017:175.13; 2018:202.75}.

Evaluation

In this study, the data spanning from 1980 to 2013 was utilized as the training set, while the most recent five years from 2014 to 2018 were reserved as the test set to assess the efficacy of various models.

To evaluate the models’ performance, three commonly employed metrics were employed: mean square error (MSE), absolute mean error (MAE), and mean absolute percentage error (MAPE). These metrics are defined as follows:

(19)${\displaystyle MSE=\frac{1}{n}\sum _{i=1}^n{\left(x_i-{\hat{x}}_i\right)}^2}$ (20)${\displaystyle MAE=\frac{1}{n}\sum _{i=1}^n\left|x_i-{\hat{x}}_i\right|}$ (21)${\displaystyle MAPE=\frac{1}{n}\sum _{i=1}^n\left|\frac{x_i-{\hat{x}}_i}{x_i}\right|\times 100\text{\%}.}$

Model implementation

The models were coded and assessed in Python 3.7, which is open-source software. Here are the implementation details for each model: