Models for tourism flash statistics integrating MNO and survey data

In this work, we rely on data provided by a single Mobile Network Operator, although several others operate within the region under examination. A brief description of how tourist overnight stays are computed from these raw MNO data is provided below; for further details on the algorithm, refer to TSS Multi-MNO project. ¹

The algorithm processes raw data to assign each user to a municipality on an hourly basis, excluding movement hours and imputing nocturnal gaps caused by device switch-offs. A valid overnight stay is recorded when a user spends at least 6 hours in a given municipality during the 20:00-08:00 window. An overnight stay is classified as tourist if the municipality falls outside the user’s Usual Environment (the set of habitually frequented locations estimated over a 12 month window using frequency and duration thresholds). The number of tourist nights is then computed as the monthly count of all such valid tourist overnight stays: if a user spends multiple consecutive nights at a destination, each night contributes independently to the total count. Finally, all counts are scaled from the single operator subscribers to the full population using proprietary expansion factors that account for market share, demographics, and device availability.

3.1. MNO data adjustment

MNO data exhibit a bias with respect to survey data. To assess whether this bias behaves as a Gaussian white noise process in log-scale, and is therefore reducible to a level correction, or whether it also presents a time-dependent structure, we analyse the log-difference

z t = log ( y t ) − log ( MNO t ) ,

(1)

We assess the white noise assumption of the bias through a set of complementary diagnostic tests: a t -test for zero mean, a Breusch-Pagan test for homoscedasticity, a Ljung-Box test for absence of autocorrelation, and a Shapiro-Wilk test for normality. The bias is considered consistent with a white noise process when none of the individual diagnostic tests is rejected at the 5% significance level. When the white noise hypothesis is not rejected according to our diagnostics, the bias is treated as purely in levels, and z t is modelled as a Gaussian white noise process with estimated mean μ ^ and variance σ ^ 2 . In this case, we apply the lognormal bias correction

MNO t * = MNO t exp ( μ ^ − σ ^ 2 2 ) .

(2)

E [ MNO ~ t ] = y t ⋅ E [ e − z t ] = y t ⋅ e − μ ^ + σ ^ 2 2 .

(3)

If the white noise hypothesis is rejected, z t exhibits time-dependent structure beyond a simple level bias. In this case, we apply a correction based on the annual lag, which is the most natural candidate given the seasonal nature of tourism data,

MNO t * = MNO t ⋅ y t − 12 MNO t − 12 .

(4)

Let y denote a response variable of interest (which in this work represents the number of nights spent by tourists according to survey data) and let x be a set of explanatory variables that become available prior to the observation of y . Within our framework, we consider the set

q = s ∪ r ,

q t = s t ∪ r t .

y ^ j = f ( x j ; s t ) , ∀ j ∉ s t .

3.3. Augmented learning

Augmented Learning (AL) techniques aim at improving predictive performance by enlarging the effective training sample through the inclusion of proxy information for units whose response variable is not yet observed. In this work, we follow the implementation of Zhang and Haug. ⁴

To address the partial unavailability of the response variable at time t , we construct an augmented sampled set s t * , defined as

s t * = s t ∪ r t * ,

3.4. Quasi-transfer learning

Transfer learning methods aim to enhance predictive performance by exploiting information from related models trained on different but comparable data sources. In this work, we adopt a quasi-transfer learning (QTL) approach, which takes advantage of auxiliary models estimated on similar populations to support the estimation of the target model.

Let f ( x ; β ) denote the target model, where β is an unknown parameter vector. We assume that an auxiliary (source) model f ( x ; θ ^ ) is available, estimated on a related population. In our application, this population coincides with the same statistical units observed at a different time point. Although the target and source models share the same functional form, they differ in their parameter values, reflecting temporal changes in the underlying data-generating process. Within this framework, the objective of transfer learning is to improve the estimation of β by exploiting the information embedded in θ ^ . The term quasi highlights that the goal is not to recover the exact target model but rather to obtain an approximation that is sufficiently close to it, as discussed in Zhang and Haug. ⁴

Let b denote a previous time period, for which complete survey data are available for all units in the population q b = s b ∪ r b , where the notation follows that introduced in Section 3.2. Since b is a past period, q b is fully observed and no proxy responses would be needed in principle. Nevertheless, we still construct the variable s b * by substituting the proxy in b for the units in r (obtaining r b * ), in order to subsequently capture its relationship with s t * , which is constructed as described in Section 3.3. The transfer mechanism is defined as

