Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
53,909 characters · 15 sections · 91 citation commands
Neural Network Modeling for Forecasting Tourism Demand in Stopića Cave: A Serbian Cave Tourism Study
\end{tabular}\hfil[0]\hfil
\hfil[0]\hfil
\hfil[0]\hfil
\hfil[0]\hfil
\hfil[0]\hfil
\hfil[0]\hfil
\hfil[4]\hfil
\hfil[4]\hfil
\hfil[0]\hfil
\hfil[0]\hfil
\hfil[4]\hfil
\hfil[4]\hfil
\hfil[0]\hfil
\hfil[0]\hfil
\hfil[4]\hfil
\hfil[4]\hfil
\hfil[0]\hfil
\hfil[0]\hfil
\hfil[4]\hfil
\hfil[4]\hfil
\hfil[0]\hfil
\hfil[0]\hfil
\hfil[4]\hfil
\hfil[4]\hfil
\hfil[0]\hfil
\hfil[0]\hfil
\hfil[4]\hfil
\hfil[4]\hfil
\hfil[0]\hfil
\hfil[0]\hfil
\hfil[4]\hfil
\hfil[4]\hfil
\hfil[0]\hfil
\hfil[0]\hfil
\hfil[4]\hfil
\hfil[4]\hfil
\hfil[0]\hfil
\hfil[0]\hfil
\hfil[4]\hfil
\hfil[4]\hfil
\hfil[0]\hfil
\hfil[0]\hfil
\hfil[4]\hfil
\hfil[4]\hfil
Modeling tourist demand includes a complex but necessary set of activities and analyses that can potentially determine market norms and directly shape tourist offers (a1_buhalis2005tourism, a2_frechtling2012forecasting). Research on tourism demand states that understanding tourist demand enables efficient allocation of resources, sustainability of revenue management, infrastructure planning, and risk management (a3_buhalis2000tourism, a4_ritchie2003competitive, a5_dwyer2009destination, a6_vanhove2022economics). Bearing in mind that this type of modeling indicates market trends, consumer behavior and preferences, management structures that manage tourist destinations can use this information to identify niche markets and new trends. Therefore, forecasting tourism demand can have significant impacts on maintaining optimal competitive markets within the tourism industry (a7_crouch1999tourism, a2_frechtling2012forecasting). Based on the prediction of demand dynamics, it is also possible to adapt competitive pricing strategies, within which prices at destinations can be increased and decreased depending on the expected tourist demand (a8_song2006tourism, a9_martins2017empirical, a10_li2019competitive, a11_abrate2019impact). In addition, a12_song2012tourism argues that understanding tourist demand can influence the development of new products and services, which are compatible with the evolving needs and preferences of tourists. Tourist demand can also dictate the efficient use of marketing resources, in order to maximize reach and impact (a3_buhalis2000tourism, a13_holloway2004marketing, a14_sigala2012social, a15_hudson2017marketing), which is crucial for branding and competitiveness. Furthermore, operational efficiency is yet another factor on which tourist demand can have a significant impact. Mandal a16_mandal2018exploring states that sustainable operational efficiency within the tourism industry largely depends on data-driven decision-making, thus exploring tourist demand is also a step towards enhanced productivity optimization. This includes managing inventory, schedules, and the number of employees (a17_lenny2007impact, a18_lovelock2013strategies, a19_shabanpour2018analysis). The applicability of tourism demand modeling is especially evident when it comes to special forms of tourism affirmation (a20_burger2001practitioners, a21_trauer2006conceptualizing, a22_xie2021forecasting). In the case of nature-based tourism (a23_dimitrov2013long, a24_aliani2018modeling, a25_rice2019forecasting, a26_abu2021sarima), forecasting tourist demand can be of great importance for adjusting carrying capacity measures in certain destinations. Numerous research (a27_o1986tourism, a28_butler1999sustainable, a29_mccool2001tourism, a30_liu2003sustainable, a31_fennell2004tourism, a32_lobo2013projection, a33_zelenka2014concept, a34_lobo2015tourist, a35_guo2017remaking, a36_carrion2021environmental, a37_cheablam2021assessment, a38_sunkar2022geotourism) indicates that carrying capacity is one of the most important indicators of sustainable and responsible tourism, especially when it comes to destinations that are highly vulnerable, both from natural processes and from anthropogenic influence. Therefore, predicting the increase in tourist demand can be of great importance for management structures, because it can indicate the need to implement certain measures to prevent overexploitation and over-tourism.
In the last few decades, there has been a development of tourism of specialized interest, which focuses on geological attractiveness. Geotourism includes the affirmation of geologically significant landscapes and places that can have a certain market value obtained through the interpretation of knowledge (a39_gordon2018geoheritage). Education and conservation of geodiversity are the primary elements of geotourism and as such have the most important role in the identification and valorization of geoheritage (a40_bentivenga2019geoheritage). Therefore, geotourism through the transfer of knowledge provides value to geologically significant areas, both for the needs of tourism development (a41_dowling2006geotourism, a42_chen2015principles, a43_dowling2018geotourism, a44_olafsdottir2019geotourism) and for the effective implementation of geoconservation efforts (a45_brilha2002geoconservation, a46_gray2005geodiversity, a47_burek2008history, a48_henriques2011geoconservation, a49_crofts2020guidelines, a50_williams2020geoconservation).
In the case of karst landscapes, which represent one of the most vulnerable areas in which tourist activities are carried out (a51_ruban2018karst, a52_telbisz2020significance, a53_zhang2023aesthetic), geoconservation is a basic indicator of ethically-responsible use of karst resources a54_taheri2021human. Within the karst areas, the sites that are mostly used for mass tourism are caves (tourist caves; i.e. show caves). A detailed study on global cave tourism that explored the number of tourist visits a55_chiarini2022global indicates a very high number of visits to show caves. China boasts the highest annual visitation rate, with 19 million tourists to its cave destinations. In the United States, 9.9 million annual visitors have been recorded and within Europe, France stands out by having 5.2 million tourists annually to its caves, followed by Spain with 2.9 million visitors. Germany and Italy contribute significantly to the global cave tourism landscape, each hosting 2.4 million and 2.3 million tourists annually. Evidently, caves are a major focus of tourists around the world. Due to geoconservation standards and protection, it is necessary to pay special attention to modeling and monitoring the global tourist demand for cave tourism. Moreover, significant challenges within cave tourism are reflected primarily in the negative consequences that arise from the very arrangement of the cave for tourist use. This includes the installation of artificial lighting, construction, and introduction of substances harmful to the underground ecosystem a55_chiarini2022global. In addition, the harmfulness of tourism for caves is reflected in the increase in subterranean temperature, CO2 levels, and changes in air humidity (a56_pulido1997human, a57_baker1988environmental, a58_sebela2015cave, a59_novas2017real, a60_constantin2021monitoring). However, caves represent important destinations for multidisciplinary education, interpretation of human history, and environmental dynamics. For this reason, it is necessary to maximize the sustainable economic affirmation of caves, so that cave tourism is compatible with geoconservation standards. The advantage of management structures is that there are significant possibilities for monitoring and control within the caves themselves. In particular, visitors cannot walk outside the marked paths and cannot visit places in the cave that are not adequately lit and arranged for visiting without specialized equipment. Thus, monitoring is in most cases at a high level and this provides the possibility of effective quality control and the protection of the subterranean ecosystem.
The aim of this paper is to model tourist demand for Stopića cave in West Serbia. In previous years, this cave had an exceptional increase in the number of tourist visits, and it became the most visited, surpassing the Resava cave, which for decades was the most visited in Serbia. This unique case represents an important local economic indicator that occurred as a result of the proximity of Zlatibor, which is a highly visited mountain center. The analysis includes a comprehensive time series dataset comprising the monthly visitation figures spanning from the year 2010 through 2023, thereby encompassing a total of 168 months of observational data. This temporal scope allows an exploration of visitor trends, facilitating forecasting methodologies to be employed effectively. Through modeling of these visitation patterns, we aim to gain insights that are essential for enhancing strategic planning and management practices in the context of Stopića cave's visitor economy.
This paper is organized as follows. In section (ref) we give an overview of different approaches used in the previous studies which aimed to forecast tourist arrivals. Section (ref) presents models we use to forecast the number of visits to Stopića cave in Serbia in this study. Next, we present the experimental setup and results in section (ref). Finally, the discussion is presented in section (ref) and we conclude in section (ref).
Methods for forecasting time series can be divided into three categories: classical statistical methods, methods based on \ac{ML}, and hybrid methods which fuse both model and data-driven methodological approaches.
The classical statistical forecasting methods were exhibiting the best performances before \ac{ML} methods started outperforming them, as demonstrated in several early time-series forecasting competitions, e.g., in M3 makridakis2000m3. These methods attempt to identify patterns, trends, seasonality, and irregularities in the data observed over different time periods. They are particularly useful for understanding the underlying structure and pattern of the data and therefore offer interpretable forecasts for stakeholders. For forecasting tourism demand, the most widely used statistical forecasting method is \ac{ARIMA} and its versions which include seasonality and/or exogenous variables, see song2019review and references therein. \ac{ES} is also used in many studies that forecast tourism demand (athanasopoulos2009hierarchical, fildes2011evaluating).
In recent years, \ac{ML} techniques became popular for forecasting tourism demand, such \ac{NN} (claveria2015tourism, chen2012forecasting), \ac{SVR} (chen2007support, chen2015forecasting) and others. The most important advantage of data-driven methods is that they do not require stationarity or specific distribution of time series. Moreover, these models can explain non-linear relationships between input and output variables without a priori knowledge about them. However, the interpretability of these models is still an open research question. Also, in some applications, the amount of available data can be still too small for \ac{ML} techniques to train well so practitioners should carefully choose model complexity in order to avoid overfitting.
Hybrid methods bridge the gap between classical statistical and scalable \ac{DL} models by uniting them. Those methods are the best performers in M4 forecasting competition makridakis2020m4. In recent years, they are also used in many forecasting applications. For the purposes of tourism demand forecasting, in nor2018hybrid \ac{ARIMA} and \ac{NN} are combined in order to forecast Malaysia's tourism demand. Similarly, abellana2021hybrid combines \ac{SARIMA} and \ac{SVR} for modeling Philippine tourism demand. In this study, we consider modeling tourism demand in Stopića cave in Serbia by NeuralPropeth triebe2021neuralprophet hybrid method. As baseline methods, we use \ac{ARIMA} as the most popular statistical/classical method and \ac{SVR} - frequently utilized \ac{ML} method for tourism demand forecasting.
Although the findings from the latest M5 time-series forecasting competition makridakis2022m5 demonstrate that modern pure \ac{ML} methods based on decision trees (such as \ac{LightGBM} method ke2017lightgbm) now outperform hybrid methods, in this study, due to the limited size of the time series, we do not consider such methods due to the risk of overfitting.
Apart from forecasting tourism demand exclusively based on its previous values, it is worth mentioning that many studies investigated how exogenous variables can help in predicting targeted time series. The most popular recently studied such explanatory variables are e.g. internet big data (e.g. Google Trend sun2019forecasting, li2020forecasting, gunter2016forecasting, park2017short, volchek2019forecasting, clark2019bringing) and social media and online reviews (e.g. TripAdvisor hu2022tourism). We also consider Google Trends for modeling our time series.
In the following subsections, a brief explanation of considered time series forecasting methods is given.
Auto-Regressive Integrated Moving Average (\ac{ARIMA}) box2015time model is one of the most frequently used models in time series analysis. The model is constructed to predict future trends of non-stationary data and represents an extension of the Auto-Regressive Moving Average (ARMA) model. It can be efficiently applied to eliminate trends and the non-stationarity of the mean using differencing between consecutive observations.
\ac{ARIMA} model is generally denoted by \ac{ARIMA}$(p,d,q)$, where $p$ represents the order (number of lags) of the auto-regressive model, $d$ is the degree of differencing and $q$ denotes the order of the moving-average model. For given time series $\data_t,$ \ac{ARIMA}$(p,d,q)$ model is given by formula
where $t$ is a positive integer, $L$ is the lag operator defined as $L^i\data_t=\data_{t-i},$ $\arcoef_i$ are the coefficients of the auto-regressive part of the model, $\theta_i$ are the coefficients of the moving average part and $\epsilon_t$ are error terms. The error terms $\epsilon_t$ are assumed to be independent with normal ${\cal N}(0,\sigma)$ distributions. There are several methods for determining values of parameters $p, d,$ and $q$ such as \ac{ADF} test, \ac{ACF}, and \ac{PACF}.
As the number of tourist visits generally depends on the period of the year, for predictions of the number of tourists, it is useful to include the seasonal component in the model. Besides regular, seasonal data require seasonal differencing to become stationary. For this purpose, the \ac{SARIMA} model is used. \ac{SARIMA} model is denoted by \ac{ARIMA}$(p,d,q)(P,D,Q,M),$ where $M$ represents the seasonal period, i.e., number of observations per year, and $P, D$ and $Q$ are auto-regressive, differencing and moving average terms for the seasonal part of the model, respectively.
\ac{SARIMAX} model represents another generalization of \ac{ARIMA} model that includes both seasonality and exogenous variables. In this paper, the Google Trends data are used as one of the most popular tools in forecasting. The model has excellent performances which will be verified through results on tested data.
\ac{SVM} is \ac{ML} model initially developed for classification and later adjusted for regression (\ac{SVR}). Here we briefly introduce \ac{SVR} SVR, one of the most powerful techniques for solving both linear and nonlinear regression problems.
The linear regression model in general is given by
where $y \in \mathbb R$ is dependent variable, $x \in \mathbb R^m$ is independent variable, $\alpha \in \mathbb R^m$ and $\beta \in \mathbb R,$ are unknown coefficients and $\innerproduct{\alpha}{x}$ denotes the inner product between $\alpha$ and $x$. Classical linear regression models are based on estimating unknown coefficients for the given training set $D=\{(x_i, y_i)\}, \ i \in \{1,2,...,n\}$ by minimizing the sum of squared prediction errors (differences between the actual and the predicted values of the dependent variable). \ac{SVR} model gives us the flexibility to define how much error is “acceptable” in finding prediction values. Instead of a simple regression line (or hyperplane in high dimensional spaces), the goal here is to find a tube (Fig. (ref)) on the distance (margin) $\epsilon$ from the line ($\epsilon$-insensitive tube). In that way, the model only cares about data outside the tube. In other words, the coefficients $\alpha$ and $\beta$, which in particular describe the relationships between $y$ and $x$, are found such that the prediction errors are minimized while the margin between the regression line and the closest data points is maximized at the same time.
More concretely, the coefficients $\alpha$ and $\beta$ are in \ac{SVR} estimated by minimizing the regularized cost function under constraints:
where $\gamma$ is the balancing parameter between the regularization term of the cost function and the training error calculated as the sum of $\xi_i$ and $\xi_i^*$, which are slack variables that represent positive and negative deviations outside $[-\epsilon,\epsilon]$ region (see Fig. (ref)). In order to solve the $(\ref{eq:SVR}),$ the dual quadratic problem is formed:
where $\lambda_i$ and $\lambda_i^*$ are Lagrange multipliers that satisfy $\lambda_i\lambda_i^*=0.$ Finally, the decision function $(\ref{eq:regression})$ has the following explicit form:
where $K(x,x_i) = \innerproduct{x}{x_i}$ in the linear case. In the non-linear case, $K$ represents the kernel function that transforms the data in a higher dimensional space to be suitable for linear separation, e.g., polynomial kernel ($K(x,x_i)=\innerproduct{x}{x_i}^d$) and Gaussian ($\displaystyle K(x,x_i)=e^{-\frac{||x-x_i||^2}{2\sigma^2}}$).
When \ac{SVR} is applied for time series forecasting, the independent variable $x$ contains time series lags, and the dependent variable $y$ is the next observation in time series.
In this study we deploy NeuralProphet triebe2021neuralprophet as a hybrid time series forecasting method. It is an extension of Facebook's Prophet taylor2018forecasting, it provides information for interpreting outputs (predictions) from internal parts of the model and therefore it belongs to \ac{XAI} methods. Interpretability of NeuralProphet is achieved thanks to the fact the model is based on an additive decomposition of time series. It combines the classic time series components with scalable \ac{NN} blocks and in that way it is able to fit non-linear relationships. Two such \ac{NN} modules are the auto-regression and covariate components and thanks to them it demonstrates better predicting accuracy in comparison to Facebook's Prophet.
More formally, the NeuralProphet decomposes the time series in multiple additive components where each produces $h$ future predictions at the same time. For a single time step forecast ($h = 1$), the model is given as:
where $\trend(t)$ is the trend at time $t$, $\season(t)$ models the seasonal effects at time $t$, $\autoreg(t)$ includes the auto-regression effects at time $t$ based on past observations of the time series of interest, $\exogenouspast(t)$ captures the regression effects at time $t$ for lagged observations of exogenous variables (covariates), $\exogenousfuture(t)$ accounts for the regression effect of future-known exogenous variables at time $t$ and $\event(t)$ represents effect of certain events and holidays at time $t$. Each of the described components can be excluded if it is not relevant to the targeted time series.
The trend is modeled in a classic way, as a piece-wise linear function with the growth rate which can change at predefined number of points, so-called changepoints (model hyperparameter).
Seasonal component is modelled by Fourier terms harvey199310, with $m$ terms for seasonality with periodicity $l$:
Number of Fourier terms is by default set to be $m = 6$ with $l = 365.25$ for yearly seasonality, $m = 3$ with $l = 7$ for weekly seasonality, and $m = 6$ with $l = 1$ for daily seasonality. Mode details can be found in triebe2021neuralprophet.
\ac{AR} predicts the future values of the target variable by using a linear combination of its past values. The auto-regressive model \ac{AR}($p$) is defined as:
where $p$ is the number of linearly combined past time steps and intercept is denoted with $s$. Coefficients $\arcoef_i$ control the direction and power/significance of included past values on the future value and $\datanoise_t$ is the noise term. Classical \ac{AR} model produces only one prediction ($h = 1$). Therefore, for prediction horizon with number of steps $h > 1$, $h$ classical \ac{AR} models have to be estimated. The \ac{AR} module in NeuralProphet is based on a modification of \ac{AR}-Net triebe2019ar, which allows single model to make $h$ forecast steps for $h > 1$. Three types of auto-regression - linear, deep and sparse, can be considered within \ac{AR} module. Linear \ac{AR} is single \ac{NN} with only one layer which has $p$ inputs, $h$ outputs, and it does not have biases nor activation functions, so it is essentially same as classic statistical \ac{AR}. Deep \ac{AR} consists of a fully connected \ac{NN} with arbitrary number of hidden layers and non-linear activation functions (such as rectified linear unit (ReLU)) after each layer apart from the final one. The first layer inputs are $p$ last observations, the outputs of the final layer are $h$ future values, whereas number of hidden layers and number of neurons in them is controlled by the user. Finally, sparse \ac{AR} allows \ac{AR} order $p$ to be chosen as higher at the beginning, and then with use of a regularization only a few past observations can be forced to have weights which are not equal to $0$. It is merely a way of selecting the most significant time series lags.
The lagged regressor component is almost same as the \ac{AR} component - the only difference is that the inputs are the past values of exogenous variable instead of the targeted time series. An individual lagged regressor component has to be made for each covariate if there are multiple.
Future regressors component is same as the lagged regressor, except that we need to know the future values of exogenous variable and not only its past values.
Two types of events and holidays can be considered: user-defined events, where the user feeds the model with information about an uncommon events, or country-specific holidays, where the user only provides the name of a country and the model automatically takes into account its national holidays. In both scenarios, events and holidays are binary variables with values $1$ when the event occurs and $0$ otherwise.
In case \ac{NN} modules are deployed within NeuralPropeth, the Huber loss function during training is optimized by PyTorch optimizers where the user can define all relevant training hyperparameters such as learning rate, number of epochs, batch size, etc.
Forecast modeling of tourism demand for Stopića cave included the use of time series with the number of visitors for each month during 2010 - 2023 (168 months in total). Fig. (ref) shows how the number of visitors changed during the entire period. It can be observed that this time series has a strong seasonal component, meaning that during the summer period (July and August) when many people go on summer vacations, yearly peaks occur, whereas during winter time the number of visits is much lower in comparison to summertime. This pattern is visible during the entire period and it repeats each year. Apart from the seasonal component, that growing trend is also, especially in the second half of the considered time frame. Another interesting event to notice is the highest number of visits that happened during the summer of 2020. It was during the COVID-19 pandemic, and Serbia was in partial lockdown in that period. Many countries had traveling restrictions during that year, so vacations were spent mainly in the country of origin. This happened to domestic tourists in Serbia, they were not able to travel abroad easily during that year so they spent holidays in Serbia and it influenced that the highest number of visits in history of Stopića cave occurred at that time.
Apart from an official number of visitors, we downloaded the Google Trend\footnote{https://trends.google.com/} index for the keyword “Stopića pećina” (pećina meaning cave in Serbian) for the considered period. Since most of the tourists who visit Stopića cave are domestic tourists (more than 95 %), we opted for the keyword Serbian name of the cave. Obtained time series has also monthly frequency, and it measures the search volume of the chosen keyword. The search volume index exhibits search interest. It has values from $[0, 100]$ where value of $100$ corresponds to highest popularity of the keyword. Fig. (ref) shows the Google Trend index together with a number of visits scaled to the same range for better visibility. As can be observed from Fig. (ref), it seems that the two time series are strongly correlated having similar seasonality, trend, and peaks. It seems that many visitors to the cave were searching the name of the cave on the web, either slightly before their planned visit or at the same time. In considered models, we try to include this series as an exogenous variable.
For evaluation of considered methods, we split the time series with the number of visitors into two parts - the first 156 months (period 2010 - 2022) were used for training \ac{ML} methods and the rest 12 months (year 2023) were used for testing all methods. For the testing phase, we predict/forecast the number of visitors and compare predicted values with actual ones.
For comparison of forecasted number of visitors with real ones, we choose \ac{RMSE} defined as:
where $T$ is the size of the data used in evaluation ($T = 12$ in our case), and $\data_t$ and $\dataforecast_t$ are actual and predicted number of visitors at time $t$. The smaller the measure is, the closer real and predicted values are. When comparing predictions of different models, the model with the smallest \ac{RMSE} is considered the best.
The first model we consider is \ac{ARIMA}. Inspecting \ac{ACF} and \ac{PACF}, we concluded that the number of lags (order of auto-regressive model) that should be included in the model equals $p=3$, whereas the degree of differencing should be $d=1$ and order of moving-average $q=0$. In order to predict the entire 12 months, we fit 12 \ac{ARIMA}$(3,1,0)$ models since a single model can only predict the number of visitors for one month ahead. The plot of the actual vs. predicted number of visitors for 12 months during 2023 is given in Fig. (ref) (a), and comparing predictions with true data gave $\text{RMSE} = 4652.32$. Further, we included in the same model also seasonal component for which we use single lag $P=1$, degree of differencing $D=1$, order of moving-average $Q=0$, and $M=12$ since we have monthly data. Fig. (ref) (b) shows that including seasonal component into \ac{ARIMA} gave more accurate predictions as \ac{RMSE} significantly decreased to $\text{RMSE} = 3254.70$. Finally, we included Google Trend as an external regressor and this led to a further decrease of $\text{RMSE} = 2873.99$ and gave the best fit among all considered \ac{ARIMA} variants.
Next, we train \ac{SVR} on monthly data from 2010-2022. As already mentioned at the end of Sec. (ref), for independent variable $x$ we use time series lags. To make a fair comparison between different models, here we also consider $3$ past values of time series to be used for predicting future ones. \ac{SVR} with Radial Basis Function (RBF) kernel, regularization parameter $C=10$ and $\epsilon = 0.05$ tube is fitted, and predictions obtained with this model are presented in Fig. (ref). Computed \ac{RMSE} $4430.33$ is a little bit better than \ac{RMSE} obtained with pure \ac{ARIMA}, but it is worse than estimated \ac{SARIMA} and \ac{SARIMAX} models.
Finally, we experimented with the NeuralPropeth model. We included in the model yearly seasonality, a growing trend with the default number of trend changepoints, $3$ lags of targeted time series, and $2$ lags of Google Trend as an external lagged regressor. For modeling non-linearity, \ac{NN} with 2 hidden layers containing $4$ and $2$ nodes, respectively, is included in the model. Here we intentionally choose \ac{NN} of small size in order to prevent overfitting since the data we have has a relatively small number of observations from \ac{ML} perspective. The model is optimized with PyTorch AdamW optimizer with a learning rate of $0.003$. Fig. (ref) presents an actual and predicted number of visitors for 2023 obtained with the NeuralPropeth model. As we can see from the plot, and also by comparing NeuralPropeth \ac{RMSE} with \ac{RMSE} of previous models, the best fit is obtained by the estimated hybrid NeuralPropeth model with the chosen parameters explained above. Computed \ac{RMSE} for this model equals to $\text{RMSE} = 1726.88$ and it is approximately $40\%$ lower than the smallest \ac{ARIMA} models \ac{RMSE} (the one which \ac{SARIMAX} gave) or more than $60\%$ lower than \ac{RMSE} obtained by \ac{SVR}.
Fig. (ref) shows estimated trend and seasonal components as well as parameters for included $3$ and $2$ lags of targeted time series and Google Trend, respectively. The possibility to extract estimated model parameters is of great importance for stakeholders and policymakers, and it is an additional advantage of the NeuralPropeth model since many \ac{ML} based models are of a “black-box” nature for experts from the field of interest.
Cave tourism includes unique opportunities and challenges that require specialized strategies for effective management. The conducted analysis of the touristic demand of Stopića cave indicates the dynamism of the demand for the most visited cave in Serbia, and gives touristic implications that may be of importance to tourist organizations and decision-makers. Similar to many other destinations, Stopića cave also has visitation patterns that are influenced by seasonality, external events, and visitor preferences. The observed peak periods of visits during the summer represent the importance of adapting to seasonal growth, which includes optimizing the visitor experience. The increase in visits during the summer months of 2020 is the result of the COVID-19 pandemic, which is associated with travel restrictions, and the emphasis on the need for adaptive management for domestic tourism. The most significant increase in visits to the Stopića cave is the proximity of the mountain/tourist center Zlatibor. During the pandemic, many visitors stayed at this tourist center, which offers tourist activities throughout the season. A visit to the Stopića cave is one of the optional trips from the Zlatibor tourist center, which are often carried out as individual trips or as part of organized group excursions offered by tourist agencies. This increase in the number of visits directly affects the sustainability of Stopića cave as a tourist-accessible cave.
The sustainability of cave tourism requires a delicate balance between visitor access and conservation. As increased tourism demand can complicate carrying capacities at destinations, the results of our analysis serve as crucial preliminary inputs for the initial assessment of Stopića cave's carrying capacity. By quantifying tourism demand patterns, seasonal variations, and visitation trends, we gain valuable insights that can inform actionable steps toward carrying capacity estimates and the development of effective management strategies. However, it is essential to translate these insights into specific, applicable actions to ensure the accuracy and reliability of carrying capacity determinations and promote sustainable cave tourism management. For this, it is necessary to conduct ecological surveys and impact assessments. Nevertheless, the data shows future peak visitation levels, thus periods of potential overcrowding. This information is crucial for managing visitor access, optimizing tour routes, and implementing crowd control measures to prevent ecological degradation in sensitive areas. Recognizing these high-impact periods, tourism authorities can implement proactive management measures, such as visitor quotas, timed entry tickets, or temporary closures, to prevent environmental degradation and preserve the integrity of the subterranean ecosystem.
Furthermore, patterns in online search activity indicate a strong connection with the actual number of visitors, which further shows evident public curiosity and potential visitation intent. Exploring online search behavior can guide targeted marketing strategies and promotional efforts aimed at attracting visitors to the cave. Therefore, by gaining insights into online search, cave management can anticipate visitor trends and thus generate sustainable adaptive management strategies.
An insight into the tourist demand for the cave enables the planning of adequate measures of environmental monitoring in order to determine the dynamics of anthropogenic influence. With the increase in demand, it is necessary to establish continuous monitoring of climatic parameters that may occur due to the increased presence of visitors. Measuring changes in temperature, air humidity, and CO2 emissions are the most important factors that indicate that the increased visitation of the cave affects its ecosystem.
Increased tourist demand also indicates the importance of security measures. Thus, it is necessary to increase the safety of visitors through the implementation of adequate infrastructure, the functionality of which should be evaluated frequently. Safety also refers to the protection of the cave itself. Increased demand reflects a greater number of visitors, therefore it is necessary to introduce precautionary measures in order to maximize the protection of fragile aspects of the cave such as speleothems and groundwater quality.
Three different methods are explored for modeling tourist arrivals in Stopića cave in Serbia on a monthly basis - classical \ac{ARIMA} with and without seasonal component and Google Trend as exogenous variable, pure \ac{ML} method \ac{SVR} and hybrid NeuralPropeth method which combines classical and \ac{ML} concepts. The best fit for the chosen test period of one year is obtained with NeuralPropeth which includes the seasonal component, growing trend, non-linearity modeled by shallow \ac{NN}, and Google Trend as an exogenous variable. The estimated NeuralPropeth model apart from giving the best predictions for considered time series, it outputs also the significance of the influence of lags for both auto-regressive and exogenous variable parts, helping policymakers in that way to better understand the model and consider using it further while taking important decisions.
The obtained research results have important implications for the touristic affirmation of caves. The implementation of advanced forecasting modeling enables management structures to make various strategic moves such as sustainable management of resources and conservation efforts. The use of such analytical techniques indicates effective modeling approaches that are crucial for the sustainability and protection of subterranean karst environments. Observed trends in tourism demand and the impact of factors such as seasonality and external events point to the fact that a continuous increase in tourist visitation to Stopića cave is to be expected. Due to this prediction, it is necessary to establish adequate protection measures that can ensure long-term subterranean environmental sustainability. This primarily involves monitoring microclimate indicators such as temperature fluctuations, air humidity, and CO2 emissions. The analysis of monitoring results can greatly contribute to the understanding of the anthropogenic impact on Stopića cave, as well as the level of its vulnerability. Establishing monitoring programs and tracking visitor trends are crucial for tourism authorities so they can assess the effectiveness of carrying capacity measures and adapt management strategies in order to ensure the long-term sustainability of Stopića cave as a tourist destination.
This research has been supported by the Ministry of Science, Technological Development and Innovation (Contract No. 451-03-65/2024-03/200156) and the Faculty of Technical Sciences, University of Novi Sad through the project “Scientific and Artistic Research Work of Researchers in Teaching and Associate Positions at the Faculty of Technical Sciences, University of Novi Sad” (No. 01-3394/1). This research was also partially funded by the Provincial Secretariat for Higher Education and Scientific Research of the Autonomous Province of Vojvodina, Republic of Serbia (Grant No. 142-451-3490/2023). We would like to thank the Tourist Organization of Zlatibor and Mr. Stojan Vuković for providing the data on tourist arrivals in Stopića cave. Author Aleksandar Antić is grateful for the postdoctoral Swiss Government Excellence Scholarships for the academic year 2023/2024.