Modelling and Forecasting Domestic Tourism. Case Study from Armenia

This paper summarizes the arguments and counterarguments within the scientific discussion on the issue of modelling and forecasting domestic tourism. During Covid-19 many countries tried to develop domestic tourism as an alternative to inbound tourism. In Armenia domestic tourism has grown recently, and in 2020 the decrease was 33% compared to last year. The main purpose of the research is to model and forecast domestic tourism growth in Armenia. Systematization of the literary sources and approaches for solving the problem indicates that many models and different variables are used to forecast tourism development. Methodological tools of the research methods were static and dynamic models, years of research were 2001-2020, quarterly data. The paper presents the results of an empirical analysis, which showed that with the static regression analysis a 1% change in GDP will lead to a change of 4.43% in the number of domestic tourists, a 1% change in the CPI will lead to a 14.55% change in the number of domestic tourists. For dynamic modelling we used 12 competing short-term forecasting models. Based on the recursive and rolling forecast simulation results we concluded that out-of-sample forecasts obtained by the small-scale models outperform forecasts obtained by the large-scale models at all forecast horizons. So, the forecasts of the domestic tourists’ growth obtained by small-scale models are more appropriate from the practical point of view. Then, in order to check whether the differences in forecasts obtained by the different models are statistically significant we applied Diebold-Mariano test. Based on the results of this test we concluded that there is not sufficient evidence to favor large-scale over small-scale models. This means that the forecast results obtained for domestic tourist growth by using the small scale models would not be statistically different from the results obtained by the large scale models. Based on the analysis, the forecasted values for domestic tourists for the future years were determined. The results of the research can be useful for state bodies, as well as private organizations, and for everybody who wants to model and forecast tourism development.


Introduction
In the recent years, tourism was growing at high rates in many countries. In the Republic of Armenia tourism is considered to be one of the main sectors of economy. During Covid-19 the sphere had a great downturn, as in the whole world. As an alternative to international tourism, many countries tried to develop domestic tourism, in order to save the sphere from collapse.
In Armenia the number of domestic tourists was 1544600 in 2019, and in 2020 it decreased to 1045756 (The socio-economic situation in the RA, 2021). The decrease was 33% compared to last year. Taking into account the severe restrictions during Covid-19, when tourism business was closed for several times, the results are not so bad. Also, because of the international restrictions, closed borders, limited flights, many Armenian tourists could not go abroad, so they had to spend their holidays inside the country.
In 2021, the international tourism is not fully restarted yet, some countries still continue to have closed borders, which mean, that domestic tourism will continue to be the main source of income for many countries.
In such situations, it is necessary to take some steps to promote domestic tourism, to make forecasts for its future development.
The aim of this research is to model and forecast domestic tourism growth in Armenia. For that purpose some static and dynamic models were used.