E { f ( x ; s t * ) } = g ( x , f ( x ; s b * ) ) ,

(5)

Our aim is to obtain the transfer schema described below, where the upper part represents the practical approximation obtained using the augmented sets s t * and s b * , while the lower part represents the ideal target, i.e., the models that would be obtained if the full populations q t and q b were completely observed. The downward arrow captures the idea that the quasi-transfer learning approximation should converge toward this ideal target.

f ( x , s t * ) ← g ^ ( ⋅ ) f ( x , s b * ) ⇓ f ( x , q t ) ← g ^ ( ⋅ ) f ( x , q b )

(6)

In the proposed implementation, the transfer function g ( ⋅ ) incorporates the prediction f ( x ; s b * ) as an additional feature in the target model: the source model is first fitted on s b * , its predictions are then computed for all units in s t * , and finally these predictions are included as an additional covariate when fitting the target model g ( ⋅ ) on s t * . In parallel, a second source model is fitted on q b , and its predictions are computed for all units in r t = q t ∖ s t , where they are included as an additional covariate when applying g ( ⋅ ) to produce the final flash estimates. The observed survey values are retained for the early respondents s t .

In summary, the models f ( x ; s b * ) and f ( x ; q b ) act as source models, whereas f ( x ; q t ) represents the target model. Here, q b denotes the full population at the past period b , which is completely observed, while q t denotes the full population at the current period t , which is the quantity we aim to approximate. All models share the same functional form but are estimated on different data subsets, resulting in distinct parameter estimates. g ( ⋅ ) differs from f only in the presence of the additional feature, and the variables chosen are described in Section 4.

Alternative substitution strategies for missing responses, as well as different transfer schemes and transfer functions, may be considered depending on the application; a comprehensive discussion is provided in Zhang and Haug. ⁴ This approach allows us to exploit structural similarities across populations while accommodating differences in their underlying distributions.

Statistical learning methods are typically developed within an inductive framework, in which a model is trained on a set of labeled units and then generalized to unseen observations. Transductive inference, ⁶ on the other hand, operates on a fixed and finite population of units with partially observed labels, where the target units are known at training time. Rather than estimating a general predictive function, transduction exploits the known characteristics of the target units directly, making inference inherently tied to the specific population under study. AL and QTL methods belong to this category: the availability of proxy information for the dependent variable allows units with unobserved survey outcomes to be included in the training process. This distinguishes our framework from standard cross-sectional prediction, where out-of-sample units are genuinely unobserved during training.

In this sense, the augmented training set s t * consists of an inductive component s t and a transductive component r t * . This hybrid structure places our approach within what might be called a hybrid inductive-transductive framework.

3.6. Testing framework

To obtain robust and comprehensive performance metrics across different configurations of early respondents, we adopt a repeated n cross-validation scheme. The purpose of this procedure is to conduct an a posteriori simulation study in all municipalities, taking advantage of the fact that the target variable y is fully observed for each unit at the evaluation time.

When partitioning the dataset into folds, units in the labeled folds are treated as early respondents and constitute the inductive component of the augmented training set s t * . Units in the held-out fold, by contrast, are treated as simulated non-respondents: their survey outcome is masked and replaced by the available proxy, so that they contribute to the training process exclusively through proxy information, constituting the transductive component r t * . Since each n -fold cross-validation rotates the role of the folds, every municipality acts as a simulated non-respondent at least once, and a forecast is produced for every unit in the region at the end of each run.

To reduce variability due to the random allocation of municipalities across folds, the procedure is repeated 50 times, each time with a different fold assignment. For each municipality i and fixed time t , a stable prediction error is obtained by averaging the absolute percentage errors (APE) across the 50 repeated runs, defined as

A P E ¯ i , t = 1 50 ∑ r u n = 1 50 | y i , t − y ^ i , t ( r u n ) | y i , t ⋅ 100.

MAPE t = 1 N ∑ i A P E ¯ i , t ,

A E ¯ i , t = 1 50 ∑ r u n = 1 50 | y i , t − y ^ i , t ( r u n ) | ,

To account for heterogeneity in the scale of overnight stays, prediction performance is also evaluated using the weighted absolute percentage error (WAPE), obtained by weighting the averaged APEs by the observed number of overnight stays:

WAPE t = ∑ i y i , t A P E ¯ i , t ∑ i y i , t .

