Improved tourism demand forecasting with CIR# model: a case study of disrupted data patterns in Italy

Purpose – This study aims to predict overnight stays in Italy at tourist accommodation facilities through a nonlinear, single factor, stochastic model called CIR#. The contribution of this study is twofold: in terms of forecast accuracy and in terms of parsimony (both from the perspective of the data and the complexity of the modeling), especially when a regular pattern in the time series is disrupted. This study shows that the CIR# not only performs better than the considered baseline models but also has a much lower error than otheradditionalmodels orapproachesreported inthe literature. Design/methodology/approach – Typically, tourism demand tends to follow regular trends, such as low and high seasons on a quarterly/monthly level and weekends and holidays on a daily level. The data set consists of nights spent in Italy at tourist accommodation establishments as collected on a monthly basisby Eurostatbeforeandduring the COVID-19 pandemic breaking regularpatterns. Findings – Traditional tourism demand forecasting models may face challenges when massive amounts of search intensity indices are adopted as tourism demand indicators. In addition, given the importance of accurate forecasts, many studies have proposed novel hybrid models or used various combinations of methods. Thus, although thereareclear benefits in adopting more complex approaches, the risk is that of dealing with unwieldy models. To demonstrate how this approach can be fruitfully extended to tourism, the accuracy of the CIR# is tested by using standard metrics such as root mean squared errors, mean absoluteerrors, meanabsolute percentage errororaverage relative meansquarederror. Research limitations/implications – The CIR# model is notably simpler than other models found in literature and does not rely on black box techniques such as those used in neural network (NN) or data science-based models. The carried analysis suggests that the CIR# model outperforms other reference predictionsintermsofstatisticalsignificanceofthe error. Practical implications – The proposed model stands out for being a viable option to the Holt – Winters (HW)model, particularly whendealingwithirregular data.