Literature Review
In the research some literature sources were analyzed about forecasting and modelling domestic tourism. Here is a brief summary of the main results of the analysis.
Some analyses show that the following factors can influence international tourism demand: income, relative prices between prices of origin and destination, exchange rates, relative prices between a destination and its competing destinations, cost of transportation, marketing expenditures, consumers' preferences, the effects of special events and other factors such as the effects of word of mouth. The main drivers of domestic tourism demand are discretionary income, consumer confidence indices and prices (Yap, 2010). Hamal (1996) argued that domestic holiday nights are strongly affected by tourists' income, prices of domestic goods and services, and prices of overseas holidays.
According to Dritsakis and Athanasiadis (2000) the main variable that definitely influences tourist movements positively is growth in Gross National Product (GNP).
According to Ghalehkhondabi et al. (2019) there is no forecasting method which can develop the best forecasts for all of the problems. Combined forecasting methods are providing better forecasts in comparison to the traditional forecasting methods.
A study by Song et al. (2008) shows that the time-varying parameter (TVP) model provides the most accurate short-term forecasts, whereas the naïve (nochange) model performs best in long-term forecasting up to two years.
The modern econometric models that have been used in tourism demand analysis are: Autoregressive distributed lag (ADL) model, Cointegration (CI) and error correction models (ECMs), Time varying parameter (TVP) model, Vector autoregressive (VAR) model, etc (Li et al., 2006;Li, 2009).
Univariate AR models are commonly used as benchmarks in the forecasting literature. Quite often AR models are perceived as advantageous compared to large multiple-equation models such as vector autoregression and traditional structural macroeconomic models (Hoffman, 2008;Arratibel et al., 2009). The univariate AR model can be estimated by using the following regression model: The unknown parameters of the model can be consistently estimated by using traditional OLS algorithm (Hamilton, 1994).
Another small-scale model that can be used for forecasting the quarterly growth of domestic tourists' is the unrestricted VAR model. A standard VAR with p lags is expressed as  (Hamilton, 1994). However, from the other side in the VAR model we very often need to estimate many parameters. This over-parametrization could cause inefficient estimates and hence a large out-of-sample forecast error. Thus, to overcome this over-parametrization we also implement the BVAR model. In order to use BVAR we first need to identify the priors. In this paper, we use the "Minnesota" type priors according to which the prior mean and standard deviation of the BVAR model can be set as follows: 1.
The parameters of the first lag of the dependent variables follow an AR(1) process while parameters for other lags are equal to zero.