The labeled and held-out folds in the cross-validation were constructed to reflect the desired proportion of early respondent units. For example, to obtain a share of 10% early respondents, we adopted a 10-fold cross-validation scheme in which, at each iteration, one fold is treated as the labeled set (early respondents), while the remaining nine folds are assigned to the held-out set, with the target variable replaced by the proxy. To obtain a 50% rate of early respondents, we use a two-fold cross-validation scheme, with one fold composed of early respondent units and the other consisting of units for which the target variable was replaced by the proxy and used as the held-out set. In each iteration, predictions and error computations are restricted to units with unobserved survey outcomes. Overall, the aggregate error values obtained through the proposed simulation framework can be considered pessimistic, as in practice learning-based estimates would be released together with the observed data from early respondents, which are unaffected by prediction error.

In the limiting case where no early respondent units are available, inference would rely exclusively on proxy-based information observed on the full population. No inductive generalization would be involved, as predictions would be restricted to the same finite set of units entering the learning process, resulting in a fully transductive inference problem.

3.7. Univariate time series model: TRAMO

The TRAMO (Time series Regression with Arima noise, Missing observations, and Outliers) model ⁷ is a widely used univariate time series methodology for preprocessing in seasonal adjustment and for forecasting. In this work, we use its implementation in the rjd3tramoseats R package. ⁸

TRAMO models the observed series through a seasonal ARIMA (SARIMA) specification, ⁹ which extends the standard ARIMA framework by including additional seasonal autoregressive and moving average components operating at lag s . In our context, s = 12 to account for the monthly seasonality of tourism data. The SARIMA ( p , d , q ) ( P , D , Q ) 12 model can be written compactly as

Φ ( B 12 ) ϕ ( B ) Δ d Δ 12 D y t = Θ ( B 12 ) θ ( B ) ε t ,

Model identification, outlier detection, and parameter estimation are handled automatically by the RSAfull procedure implemented in TRAMO, which simultaneously tests for the need of a log-transformation, accounts for calendar effects (Easter, trading day, and leap year), selects the appropriate SARIMA orders, and detects structural anomalies such as additive outliers, temporary changes, and level shifts. No external regressors are included in the models estimated in our experiments.

The choice of a univariate time series model for this work fell on TRAMO due to its simplicity in estimation, manipulation, and interpretation, combined with the level of detail it provides in modeling (including automatic specification).

4.1. Experiment 1: Learning-based approaches and their variants

All analyses presented here were carried out on municipalities with a number of overnight stays above the median value computed across all municipalities and all available months (577 nights), resulting in a subset of 175 units. This selection excluded municipalities without a touristic vocation, which would mainly introduce noise into the estimates and were less likely to benefit from the production of flash indicators.

Rather than defining a strict a priori target, the selection process followed an initial exploratory phase. In the preliminary explorations, we observed that predicting the tourist’s overnight stays for municipalities with low tourism volumes was highly volatile and challenging. Consequently, we chose to focus our analysis on the most significant ones to ensure robust results.

At this stage, we fixed the proportion of early respondent municipalities to 50% of the units by adopting the two-fold repeated cross-validation testing framework, applied to multiple model variants for each approach. The observed share of early respondents is not systematically recorded by our institute at present. However, it is known that data from some municipalities could be made available earlier by reorganising the data collection process appropriately. In the absence of a precise estimate of this share, the 50% threshold was adopted as the baseline configuration. This threshold also corresponds to the share beyond which the MAPE of the baseline approach no longer improves with the addition of further early respondent units, as shown in the scalability analysis of Section 5.3. As proxy variables, we tested the number of overnight stays derived from MNO data in its adjusted form ( MNO t * ) and the survey-based value from the previous year ( y t − 12 ). When a variable was used as a proxy, it was excluded from the set of model features in order to avoid using the same information to predict its own value.

Regarding the features, we always included the one-month lagged survey variable ( y t − 1 ), combined with several alternatives: the adjusted MNO variable ( MNO t * ), the unadjusted MNO variable ( MNO t ), and the pair composed of MNO t and MNO t − 12 . Due to the limited length of the available MNO time series, it was not possible to compute MNO t − 12 * .We also tested the inclusion of the tourism municipality typology defined by Regione Emilia-Romagna ⁵ as an additional feature. However, preliminary experiments showed no improvement in predictive performance, suggesting that the other predictors already captures the information associated with tourism profiles.