Introduction
Tourism is one of the fastest growing industries in the world, the third largest export category and is crucial for many small developing countries which account for up to 50% of total exports (UNWTO, 2021a(UNWTO, , 2021b. Tourist arrivals in 2019 were 1.5 billion, confirming sustained growth for the tenth consecutive year. Tourism as a fair weather sector is heavily dependent on friendly framework conditions such as crime, hospitality and health (Burini, 2020). Frechtling (2012) classified factors affecting tourism into push, pull and resistance. According to Meleddu and Pulina (2016), pull factors pertain to the tourist destination, such as the quality of the natural resources and awareness of ecotourism. For Martins et al. (2017), push factors reflect both macroeconomic growth and specific conditions of the source market, such as leisure time, per capita income and relative prices. Conversely, resistance factors are those constraining travel from the source market to the destination such as corruption (Poprawe, 2015;Saha and Yap, 2015) and purchase power (Martins et al., 2017). The motivation to travel can be categorized into intrinsic and extrinsic motivations, each associated with different emotional experiences (Marino and Pariso, 2021).
Italy, with 64.5 million international tourist arrivals in 2019, was the fifth most visited destination in the world (see UNWTO, 2020). Tourism is a highly volatile industry depending on seasonality, social trends, connectivity and infrastructure. This volatility might be further exacerbated by both internal and external shocks. Regarding the latter, the industry was badly hit by the COVID-19 pandemic to the point that the World Tourism Organization (UNWTO) expected a decline of over 70% in 2020, back to levels of 30 years ago (UNWTO, 2020). Still in 2021, according to the latest UNWTO data, international tourist arrivals are expected to remain 70%-75% below 2019 levels (UNWTO, 2021a(UNWTO, , 2021b. Despite current difficulties, the tourism industry has experienced phenomenal growth over the last 30 years generating demand in modeling and forecasting. In fact, "accurate forecasts are critical for destinations where the decision-makers try to capitalize on developments in the tourism market and/or to balance their local ecological and social carrying capacities" Song et al. (2019). Song and Witt (2012) provided an account of econometric modeling methodologies, and Song et al. (2019) reviewed 211 key papers published between 1968 and 2018 with the intent to explain how the methods of forecasting tourism demand have evolved over time. The authors, therefore, concluded that there is no single method that works well for all situations and that the evolution of forecasting methods is still ongoing. That is why, for benchmarking, we decided to take six different baseline models, from time series to econometric and artificial intelligence models.
As mentioned by Law et al. (2019), "traditional tourism demand forecasting models may face challenges when massive amounts of search intensity indices are adopted as tourism demand indicators." In addition, given the importance of accurate forecasts, many studies have proposed novel hybrid models or used various combinations of methods. Thus, although there are clear benefits in adopting more complex approaches, the risk is that of dealing with unwieldy models. Among them, we mention data requirements and related issues (e.g. availability, cleansing and validation), estimation issues and model risk. In addition, typically, tourism demand tends to follow regular trends, such as low and high seasons on a quarterly/ monthly level and weekends and holidays on a daily level (Hu et al., 2021). This has been disrupted by the COVID-19 pandemic (Xie et al., 2020). In fact, unexpected excessive tourist demand and consumption place significant strain on resources and infrastructure. This strain can pose challenges for local authorities in terms of investment management and workforce recruitment, ultimately disrupting the stability of community economics. As a result, precise forecasting of tourist demand across different resources can greatly benefit researchers, industry workers and local authorities responsible for decision-making (Yao and Cao, 2020).
As a solution to the above-mentioned issues, in the present work, we propose a single factor, parsimonious stochastic model as a practical means of tourism forecasting. The suggested approach is robust to the breaking of regular patterns which is one of the main contributions to the current literature. Over the years, Orlando et al. (2018Orlando et al. ( , 2019aOrlando et al. ( , 2019bOrlando et al. ( , 2020 and Orlando and Bufalo (2021a) have developed the so-called CIR# model for forecasting time series in finance. The model, while parsimonious in terms of data and complexity, has proven robust when it comes to modeling cluster volatility, regime changes and spikes in time series that can typically be observed during turbulent periods such as the COVID-19 pandemic (Yonar et al., 2020). To demonstrate how this approach can be fruitfully extended to tourism, the accuracy of the CIR# is tested by using standard metrics such as root-mean-squared errors (RMSE), mean absolute errors (MAE), mean absolute percentage error (MAPE) or average relative mean squared error (AvgReIMSE). We show that the CIR# not only performs better than the considered baseline models but also has a much lower error than other additional models or approaches reported in the literature such as those by Gunter and Ö nder (2015), Kourentzes and Athanasopoulos (2019), Di Fonzo and Girolimetto (2022) and Wu et al. (2022). This holds true for the entire data set, as well as during the COVID-19 pandemic.
The main novelty of the proposed approach is its effectiveness as an alternative to the Holt-Winters (HW) model and artificial intelligence techniques, particularly in cases where data irregularity is a concern. Compared to many advanced existing models in the literature, the CIR# model is notably simpler and more transparent, avoiding the black box nature of neural network (NN) and data science-based models. In summary, we demonstrate that the CIR# model effectively handles the challenges posed by tourism time series when they exhibit characteristics such as positive kurtosis, nonnormality, autocorrelation and heteroscedasticity. This suitability for future purposes makes it valuable for tourism stakeholders who require reliable forecasts, particularly during periods of volatility in data patterns (Xie et al., 2020). The timing of the research is significant due to the increased need for accurate tourism demand forecasting, especially in the context of tourism disruptions caused by sudden travel restrictions, changes in consumer behavior and economic crises.
The remainder of the article is organized as follows. Section 2 provides a sketch of the existing literature. Section 3 presents the data set as well as its statistical characteristics. Section 4 is divided into three parts: proposed model, baseline models and measures of accuracy. Section 5 reports the obtained out-of-sample results. Section 6 provides a summary on the applicability of models and discusses results and the implications of the research. The last section concludes.

Literature review
As mentioned, Song et al. (2019) provided a comprehensive review of studies published on tourism demand forecasting. Essentially, forecasting methods fall into three categories: time series, regression and artificial intelligence models.

Time series models
Time series models such as the autoregressive integrated moving average (ARIMA) model have been used by Park et al. (2017) to internet search data from Google Trends to provide short-term forecasts for the inflow of Japanese tourists to South Korea. Generalized autoregressive conditional heteroscedastic (GARCH) model can be used in combination to forecast variances. Chan et al. (2005) applied ARIMA-EGARCH (exponential generalized autoregressive conditional heteroscedasticity) model from Japan, New Zealand, UK and USA to Australia confirming interdependency in the conditional variances between the four leading countries. Similarly, Coshall (2009) used an ARIMA-GARCH model for the UK demand for international tourism showing that volatility is relevant for tourism demand and that effects are asymmetric. The exponentially weighted moving average (EWMA) predicts future values based on appropriately weighted past observations, giving more importance to recent observations (Holt, 2004). Exponential smoothing methods have been developed by Holt (2004) and Winters (1960). The HW approach is an extension of EWMA incorporating linear trend and seasonality into the exponential smoothing (Rosy and Ponnusamy, 2017). Guojun and Ningning (2021) suggested exponential smoothing and ARIMA predict the number of tourists and the total amount of tourism consumption in Guangxi. Abdulmajeed et al. (2020), to forecast COVID-19 cases in Nigeria, devised an ARIMA model and a HW exponential smoothing model combined with a GARCH.

Regression models
Gunter and Ö nder (2015) compared the predictive accuracy of various uni-and multivariate models in forecasting international city tourism demand for Paris from its five most important foreign source markets (Germany, Italy, Japan, UK and USA). Specifically, they used seven different forecast models, i.e. error correction-autoregressive distributed lag (EC-ADLM), classical and Bayesian vector autoregression (BVAR), time-varying parameter (TVP), autoregressive moving average (ARMA) and exponential smoothing (ETS), as well as the naïve-1 model. However, the outcome is that there is no clear indication of the best model across countries and forecast horizons. Tratar and Strm cnik (2016), while studying heat load, found that multiple regression is the best for daily and weekly short-term forecasting, whereas HW method suits better for monthly and longer horizons forecasting.

Artificial intelligence models
For artificial intelligence approaches, a popular choice is machine learning based on NNs. Silva et al. (2019) claimed that hybrid singular spectrum analysis (SSA) combined with a neural network autoregression model (DNNAR) outperforms ARIMA forecasts in most cases. This is because "it is noise and not seasonality which hinders neural network autoregression (NNAR) forecasting capabilities." In this sense, SSA can be used for removing noise from data so that the NNAR model can be used on "smooth" data.
To tackle the problem of insufficient interpretability in tourism demand forecasting, Wu et al. (2022) proposed the use of a temporal fusion transformer (TFT) model, optimized using an adaptive differential evolution algorithm (ADE). TFT is a newly developed attention-based deep learning model that offers both high-performance prediction and time-dynamic interpretable analysis.

Alternative approaches
As the "no free lunch" theorem holds (see, e.g. Wolpert, 1996;Wolpert and Macready, 1997), one has to recognize that there is no single method outperforming the others on all scenarios in terms of accuracy and that all methods have their own limitations. For example, Law et al.
(2019) stated that "time-series and econometrics models rely on the stability of historical patterns and economic structure, while artificial intelligence models are dependent on the quality and size of available training data." Thus, the choice between different models depends on the type and quality of data having in mind the trade-off between complexity and improved modeling accuracy. Along these lines, Orlando et al. (2018Orlando et al. ( , 2019aOrlando et al. ( , 2019bOrlando et al. ( , 2020 and Orlando and Bufalo (2021b) have suggested an effective parsimonious nonlinear stochastic model called CIR# for modeling time series in presence of skewed distributions, jumps and cluster volatility. This, through an appropriate partitioning, shifting and calibration of the time series, has proved effective in terms of directionality of data and forecasting error even in case of disruption of data patterns.

Data set
The data set consists of nights spent in Italy at tourist accommodation establishments as collected on a monthly basis by Eurostat (2021). Figure 1 shows, for the Eurostat code NACE_R2, the total nights reported by hotels, holiday and other short-stay accommodations, camping grounds, recreational vehicle parks and trailer parks. The data exhibits a strong dependence on seasonality and displays an increasing trend. However, the outbreak of the COVID-19 pandemic caused a significant disruption in the pattern, as seen from record 260 onwards. Subsequently, there was a partial recovery in the data. This observation is further supported by Figure 2, which displays the seasonally adjusted data. Table 1 displays positive kurtosis and non-normality. The Kolmogorov-Smirnov test (normal distribution) and the Ljung-Box Q and Engle (ARCH) tests were conducted to evaluate the nights spent in Italy time series. The rejection decision, along with the corresponding pvalue, was obtained from these tests. Table 2 confirms that data are not normally distributed, exhibit autocorrelation and heteroscedasticity. When a time series exhibits bias, positive kurtosis, nonnormality, autocorrelation and heteroscedasticity, it can be challenging to model and predict using traditional statistical models like linear regression. In such situations, models specifically designed to handle these characteristics, such as the CIR# model, would be more appropriate for fitting the data and producing accurate predictions. The analysis of Tables 2 and 1 shows that the time series of nights spent in Italy demonstrates these characteristics, making the CIR# model a suitable choice for accurately modeling and predicting the data.

Methodology
As mentioned, we aim at predicting with a parsimonious model overnight stays in Italy at tourist accommodation facilities when a regular pattern in the time series is disrupted. To this end, we have selected a nonlinear, single factor, stochastic model called CIR#, and we compare it to others ranging from simple EWMA to HW; from Note: Total number of observations = 332 (monthly data) Source: Eurostat (2021) j TOURISM REVIEW j integrated ARIMA combined with EGARCH models, to seasonal integrated autoregressive moving average (SARIMA) and DNNAR models (see Appendix). As already mentioned, multiple regression is found to be effective for short-term forecasting on a daily and weekly basis, while the HW method is more suitable for monthly data and longer-term forecasting (Tratar and Strm cnik, 2016). Therefore, we have not included multiple regression from our set of baseline models. Finally, in Section 5, we also offer a comparison with other models used in the literature, including EC-ADLM, VAR, BVAR, TVP models and TFT model, optimized using an ADE as explored by Gunter and Ö nder (2015) and Wu et al. (2022), respectively.

The proposed model: the CIR#
The CIR process is a stochastic differential equation describing the evolution of the stochastic variable V(t) as introduced by Cox et al. (1985).
where V(0) ¼ V 0 > 0 is the initial condition and W(t) a Wiener process. The CIR process is said to be a single-factor, time-homogeneous model, because s, k and u are timeindependent, and V(t) is the only variable. The CIR model has been devised for fixed income pricing, but in the context of this work, the stochastic variable is tourist accommodations described in Section 3.
The proposed CIR# preserves the structure of the original CIR model by Cox et al. (1985) and involves partitioning the market data into subsamples to capture statistically significant changes in variance and jumps. The second step involves fitting an "optimal" ARIMA model to each subsample of market data and applying Johnson's transformation (Johnson, 1949) to the standardized residuals to ensure they resemble Gaussian white noise. This is necessary because empirical excess returns of financial assets often have more kurtosis and positive serial correlation (Orlando and Bufalo, 2021b). For more details, a reader can refer to Orlando et al. (2019aOrlando et al. ( , 2019b, Bufalo (2021b, 2023).

Baseline models
Below is the description of some basic models to compare their performance with those of the CIR# model. Notice that, throughout this section, 1,n] are the observations related to the explanatory variable V(t) over n periods.
4.2.1 Autoregressive integrated moving average exponential generalized autoregressive conditional heteroscedasticity model. An ARIMA(p, i, q) model is described by: where L is the lag operator, f j and u j indicate the parameters of the autoregressive part and of the moving average part of the model, respectively, and « are white noise terms.
The EGARCH(a, b) can be expressed as: where the variable« follows a generalized Gaussian distribution, s 2 is the conditional variance process and v, g, d, and j are real parameters. The ARIMA combined with EGARCH is suitable for time series data with heteroscedasticity, volatility clustering and conditional skewness or kurtosis. The ARIMA component models the autocorrelation and trend, while the EGARCH component models the volatility and asymmetry in the error terms. This combined model can provide more accurate forecasts and account for the nonlinear patterns and dynamics present in financial time series data.
4.2.2 Seasonal integrated autoregressive moving average model. The ARIMA model is an extension of the classical ARIMA, which is capable of modeling a wide range of seasonal data. SARIMA models are capable of capturing complex patterns in the data and can be used to forecast future values with a certain level of accuracy. However, they may not be suitable for nonstationary data or data with irregular patterns.

4.2.3
Exponentially weighted moving average model. The EWMA (see e.g. Perry, 2010) is a weighting scheme to simulate future values averaging on a historical data set. Holt (2004) suggested that the EWMA address trends and seasonality in forecasts. The EWMA is a means of smoothing out random fluctuations with some desirable properties such as: j TOURISM REVIEW j decrease in weight attributed to older data; ease of computation; and minimum data requirement.

Holt-Winters model.
The HW model is a generalization of the EWMA model in presence of seasonal data (Chatfield, 1978). There are two versions of such a method, the so-called additive version and the multiplicative one. The additive method is best suited when seasonal variations are roughly constant across the time series, as observed in our data set (see Subsection 3). The multiplicative method is preferable when seasonal variations change proportionally to the level of the time series. The method is useful when the data has a strong seasonal pattern that repeats over time. It can also handle trend and level changes in the time series. However, the method assumes that the time series is stationary, which means that its statistical properties do not change over time. It may not be suitable for nonstationary data with trends and seasonality changes.

Neural network autoregression model. A NN is a network of neurons structured in
layers where the predictions (outputs) form the upper layer and the predictors (inputs) form the lower layer. Often there are also some layers of neurons between predictors and predictions and which are called hidden layers. This configuration is called a multilayer feed-forward network.
As explained in Silva et al. (2019), the NNAR model has to be implemented jointly with the SSA to denoise data. SSA is a nonparametric estimation method where the covariance matrix is decomposed into a spectrum of eigenvalues. The resulting (denoised) time series is then used as the input for a NNAR model. This hybrid approach is called DNNAR, and in the remainder of this work, the "optimal" DNNAR model is denoted by DNNAR Ã . The details of the performed SSA are available in the Appendix. The DNNAR model is particularly useful for time series with high volatility and nonlinearity, where traditional statistical models may struggle to provide accurate forecasts.

Measures of accuracy
Here, we list the forecast accuracy measures used to compare the results between the considered models: The RMSE measures the closeness between the observed data and their predictions.
The MAPE is another measure of the prediction accuracy of a forecasting model.
The AvgReIMSE compares two different forecasts.
The AvgRelMSE is symmetric to over and underforecasting (Davydenko and Fildes, 2013) and is used in the literature to test different scenarios (Kourentzes and Athanasopoulos, 2019). As a reference, a given model A is more accurate than the selected benchmark B when the AvgReIMSE is smaller than 1. In particular, the considered model is better than the benchmark by (1-AvgRelMSE)100%. A variant of the AvgReIMSE is the AvgReIRMSE where the RMSE is taken instead of the MSE.

Out-of-sample results on the full data set
In this section, we show the out-of-sample performance of the proposed model against the results obtained with the five baseline models listed in Section 4.2. Concerning the CIR# model, to predict the future value, first, we calibrate the parameters (k, u, s, p, i, q) through a rolling window W of length M ¼ 12 so that W ¼ {v h , . . .., v hþMÀ1 (h ! 1). Then, the forecasted values v F hþMþs s ! 0 ð Þ, are determined with the procedure explained in Orlando et al. (2019b). Similarly, for the baseline models, we first calibrate the parameters on the same window and then we compute the next value (i.e. the forecast for the next month). In this regard, Figure 3 shows the out-of-sample forecasts versus actual data, while Figure 4 displays the error. As can be seen: the EWMA is not a suitable choice for this type of data; the SARIMA, HW and DNNAR models are quite good until the COVID-19 pandemic where they perform poorly; and CIR# and ARIMA-EGARCH seem to perform better even during the pandemic.
In particular, there are two jumps from record 260 onward that are captured by CIR# and ARIMA-EGARCH models.
To better appreciate the performance of the CIR# versus the baseline models, Table 3 reports the results of the measures listed in Section 4.3. In agreement with the above graphical analysis, the CIR# appears to be the best model in terms of all measures (MAE, RMSE and MAPE). Concerning the literature, we refer to Gunter and Ö nder (2015) who used multiple benchmarks such as EC-ADLM, VAR, Bayesian VAR, TVP models (multivariate or econometric models), ARMA and the ETS models (univariate or time-series models). In their tests, Gunter and Ö nder (2015) reported that the minimum prediction error (RMSE and MAE) of the considered models over a month horizon is about 7% and that the average is 10% (i.e. at least 3.5 times higher). Other models, like the TFT that uses an ADE, achieve a best MAPE of 3.02% (Wu et al., 2022), which is twice as high as that of CIR#.
Finally, Table 4 shows the high accuracy in terms of AvgReIMSE (as well as the AvgReIRMSE) of the CIR# model with respect to the selected baseline models. For reference, Kourentzes and Athanasopoulos (2019), when comparing the forecasts for Australian tourism against the   (2022), when performing a test on the accuracy of the proposed forecast combination-based approach, found an AvgReIMSE equal to 0.9618 (i.e. much higher than any value reported in Table 4).

Out-of-sample results amid COVID-19
In the previous section, we have shown the out-of-sample performance on the full data set of the proposed model against the results obtained with the five baseline models and   reported prediction errors in other studies of the literature. In this section, we focus on the part of the data set where the data pattern was disrupted because of COVID-19 pandemic. Figure 5 presents the out-of-sample predictions compared to the actual data, whereas Figure 6 illustrates the corresponding errors. As previously seen in Section 5.1 during the COVID-19 pandemic, it can be observed that: The EWMA model is not appropriate for this type of data.
The SARIMA, HW and DNNAR models perform well until the data pattern is significantly disrupted.
The CIR# and ARIMA-EGARCH models demonstrate superior performance in any condition, with the former outperforming the latter. Table 5 presents the results of the measures listed in Section 4.3 to provide a better understanding of the performance of the CIR# model compared to the baseline models. As confirmed by the graphical analysis discussed earlier, the CIR# model outperforms all other models in terms of all measures (MAE, RMSE and MAPE). It is worth noting that these results are superior to those reported in other studies such as Gunter and Ö nder (2015) and Wu et al. (2022).
Finally, Table 6 shows the high accuracy in terms of AvgReIMSE (as well as the AvgReIRMSE) of the CIR# model with respect to the selected baseline models which, once again, are better than those found in the mentioned literature (Kourentzes and Athanasopoulos, 2019;Di Fonzo and Girolimetto, 2022).

Model selection and distribution of forecast errors
The analysis presented above indicates that the CIR# model consistently has a small error compared to other reference predictions. However, it is important to determine whether this Note: Out-of-sample results over monthly data (changes) during the COVID-19 period difference is statistically significant or if it could be attributed to the specific sample of data used. Table 7 demonstrates that the differences between the CIR# forecasts and the reference predictions are statistically significant. Furthermore, Table 8 shows that, except for the SARIMA model, all the models produce homoscedastic errors indicating that the variance of errors is constant across all levels of the independent variable. This allows for more accurate statistical modeling and estimation of parameters, as well as valid inference and hypothesis testing. However, it is worth noting that the forecasting errors of all the models are not normally distributed.

Discussion and implications of the research
Analysis of data shows that the nights spent in Italy time series have a positive kurtosis and exhibit nonnormality. When a time series presents bias, positive kurtosis, nonnormality, autocorrelation and heteroscedasticity, it can be challenging to model and predict using conventional statistical models like linear regression. In such cases, models specifically developed to handle these characteristics, such as the CIR# model, would be more appropriate for fitting the data. Tables 1 and 2 confirm that the nights spent in Italy time series exhibit these characteristics, reinforcing the need for models that can handle such challenges. Therefore, it is not a case that results in Section 5 prove that the CIR# model performs better than the baseline models in terms of all measures. Furthermore, the carried analysis suggests that the CIR# model outperforms other reference predictions in terms of error the statistical significance of this difference. The proposed model has demonstrated superiority even when compared to other models in the literature and can be especially useful for tourism stakeholders in making decisions when there are disruptions in data patterns. Although all models, except for SARIMA, produce homoscedastic errors, it is worth noting that the forecasting errors of all models are not normally distributed. This assumption is important for statistical modeling, and violating it can lead to biased estimates, invalid inferences and inaccurate predictions. Such a result is not unexpected because, by design, the CIR# model is built in a way that standardized residuals resemble Gaussian white noise (Orlando et al., 2019a(Orlando et al., , 2019b.
In summary, the CIR# model offers a valuable addition to the body of knowledge regarding tourism demand forecasting, and its use could lead to more accurate predictions and better-informed decisions for stakeholders in the tourism industry. In fact, accurate forecasting of tourist demand is invaluable for researchers, industry workers and decisionmakers as it helps mitigate the challenges posed by excessive unexpected demand during peak periods and the underutilization of capacity during low-demand periods.

Conclusion
In this article, we have explained the relevance of tourism both in economic and social terms. Then, we have presented the case of Italy and we have walked through relevant literature for forecasting tourism demand. As mentioned by Song et al. (2019), the evolution of forecasting methods is still ongoing. All methods have their limitations, and there is no single method outperforming the others (Law et al., 2019). Moreover, because there is no free lunch (Wolpert, 1996;Wolpert and Macready, 1997), there is a trade-off between complexity and accuracy. Ever-increasing sophisticated models may suffer from model risk such as incorrect specification, wrong implementation, lack of sufficient data and calibration errors. The CIR#  model is parsimonious as it requires a single time series, and its forecasting power is tested against numerous benchmarks common in the literature. The result is that, for the case of Italy at the time of the COVID-19 pandemic, the proposed stochastic approach compares well against all the other considered baseline models providing a forecast error reduced by 70%. The same applies when comparing the obtained results to those available in the literature (e.g. see Gunter and Ö nder, 2015;Kourentzes and Athanasopoulos, 2019;Di Fonzo and Girolimetto, 2022;Wu et al., 2022).
The proposed model stands out for being a viable option to the HW model, particularly when dealing with irregular data. In addition, the CIR# model is notably simpler than other advanced models found in literature and does not rely on black box techniques such as those used in NN or data science-based models. In this sense, we have presented a simple but not simplistic model that may be added to the range of tools for predicting tourism volumes. Due to the weight of tourism in Italy, this has important implications in terms of development policies and regulations. Next research will explore the suitability of the CIR# model to forecasting tourism time series in other countries.
The limitations of the proposed approach arise when the underlying data of the tourism time series do not exhibit significant disruptions. In such cases, the CIR# model may not provide any significant advantages over conventional models such as the HW. Therefore, when interpreting the results, it is crucial to consider the specific characteristics of the data, the context of the research and the timing of the analysis. These factors can influence the suitability and performance of the hypothesized models, offering insights into the limitations and applicability of the findings.

Appendix. Singular spectrum analysis
In this appendix, we show the details of the singular spectrum analysis procedure, applied to our time series to denoise data as suggested by Silva et al. (2019): First of all, we have to normalize our time series, so we introduce the series z h ¼ v h Àm s , wherem;ŝ are the sample mean and sample standard deviation of the observations (v h ) h [[1,n] , respectively.
Then, we compute the covariance matrix C (see Figure A1). To do this, we determine C by the following scalar product: where the matrix Y is the time-delayed embedding of (z h ) h [[1,n] , i.e.: Y h;k ¼ z hþsÀ1 h 2 1; n À s þ 1 ½ ; k 2 1; s ½ ð Þ : Notice that, even though the estimated C matrix does not have a Toeplitz structure (with nonsymmetric or antisymmetric eigenvectors), it at least guarantees that C is positive semidefinite: Figure A2 shows the eigenvalues and eigenvectors of matrix C. Figure A3 shows the principal components as obtained by the scalar product between Y, the time-delayed embedding of (zh)h [[1,n], and the eigenvectors P. Figure A4 shows the reconstructed components, as obtained by inverting the projecting PC ¼ Y s P.
The upper panel of Figure A5 shows the original time series (vh) h [[1,n] as obtained by the sum of all reconstructed components. The lower panel of Figure A5 displays the denoised time series obtained with the first pair of reconstructed components.
j TOURISM REVIEW j