2.
The variances of the priors can be specified as follows:  are set by the researcher and control the tightness of the priors. Thus, having "Minnesota" type priors it is possible to calculate the posterior parameters using the Bayesian approach to estimation (Poghosyan K., Poghosyan R., 2021).
Unlike small-scale benchmark models (AR, VAR and BVAR), the large-scale factor-augmented models include static or dynamic factors. As a rule, the factor-augmented models are estimated in two steps. First, we estimate the dynamics of unobservable factors using static and dynamic approaches and then we employ these extracted factors to forecast quarterly growth of the domestic tourists. In the modern time series econometrics literature there are three main algorithms for extracting factors, namely the static principal components as in Stock and Watson (2002), the dynamic principal component (frequency domain) approach as in Forni et al., (2000Forni et al., ( , 2005 and the dynamic principal component approach (time domain) as in Doz et al., (2011Doz et al., ( , 2012. In this working paper we using Stock and Watson (2002) and Doz et al., (2011Doz et al., ( , 2012 approaches to extract the factors dynamics. There are a number of papers that present the computational steps of these factor models in great detail (Doz et al., 2011(Doz et al., , 2012Barhoumi et al., 2014).
Second, we add the extracted factors into the small-scale benchmark models as additional explanatory variables. We can present the factor augmented model as follows: Where t X is the vector of observable variables, is the vector of unobservable variables estimated using the three previously mentioned methods, ,..., , 2 1 are ( r r  ) matrices of estimated parameters, t  and t u are the error terms with zero mean and diagonal variance-covariance matrices. The above presented model can be consistently estimated by OLS and Bayesian approach (Hamilton, 1994).

Static regression model of domestic tourism development
For static regression model the following variables were used: the number of domestic tourists in Armenia, the real GDP index, the consumer price index, the average cost of a tourist tour in other countries, which can be considered as an alternative to domestic tourism for the citizens of Armenia, and false variables were taken such as the global financial and economic crisis in 2008, the COVID-19 epidemic, and the war in Artsakh Republic (second Armenian country).
The indicators of the value of the tourist package from Armenia to other countries were obtained on the basis of surveys conducted among different travel agencies. Georgia, Egypt, UAE (Dubai), Greece, Spain and Italy were taken as destinations. Proxy approximation of their values was done.
The existence of a social package since 2012, which also includes a component of domestic tourism, has not been considered separately as a false variable, as it is already included in GDP.
The above database for the regression model was built for years 2005-2020, based on quarterly data. The time series used in the model were subjected to preliminary statistical processing, in particular, the natural logarithm of their absolute values was calculated, the series having seasonality were subjected to seasonal smoothing, and then the first differences of the smoothed series were calculated. The result is a stationary time series with quarterly percentage changes.
As the consumer price index also includes the prices of services, including the cost of hotels, so to avoid duplication, the cost of holidays in Armenian hotels was not included in the model. The results of the regression analysis show R 2 is 0.32, i.e. the change in the number of domestic tourists is explained by 32% by the model variables.
The cost of the tourist tour of other countries is not significant in the model, that is, it can be said that the latter does not affect the domestic tourism of the Republic of Armenia, and regardless of its price fluctuations, the t F domestic tourism will not change. This can be explained by the fact that a tourist who wants to spend his vacation abroad, regardless of price fluctuations, will still go on vacation, besides, in recent years the cost of that vacation has not changed significantly, it is not very high. On the other hand, as already mentioned, the above data were received from travel agencies and may contain averaged approximate data.
False variables 1 and 2, the global financial and economic crisis and the epidemic, are not significant, and false variable 3, the Artsakh war, is significant.
Thus, a regression equation whose unknown parameters have been estimated using the traditional least squares method will look like this: where DTRST is the number of domestic tourists, DGDP is real GDP of the country, DINFL is CPI, DPRICE is the price of tour packages in different countries, Dummy1 is the global financial and economic crisis, Dummy2 is COVID-19 pandemic, Dummy3 is war in Artsakh.
It turns out that a 1% change in GDP relative to the average value will lead to a change of 4.43% in the number of domestic tourists, a 1% change in the CPI to the average value will lead to a 14.55% change in the number of domestic tourists.
The results are logical and in line with the results of other studies, as the increase in income will lead to an increase in domestic tourism, and an increase in prices will lead to a decrease in domestic tourism.
In addition, dynamic models for domestic tourism forecasting have been used, which are described in more detail below.

Data descriptions of dynamic models used for domestic tourism forecasting
In this paper our main purpose is to forecast the growth of domestic tourists in Armenia. For that we consider three alternative approaches. First, we forecast domestic tourist's growth as univariate time series with using traditional AR model. Second, we forecast domestic tourist growth as a small-scale multivariate VAR and BVAR models including in the models additionally other key macroeconomic variables as an GDP growth, inflation and short-term nominal interest rates for loans. Third, we also forecast domestic tourist's growth with using large scale factor augmented models, including in the model besides of the main macroeconomic variables also extracted factors. Thus, our dataset includes four key macroeconomic variables, which we mainly use in the smallscale benchmark models, and 24 additional macroeconomic variables, which we use to extract the dynamics of unobservable factors. According to Barhoumi et al., (2014) our additional set of variables is medium-sized. Some studies have shown that the usage of smaller datasets which include about 10-30 series outperform the usage of larger datasets with disaggregated data with more than 100 series (Alvarez et al., 2016). Our dataset is balanced and in quarterly terms from 2001Q2 to 2020Q4 (79 observations in total for each variable).
Thus, as it was mentioned our aim is to forecast the growth of domestic tourists with using three alternative models and different dataset. Therefore, it is reasonable to consider the dynamics of domestic tourist's data. In order to calculate the growth of domestic tourist's data we have done the following calculations. First the absolute values of the domestic tourist's data were transformed to natural logarithms. Then transformed values were seasonally adjusted (X12 ARIMA method) and first differenced. As a result, we obtain the growth of domestic tourists ( Figure 1). As we can see, our dataset includes both 2008 world financial crisis and the period of COVID-19 pandemic. Including these two periods in our dataset increases the volatility of time series and therefore we can check the applicability of the above-mentioned models in the case of relatively large uncertainty environment.

Figure 1. Domestic tourist's growth (%-th change to the previous quarter)
Source: Compiled as a result of the above-mentioned analysis, the source of the data is the quarterly publications of the RA Statistical Committee.
The next important macroeconomic variable is the real GDP growth rate with respect to the previous quarter. In order to calculate the values of the growth rate we made the following calculations. First the absolute values of the real GDP in terms of average 2005 prices were transformed to natural logarithms. Then transformed values were seasonally adjusted (X12 ARIMA method) and first differenced. As a result, we obtain the real GDP growth rates ( Figure 2).

Figure 2. Real GDP growth (%-the change to the previous quarter)
Source: Compiled as a result of the above-mentioned analysis, the source of the data is the internal use databases of the Central Bank of Armenia The third macroeconomic variable is the CPI with respect to the previous quarter. It was first transformed to month-to-month indices to the base month indices, then these values were transformed to natural logarithms, the results were seasonally adjusted and the first differences were calculated. Figure 3 presents the seasonally adjusted CPI rates with respect to previous quarter. The fourth macroeconomic variable is the short-term nominal interest rate for loans in national currency. This time series is much more persistent than real GDP growth or CPI inflation. The preliminary treatment for this variable includes only first differences (in percentage points). As we see in Figure   Source: Compiled as a result of the above-mentioned analysis, the source of the data is the website of the Central Bank of Armenia: https://www.cba.am/am/SitePages/statmonetaryfinancial.aspx.
Besides of the above-mentioned variables, our dataset also includes additional 24 variables, as mentioned above. This set of additional variables can be grouped into the following categories: national accounts, consumer and producer prices, labor force and unemployment, financial variables. Most of the time series are obtained from CBA internal databases. A detailed data description is provided in Appendix 1. The vector of time series presented in Appendix 1 was preliminary treated. First, the time series were corrected for outliers and then seasonally adjusted as explained in Appendix 1. All nonstationary time series were transformed to be stationary, by taking first differences. We formally check the stationarity of the transformed additional series by using ADF (augmented Dickey-Fuller test) unit root tests. Based on the results of the tests we have concluded that all transformed series are stationary.

Experimental design
To conduct out-of-sample forecast evaluation we use both recursive and rolling regression schemes. For out-of-sample forecast evaluation we divide the whole sample into two subsamples, namely in-sample and out-of-sample periods. The first period is the training sample (in-sample), and the second period is the forecasting sample (out-of-sample). In our experiments the in-sample period includes 70% of observations, while out-of-sample period 30% of observations. The 70/30 proportion is a good compromise among the standard in-sample and out-of-sample proportions of 50/50, 70/30 and 90/10 broadly employed in modern machine learning algorithms (Brownlee, n/d). After choosing the proportion between in-sample and out-of-sample periods the recursive simulation scheme proceeds as follows.
First, we estimate the models using subsample 2001Q2-2015Q1 (56 observations). Using estimated model, we generate and then to store 1 to 4 steps-ahead forecasts results. Then we increase the sample size by one (57 observations, 2001Q2-2015Q2) and generate again 1 to 4 steps-ahead forecasts and then we store the forecast results. We continue increasing the sample size by one and generating 1 to 4 steps-ahead forecasts until the sample size 75 (2001Q2-2019Q4). Then we increase the sample size by one but only generate 1 to 3 steps-ahead forecasts (since we only have 76 observations in total). We continue increasing the sample size until we have 78 observations in the sample, in which case we can only compute the 1-step-ahead forecast. In this way, we obtain 23 1-step-ahead forecasts, 22 forecasts for 2-steps-ahead, 21 for 3-steps-ahead and finally 20 forecasts for 4-steps-ahead.
For the rolling forecast scheme, the initial sample is the same as in recursive scheme, but when the additional observation is added after the first forecast, the first values of the initial estimation sample are also deleted. Hence, while in the recursive scheme the sample size increases by one quarter at each step, in the rolling scheme the estimation size remains constant. The rolling regression scheme proceed as follows.
First we fix the sample size at 56 observations. As in the recursive regression scheme the forecasts are computed with a forecast horizon from 1 up to 4 and the results are stored. Then we add one observation to the sample and delete the first observations (in total we have 56 observations). Then we generate again 1 to 4 steps-ahead forecasts and the results are stored. Continuing in this manner we obtain the same number of forecasts as in the case of recursive regression. The recursive simulation scheme has an advantage of using all available data at a certain point of time, while the rolling scheme skips available information. But considering that our time series includes the episodes of global financial crisis and COVID-19 pandemic we have assumed that rolling regression scheme may be useful to use.
Next, we use the out-of-sample forecasts from recursive and rolling regression to compute the corresponding root mean squared forecast error (RMSFE) indices for each of the fourth forecasting horizons. More formally, we denote the outof-sample period by * T and forecast horizons h =1,2,3,4. Then the RMSFE index is calculated by the following formula: where, RMSFEis the root mean squared forecast error for the h-th forecast horizon, t ŷ is the forecasted value of the domestic tourists' growth, t y is the actual value of domestic tourists' growth (Kočenda, Poghosyan, 2020).

Empirical results
In this section we present the out-of-sample forecast results for 12 competing models. To keep robustness of our results we have estimated models with different lags length and different combinations of dynamic and static factors. Following (Pirshel and Wolters, 2014; Jos Jansen et al., 2016), we vary the number of lags from 1 up to 4 lags. In addition, we vary the number of static and dynamic factors, enabling us to consider all possible combinations. Thus, varying both the number of lags and the number of static and dynamic factors we compare the estimated models to each other. Finally, we select the lags length and number of static and dynamic factors by looking at the pseudo outof-sample forecast performances, particularly we select the combination that minimizes the RMSFE, evaluated over the entire out-of-sample period. The model with the smallest RMSFE index is selected as a model for forecasting at all horizons. Now, we will explain how we determine the number of static and dynamic factors in more detail for each model separately.
For the standard AR(p) model we only vary the number of lags from 1 up to 4. We conduct an out-of-sample forecast evaluation for four models, particularly models with 4 ,..., 1 = p lags and choose the model with the smallest RMSFE. We use the same approach for the VAR model. The only difference is that the VAR model includes four explanatory variables (domestic tourist's growth, real GDP growth, inflation, short-term nominal interest rate for loans). Again, we conduct an out of sample forecast evaluations for four models with 4 ,..., 1 = p . We chose the model with the smallest RMSFE. We use approximately the same approach for small-scale BVAR model. The only difference is that in addition to varying lags, we also vary overall tightness and lag decay. Following the paper by (Rangan, Kabundi 2011) we set the overall tightness 0.1 to 0.3 with increment equal to 0.1. The decay factor takes values of 1 and 2. Then we select the optimal combinations of lags and parameters by looking at the pseudo out-of-sample forecast performances and we select the model with the smallest value of RMSFE. Selected models with optimal lags and number of dynamic and static factors are presented in the Tables 1 and 2. For FAAR_SW model, in contrast above models, we vary both the number of lags and the number of static factors. Now we will explain how we estimate the number of static factors. Taking into account that additional set of variables is medium-sized, we estimate the number of static factors as follows. To estimate the number of static factors we retain in the analysis only the factors with eigenvalues more than 1 1  In other words, we can construct the following combinations for dynamic and static factors: 1 dynamic and 3 static, 2 dynamic and 3 static and finally 3 dynamic and 3 static factors. We use this idea to construct of all possible combinations of dynamic and static factors. Thus, if we estimate 8 static factors then according to our approach, we can have 36 combinations for dynamic and static factors in total. Taking into account that we also vary the number of lags ( 4 ,..., 1 = p ) then all possible combinations yield 144 FAAR_2S and 144 FAAR_QML models. We should note that the QML algorithm is an iterative procedure and for each model we run 100 iterations. Finally, we choose the model with the smallest RMSFE (Tables 2 and 3).
For the FAVAR_SW, FAVAR_2S and FAVAR_QML the selection procedures are the same as in the case of FAAR, with the only difference being that in the FAVAR models there are 4 target variables (domestic tourists' growth, real GDP growth rate, inflation rate, short-term interest rate). For BFAVAR_SW model we also vary additional hyperparameters (overall tightness and lag decay). As mentioned above, we vary overall tightness from 0.1 to 0.3 with increments equal to 0.1. The decay takes a value of 1 or 2. As in the case of previous models we choose the model with smallest RMSFE. For BFAVAR_2S and BFAVAR_QML the number of possible combinations is much higher, because for these models we also vary the number of dynamic factors. Again, we choose the model with smallest RMSFE and store the results in the Tables 2 and 3.   Source: The results were calculated using ForecastXL software package running in MS Excel, which was made by C # and VBA (Visual basic for application) programming languages. The package can be downloaded from the following website https://github.com/KarenPoghos/ForecastXL/.
The next question that arises from Tables 2 and 3 is whether the differences between forecasts generated by the large-and small-scale models are significantly different.
To give an answer on this question we have to conduct the equal forecast accuracy test, particularly Diebold-Mariano (1995) test. Before presenting the results of the tests we present some explanations related to this test.
In this paper we calculate the Diebold-Mariano statistic by regressing the loss differential on an intercept, using heteroscedasticity autocorrelation robust (HAC) standard errors (Diebold, 2015). Let AR t  denote the forecast errors in the benchmark AR(p) model and i t  denote the forecast errors in the competing i-th short term forecasting models (i = FAAR_SW, FAAR_QML, FAAR_2S). Then the loss differential t l can be calculated as Thus, we regress the loss differential on an intercept using HAC standard errors. The null hypothesis is that the loss differentials equal to zero    Note: *, ** indicate 10 and 5% level of significance.
Source: The results were calculated using ForecastXL software package running in MS Excel, which was made by C # and VBA (Visual basic for application) programming languages. The package can be downloaded from the following website https://github.com/KarenPoghos/ForecastXL/.
The statistics presented in Tables 4 and 5 indicates whether the performance results of large-scale and small-scale benchmark models are significantly different. From Table 4 and Table 5 we see that when we compare the predictive accuracy of the large-scale models with small-scale models then for both recursive and rolling regressions the differences are not statistically significant (with exception of only four cases for one step forecast horizon). In other words, there is not sufficient evidence to favor large-scale models over small-scale benchmark models. This means that the forecasting results for domestic tourist's growth rate obtained by the small-scale benchmark models could be just as good as the results obtained from models based on large data set.
Taking into account the analysis above, and the predictions made, the forecasts for the growth of domestic tourists for the quarters in 2021-2023 will be as follows. Source: The results were calculated using ForecastXL software package running in MS Excel, which was made by C # and VBA (Visual basic for application) programming languages. The package can be downloaded from the following website https://github.com/KarenPoghos/ForecastXL/.
According to calculations, in 2021 the total number of domestic tourists will be about 1,238 thousand people, in 2022 -1,369 thousand people, in 2023 -1,535 thousand people. Note that in 2020 the total number of domestic tourists was about 1 million people, and the projected figures seem realistic during Covid-19 pandemic.

Conclusion
We analyze the forecast performances of the 12 competing short-term forecasting models. In our analysis we generate ex-post out-of-sample forecasts based on the quarterly actual Armenian time series. For the ex-post outof-sample simulations we use both recursive and rolling regression schemes. Based on the recursive and rolling forecast simulation results we conclude that out-of-sample forecasts obtained by the small-scale models outperform forecasts obtained by the large-scale models at all forecast horizons. Based on these results we conclude that the forecasts of the domestic tourists' growth obtained by small-scale models are more appropriate from the practical point of view. Then, in order to check whether the differences in forecasts obtained by the different models are statistically significant we apply Diebold-Mariano test. We conduct this test both for recursive and rolling regression schemes. Based on the results of this test we conclude that there is not sufficient evidence to favor large-scale over small-scale models. This means that the forecast results obtained for domestic tourist growth by using the small scale models would not be statistically different from the results obtained by the large scale models.