It is worth noting that the inclusion of the lagged survey variable y t − 1 , which is expected to have a strong impact on short-term forecasting performance, implicitly assumes that survey data for the previous month are available at the time of estimation. This requirement represents the main operational bottleneck of the system. As a possible workaround, when y t − 1 is not available, one could consider replacing it with an MNO lagged variable. This strategy would allow estimates to be anticipated to the release time of MNO and early respondents data, for example about 15 days after the end of the reference month, rather than one month later. To evaluate learning-based approaches that can operate in such an earlier setting, we specified additional model variants including the two-month lagged survey variable ( y t − 2 ) and the one-month lagged MNO variable ( MNO t − 1 ) for short-term adjustment. The current adjusted MNO value was also included as a feature, while y t − 12 was used as the proxy variable.

In a fully transductive inference scenario, capable of operating without early respondent units, the learning-based models described above could produce a flash estimate as early as one day after the end of the target month, provided that MNO data are available in real time.

4.2. Experiment 2: Learning-based vs time series approaches

This experiment focused on comparing the adjusted MNO proxy, the learning-based approaches, and predictions obtained from univariate time series models estimated solely on survey data using TRAMO. The latter, hereafter referred to as the “time series approach”, relied on models fitted to approximately nine years of monthly survey data (starting in January 2015) for a limited set of ten municipalities.

In contrast, as in the previous experiment, the AL and QTL approaches were trained on much shorter time series, consisting of 12, 13 and 14 monthly observations for August, September and October, respectively, while exploiting cross-sectional information from a much larger set of 175 municipalities. This comparison allowed us to assess whether adjustment and learning-based methods leveraging rich cross-sectional information could compensate for the limited temporal depth of the available data. ARIMA-TRAMO models therefore provided a natural benchmark for evaluating the performance of more flexible learning-based methods in contexts where long time series are available.

MAE, MAPE, and WAPE for the learning-based approaches were derived from the A E ¯ i , t of the 10 municipalities from Experiment 1, averaged over the newly considered set.

The baseline time series specification, whose results are denoted by y ^ t + 1 (forecast horizon h = 1 ), corresponds to a TRAMO model estimated on data including the observation of the month preceding the target, that is, assuming full availability of y t − 1 . Several variants designed for the absence of lag-1 terms were also evaluated. First, a pure time series application with a longer forecast horizon ( h = 2 ) was tested, denoted as y ^ t + 2 . Second, an “augmented time series” approach (forecast horizon h = 1 ) was introduced, with y t − 1 replaced by MNO t − 1 * only in the prediction phase, while the model is still trained on survey data. The resulting forecasts are denoted by y ^ t + 1 * . Since the rjd3tramoseats implementation does not allow separating the estimation and forecasting phases for conditional inputs, this hybrid augmented variant was implemented using the tsibble ¹³ and fable ¹⁴ packages.

4.3. Experiment 3: Scalability analysis

A scalability analysis was performed for the baseline, AL, and QTL approaches by selecting a specific model variant from Experiment 1, namely the one using y t − 12 as proxy and { y t − 1 , MNO t * } as features. The 50-times-repeated n -fold cross-validation was then repeated under different partitions of the folds between training and test sets, yielding the following percentages of early respondents: 10%, 25%, 33%, 50%, 66%, 75%, and 90%. The MAPE was used as the performance evaluation metric.

Results for September and October 2023 are reported in Tables 2 and 3. In each table, the best-performing variant within each approach (Baseline, AL, and QTL) is underlined, while the overall best value for each metric across all approaches is highlighted in bold. When the two coincide, both formatting conventions are applied. In general, WAPE was consistently lower than MAPE across all approaches and variants, indicating that prediction accuracy tended to improve in municipalities with more intense tourist activity. Excluding cases in which y t − 2 was used instead of y t − 1 for short-term adjustment, a striking result is that the best-performing variant, in terms of MAPE, of the baseline approach performed worse than all AL and QTL variants, with only two really marginal exceptions. In particular, for October, the baseline specification using the adjusted MNO proxy ( MNO t * ) achieved a MAPE of 34.9%, virtually identical to the corresponding AL result (35.0%), and comparable to the third QTL variant (38.3%), the latter being the lowest-performing among those incorporating y t − 1 as a feature.

OPEN IN VIEWER

Table 3.

Performance indicators for different approaches and model variants, including the adjusted MNO proxy ( MNO * ), October 2023.

Approach	Proxy	Features	MAPE	WAPE	MAE	Tot. AE
MNO *	MNO t *		38.0	20.6	2 311	369 781
Baseline		{ y t − 1 , y t − 12 }	34.9	28.1	5 054	808 695
		{ y t − 1 , MNO t * }	38.6	28.6	5 137	821 993
		{ y t − 1 , MNO t }	46.7	31.9	5 738	918 050
		{ y t − 1 , MNO t , MNO t − 12 }	49.0	31.5	5 662	905 985
		{ y t − 2 , MNO t * }	66.1	31.0	5 565	890 389
		{ MNO t − 1 , MNO t * }	42.8	25.1	4 517	722 660
AL	MNO t *	{ y t − 1 , y t − 12 }	35.0	14.0	2 509	401 490
	y t − 12	{ y t − 1 , MNO t * }	25.1	11.0	1 974	315 796
	y t − 12	{ y t − 1 , MNO t }	28.4	12.8	2 301	368 208
	y t − 12	{ y t − 1 , MNO t , MNO t − 12 }	27.7	10.5	1 880	300 736
	y t − 12	{ y t − 2 , MNO t * }	44.2	18.8	3 387	541 869
	y t − 12	{ MNO t − 1 , MNO t * }	27.2	10.2	1 837	293 913
QTL	MNO t *	{ y t − 1 , y t − 12 }	30.0	11.7	2 110	337 603
	y t − 12	{ y t − 1 , MNO t * }	30.3	12.5	2 255	360 811
	y t − 12	{ y t − 1 , MNO t }	38.3	21.6	3 879	620 611
	y t − 12	{ y t − 1 , MNO t , MNO t − 12 }	31.9	14.2	2 558	409 282
	y t − 12	{ y t − 2 , MNO t * }	50.9	30.7	5 528	884 417
	y t − 12	{ MNO t − 1 , MNO t * }	33.1	14.2	2 552	408 278
Total			y tot			1 797 529

The performance of MNO * used directly as a flash estimate was also analyzed. Its MAPE was of the same order of magnitude as that of the baseline variants, while outperforming most of them. It achieved a substantially lower WAPE, only slightly higher than that of the learning-based approaches, indicating particularly good performance in municipalities with a larger number of overnight stays. This situation highlighted the need for a model-based approach to fully exploit MNO data, especially in municipalities with lower values of the target variable, where noise played a more substantial role. In September, the best performance of the baseline approach was obtained when using the adjusted MNO as feature, together with y t − 2 , yielding a MAPE of 28.9%. By contrast, the advanced approaches reached substantially lower errors, with the best MAPE equal to 21.1% for the AL model using y t − 12 as proxy and the adjusted MNO among the features, together with y t − 12 . The AL model configuration using y t − 12 as proxy and { y t − 1 , MNO t * } as features was optimal in October in terms of MAPE, albeit with a wider performance gap. Remaining on October, the best baseline specification reached a MAPE of 34.9%, whereas the top-performing AL model attained a substantially lower value of 25.1%. Differences were even more pronounced when considering WAPE. In September, the best WAPE for the baseline approach was 15.8%, which was substantially higher than the values achieved by the advanced methods, with the worst AL specification still below 9% and the less precise QTL specification reaching 10.8%. This gap further widened in October, with the lowest WAPE for the Baseline approach reaching 25.1%, compared with 10.2% and 11.7% for the best-performing AL and QTL variants, respectively.

In summary, considering both MAPE and WAPE, AL was the best-performing approach. The Total Absolute Error, which is related to both metrics, confirmed this result. A strong model variant for forecasting, in terms of both MAPE and WAPE, in September and October, was the one using y t − 12 as a proxy and { y t − 1 , MNO t * } as features. This variant was therefore selected for further investigation in Experiments 2 and 3. Considering the variants capable of producing estimates with only a 15-day delay (the last two variants for each approach), the relative performance differed between the analyzed months. In October, variants using MNO t − 1 clearly outperformed those relying on y t − 2 , whereas the opposite was observed in September. Overall, this seemed to favor the use of MNO t − 1 .

These tests provided an example of how timeliness requirements and data availability can influence the choice of model specifications. It was noted that, in these examples, although adjusting MNO data was essential when they were used as proxies, this adjustment was not necessarily beneficial when MNO variables were employed as model features. This may have been due to the use of random forests, which are non-linear models, whereas the same effect might not hold for linear models.

Figure 3 shows the distribution of MAPE across municipalities for September 2023. To ensure a fair comparison, the MAPE scale was clipped at 50% in all scenarios. The maps clearly show that, in coastal municipalities and major urban areas, which were characterized by higher levels of tourist activity, both AL and QTL outperformed the baseline approach. In contrast, municipalities located in Apennine areas exhibited higher errors across all approaches.

OPEN IN VIEWER

Figure 3.

Spatial distribution of Mean Absolute Percentage Error (MAPE), September 2023. MNO * and model variants with features: { y t − 1 , MNO t * } , proxy: y t − 12 , 50% of early respondent units. (a) Adjusted MNO data, (b) Baseline approach and (c) Augmented learning and (d) Quasi-transfer learning

The results for September and October indicate that AL and QTL exhibit comparable performances, especially for some variants. Due to the shortness of the available survey time series, it was not possible to apply the QTL approach for August: doing so would have required observations of y at the past time point b , and even considering b equal to July 2023 (the month preceding the target), the implemented variants would need y for July 2022 ( y t − 12 ), which falls outside our observation window. Considering this, we realized the experiment for August without the QTL, but we can hypothesize that it would reach similar results to the AL. We report in Table 4 the error metrics for the baseline and AL approaches only. The results obtained are consistent across all model variants, allowing us to focus on one specific configuration and generalize the findings. We select the variant utilizing y t − 12 as the proxy and { y t − 1 , MNO t * } as features, which achieves the best MAPE and exhibits a similar behavioral pattern in both AL and QTL for the subsequent months. The results are highly significant and reveal that predicting August appears less challenging than September and October; this holds true even for the baseline model, which records a MAPE and WAPE of 27.5% and 14.8%, respectively. Particularly noteworthy is the capacity of the adjusted MNO data to achieve highly precise approximations of the target variable in municipalities with a high volume of tourist overnight stays, yielding a WAPE of 5.9% against a MAPE of 35.5%. The AL approach further enhances these results, lowering the metrics to a 20.5% MAPE and a 4.6% WAPE. However, in terms of absolute errors, the adjusted MNO data still exhibits a lower MAE (and, consequently, a lower Total AE). This implies that while AL smooths the predictions to optimize performance across the entire set of municipalities, it slightly sacrifices its accuracy on the largest destinations.

OPEN IN VIEWER

Table 4.

Performance indicators for different approaches and model variants, including the adjusted MNO proxy ( MNO * ), August 2023.

Approach	Proxy	Features	MAPE	WAPE	MAE	Tot. AE
MNO *	MNO t *		35.5	5.9	2 807	516 420
Baseline		{ y t − 1 , y t − 12 }	21.3	15.5	13 560	2 495 085
		{ y t − 1 , MNO t * }	27.5	14.8	12 966	2 385 739
		{ y t − 1 , MNO t }	38.2	16.7	14 634	2 692 730
		{ y t − 1 , MNO t , MNO t − 12 }	45.2	15.6	13 646	2 510 932
		{ y t − 2 , MNO t * }	34.3	15.2	13 283	2 444 064
		{ MNO t − 1 , MNO t * }	41.3	15.0	13 108	2 411 916
AL	MNO t *	{ y t − 1 , y t − 12 }	26.3	4.0	3 542	651 664
	y t − 12	{ y t − 1 , MNO t * }	20.5	4.6	3 980	732 351
	y t − 12	{ y t − 1 , MNO t }	25.5	3.3	2 858	525 857
	y t − 12	{ y t − 1 , MNO t , MNO t − 12 }	28.8	6.2	5 390	991 711
	y t − 12	{ y t − 2 , MNO t * }	23.1	4.8	4 222	776 836
	y t − 12	{ MNO t − 1 , MNO t * }	27.3	5.3	4 658	856 984
Total			y tot			8 755 264

Considering the quality of the results obtained using the adjusted MNO data as predictors, we also evaluated the performance of the raw MNO data instead. While the resulting WAPE was a reasonable 26%, the MAPE skyrocketed to around 330%, confirming that the adjustment procedure described in Section 3.1 is strictly necessary to make these data operationally useful when used as a standalone proxy.

Taken together, all these analyses suggest that the MNO data, when properly treated, are more effective in describing municipalities with high touristic activity, i.e., where the target value they aim to approximate exhibits a higher magnitude. In fact, the WAPE is consistently lower than the MAPE across all the considered months. This is particularly evident for August, the month with the highest tourist volumes, where the WAPE drops drastically and the adjusted MNO data alone achieve a lower Total AE than any other modeling approach.

The scenarios reported in Table 5 assume that survey data at t − 1 relative to the target month were available, whereas Table 6 presents the solutions adopted when such data were not. The tables report MAE, MAPE, and WAPE (computed by weighting errors only over the ten selected municipalities) for the different approaches applied to this subset.

OPEN IN VIEWER

Table 6.

Accuracy measures for flash estimates of tourist overnight stays.

	August			September			October
Approach	MAE	MAPE	WAPE	MAE	MAPE	WAPE	MAE	MAPE	WAPE
Adjusted MNO proxy (MNO * )	10 296	8.5	2.7	17 774	18.3	9.1	16 637	25.4	15.5
Augmented learning (AL)	54 098	16.6	10.6	22 773	13.6	11.6	9 678	13.3	9.0
Quasi-transfer learning (QTL)	–	–	–	22 629	10.8	11.5	19 027	17.8	17.7
Time Series: y ^ t + 2	12 150	6.0	3.3	21 067	9.0	11.0	10 734	11.5	10.7
Augmented Time Series: y ^ t + 1 *	305 647	39.8	81.8	121 501	47.1	63.2	38 460	41.0	38.3

Comparison of time series models, learning-based approaches, and the adjusted MNO proxy for August, September, and October in selected municipalities, assuming that data for the month preceding the target period are not available. Learning-based approaches use features { MNO t − 1 , MNO t * } and proxy y t − 12 .

If we consider all the accuracy measures together across the three months analyzed, the time series approach achieved the best performance. The time series approach did not significantly prevail in August when y t − 1 was available, as the adjusted MNO proxy reached the exact same WAPE as the time series model (2.7%), albeit with a slightly higher MAPE. When the t − 1 data were not available, the Adjusted MNO proxy registered a lower MAE and WAPE than the TRAMO variant with forecast horizon h = 2 ( y ^ t + 2 ) in August and September, while AL achieved the best performance in October, although the time series approach always maintained the lowest MAPE. Nevertheless, despite being slightly higher than those of MNO * and AL in these specific cases, the MAE and WAPE of the time series approach remain close, rendering the performance of these approaches remarkably similar.

Another consideration is that the “augmented” time series variant y ^ t + 1 * , which replaces the unavailable y t − 1 with MNO t − 1 * , did not perform satisfactorily.

Although the traditional time series approach generally appears to prevail when y t − 1 is not available, a municipality-level analysis reveals a richer picture. Tables 9, 7 and 8 show that, in several cases, the learning-based approaches or the adjusted MNO ( MNO * ) substantially outperformed the time series forecasts in terms of unit-level Absolute Percentage Error (APE).

OPEN IN VIEWER

Table 7.

Absolute percentage error (APE) by municipality and approach, including the adjusted MNO proxy (MNO * ).

Municipality	y	APE Time series y ^ t + 1	APE AL	APE QTL	APE MNO *	Best approach
Alto Reno Terme	4 470	16.0	18.1	16.4	34.9	Time series
Forlì	25 635	13.0	13.4	14.2	3.1	MNO *
Piacenza	34 877	5.2	9.1	10.7	28.1	Time Series
Ferrara	45 242	11.9	14.1	4.5	40.7	QTL
Modena	65 063	9.9	3.8	9.0	32.2	AL
Parma	86 043	4.6	16.2	16.5	18.4	Time Series
Ravenna	293 864	12.8	4.5	0.4	0.1	MNO *
Bologna	345 342	3.4	10.8	12.4	16.0	Time Series
Riccione	369 663	11.6	7.9	2.3	2.8	QTL
Rimini	692 827	2.0	22.4	20.8	6.5	Time Series

Learning-based approaches use features { y t − 1 , MNO t * } and proxy y t − 12 . September 2023.

OPEN IN VIEWER

Table 8.

Absolute percentage error (APE) by municipality and approach, including the adjusted MNO proxy (MNO * ).

Municipality	y	APE Time series y ^ t + 1	APE AL	APE QTL	APE MNO *	Best approach
Alto Reno Terme	2 891	2.4	16.1	5.4	37.9	Time Series
Forlì	31 495	0.5	22.8	20.7	6.7	Time Series
Piacenza	33 386	3.5	1.8	2.3	30.1	AL
Ferrara	46 868	15.2	23.3	29.1	68.9	Time Series
Modena	70 430	3.0	5.3	1.5	24.5	QTL
Ravenna	81 862	9.0	8.2	2.7	15.8	QTL
Riccione	83 551	18.8	10.8	23.1	16.5	AL
Parma	84 450	10.8	0.8	0.2	42.9	QTL
Rimini	255 793	4.2	21.1	32.7	1.5	MNO *
Bologna	382 827	0.7	12.7	3.4	9.7	Time Series

Learning-based approaches use features { y t − 1 , MNO t * } and proxy y t − 12 . October 2023.

OPEN IN VIEWER

Table 9.

Absolute percentage error (APE) by municipality and approach, including the adjusted MNO proxy (MNO * ).

Municipality	y	APE Time series y ^ t + 1	APE AL	APE MNO *	Best approach
Alto Reno Terme	4 251	3.8	17.4	16.4	Time Series
Forlì	21 543	6.8	3.8	11.5	AL
Piacenza	25 739	6.7	7.8	11.8	Time Series
Ferrara	38 344	12.0	23.0	14.9	Time Series
Modena	46 540	13.9	22.6	20.6	Time Series
Parma	64 554	2.3	4.9	1.6	MNO *
Bologna	303 727	4.4	2.8	1.8	MNO *
Ravenna	789 693	2.1	3.6	1.0	MNO *
Riccione	885 171	0.6	5.4	3.1	Time Series
Rimini	1 650 831	3.1	15.6	2.4	MNO *

Learning-based approaches use features { y t − 1 , MNO t * } and proxy y t − 12 . August 2023.

This occurred in eleven municipalities across September and October, typically when the Absolute Percentage Error of the time series approach exceeded 10%. The only notable exceptions were Ferrara in October and Alto Reno Terme in September (where QTL exhibited performance comparable to TRAMO). For instance, in Ravenna in September, the time series approach yielded an APE of 12.8%, whereas QTL achieved a very small error of 0.4% and the adjusted MNO reached an even lower value of 0.1%. Similarly, in October in Riccione, the time series approach recorded a MAPE of 18.8%, compared with 10.8% for the AL approach. Although there were no municipalities in which a single approach consistently dominated the others in both September and October, these results opened up perspectives for hybrid solutions. In particular, when the time series approach struggles, for instance during periods of high uncertainty, MNO-based solutions, which incorporate near-real-time information, may represent effective alternatives. In August, when tourist overnight stays reach their overall peak, the adjusted MNO proxy always prevails in the most visited municipalities (where the number of overnight stays surpasses 50 000), with the sole exception of Riccione (despite its excellent APE of 3.1%). In these high-volume destinations, the AL approach also performs remarkably well.

Taken together, these findings further confirmed the informative value of MNO data, which could be exploited as external regressors in time series models such as TRAMO when longer MNO time series become available.

It is worth noting that a limitation of the present comparison between learning-based and time series approaches is its restriction to ten municipalities. In our case, the difficulty in retrieving long time series for all the municipalities is due to operational constraints linked to the volume of data involved and to changes in the administrative boundaries of municipalities over time, such as mergers and splits. This limitation, however, further highlights one of the advantages of the proposed MNO-based learning approaches when estimates are needed for a large set of municipalities. By combining MNO data with cross-sectional information, they can produce competitive estimates for all the municipalities while requiring a considerably smaller volume of historical data. This suggests that learning-based methods could represent a simpler yet effective alternative to traditional approaches in contexts where long time series are difficult to obtain or to handle, as well as a complementary source of estimates to be considered alongside traditional forecasts.

From the scalability analysis described in Section 4.3, we observed that MAPE remained stable for both AL and QTL approaches, while it deteriorated for the baseline approach, as shown in Figure 4. The results confirmed the effectiveness of the learning-based approaches and highlighted the relevance of the information contained in the proxy variables. Moreover, since predictive performance did not degrade even when only 10% of units provided early responses, these findings opened up perspectives for fully transductive learning solutions, where predictions could be based on proxy-dependent variables only, without relying on early respondent observations.

OPEN IN VIEWER

Figure 4.

Performance of baseline approach, AL and QTL by varying the number of early respondents. Features: { y t − 1 , MNO t * } , proxy: y t − 12 .