Coronavirus disease (COVID-19) is a severe ongoing novel pandemic that has officially emerged in Wuhan, China, in December 2019. As of October 13, 2020, it has affected 215 countries and territories, with over 38 million cases and approximately 1.1 million deaths (Worldometer, 2020) . At the time of writing, the most affected countries are both advanced and developing countries, such as Argentina, Brazil, Colombia, France, India, Peru, Russia, Spain, and the USA. In the last two weeks, several European countries, including Italy, saw a worrying surge of COVID-19 infections. Italy was the first European country to be severely hit by COVID-19, and it has been one of the main epicenters of the pandemic for about two months, i.e. from mid-February 2020 to mid-April 2020, when the outbreak reached the first peak. Then, the epidemic curve progressively decreased until mid-August 2020, and after that the spread of infection re-accelerated again until today. As of October 13, Italy has suffered 36,246 deaths and 365,467 cases. The likelihood of a second wave is real and makes it necessary to predict the future epidemic evolution, in plan to buy medical devices and healthcare facilities, and to manage health centers, clinics, hospitals, and ordinary and intensive care beds. Thus, the major goal of this paper is to provide short and mid-term forecasts of the patients hospitalized from COVID-19 over the period after October 13, 2020. In fact, the COVID-19related hospitalizations trends allow to have a clear picture of the overall stress and pressure on the national health care system. In particular, I implemented and compared three different time series forecast techniques and their feasible hybrid combinations: autoregressive moving average (ARIMA) model, innovations state space models for exponential smoothing (ETS), neural network autoregression (NNAR) model, ARIMA-ETS model, ARIMA-NNAR model, ETS-NNAR model, and ARIMA-ETS-NNAR model. I organized the rest of this paper as follows. In section 2, I quickly discussed the relevant literature. In section 3, I presented the data used in the analysis and discussed the empirical strategy. In section 4, I discussed the main findings and policy implications. Finally, in section 5, I provided some conclusive considerations.
From the beginning of the 2020, an increasing body of literature has attempted to forecast the spread of the COVID-2019 outbreak using different approaches (Bhardwaj, 2020; Fanelli and Piazza, 2020; Giordano et al., 2020; Nesteruk. 2020; Tuite et al., 2020; Wu et al., 2020; Xu et al., 2020; Zhao et al., 2020; Zhou et al., 2020) . The most used are ARIMA (Alzahrani et al. 2020; Benvenuto et al., 2020 , Ceylan, 2020 Perone, 2020a & b) , ETS (Bhandary et al., 2020; Cao et al., 2020; Joseph et al., 2020) , artificial neural network models (Melin et al., 2020; Wieczorek et al., 2020) , models derived from the susceptible-infected-removed (SIR) basic approach (Fanelli and Piazza, 2020; Giordano et al., 2020; Nesteruk, 2020; Wu et al., 2020; Zhou et al., 2020) , and hybrid models (Chakraborty and Ghosh, 2020; Hasan et al., 2020; Singh et al., 2020; Swaraj et al., 2020) .
The implementation and the comparison of them, except for SIR models, represent the core of this paper. Thus, in table 1, I reported 18 international studies that used single or hybrid ARIMA, ETS, and neural network to forecast the transmission patterns of COVID-19 across the world. 
The data used in this article refer to the real-time number of COVID-19 patients with mild symptoms and in ICU in Italy from February 21, 2020 to October 13, 2020, for 236 observations. I extracted the data from the official Italian Ministry of Health's website (www.salute.gov). The confirmed COVID-19-related hospitalizations trends are shown in Figure 1 . Data show that both patients hospitalized with mild symptoms and in ICU reached a first peak on April 4, 2020; after that, they followed a downward trend until mid-August, before re-accelerating again from the end of September 2020 to mid-October 2020. I compute the forecasts by using different statistical techniques and their combination. Specifically, I use a univariate linear autoregressive integrated moving average (ARIMA) models, univariate exponential smoothing state space model (ETS), univariate linear and nonlinear neural network autoregression (NNAR) models, hybrid ARIMA-ETS model, hybrid ARIMA-NNAR model, and hybrid ETS-NNAR model. The combination of different times series forecast methods should allow to maximize the chance of capturing both the linear (if any) and nonlinear epidemic patterns (Zhang, 2003; Pannigrahi and Behera, 2017) , and it is useful with phenomena like COVID-19 epidemic, which seems characterized by linear and nonlinear dynamics and components (Batista, 2020, Gupta and Pal, 2020) . As well established from the seminal work of Bates and Granger (1969) , combining techniques with unique properties may allow to achieve better performance and forecast accuracy. 2 The models are calculated as follows:
ARIMA models are detected by applying the "auto.arima" function included in the package "forecast" (in R environment) and developed by Hyndman and Khandakar (2008) . 3 This function follows sequential steps to identify the best ARIMA models, i.e. the number of p parameters of the autoregressive process (AR), the order i of differencing (I), and the number of q parameters of the moving average process (MA). 4 It combines unit root tests, 5 and the minimization of the following estimation methods: the biascorrected Akaike's information criterion (AICc), 6 and the maximum likelihood estimation (MLE). The unit root tests allow to identify the order of differencing; while the AICc and the MLE methods allow to identify the optimal parameters of the AR and MA processes. Finally, the overall goodness of fit is tested by using five common forecast accuracy measures: mean absolute error (MAE), mean absolute percentage error (MAPE), mean absolute scaled error (MASE), and root mean square error (RMSE); ETS models are identified by using the "ets" function included in the package "forecast" (in R environment), developed by . 7 ETS simple models comprise two main equations: a forecast equation and a smoothing equation. By implementing the last two equations into an innovation state space model, it is possible to get an observation/measurement equation and a transition/state equation, respectively. The first equation allows to describe the observed data; while the second equation allows to describe the behavior of the unobserved states. The states refer to the level, trend, and seasonality. In particular, I use the AICc metric for choosing the best ETS model. The goodness of fit is tested by using MAE, MAPE, MASE, and RMSE metrics; NNAR models are identified by using the "nnetar" function included in the package "caret" (R environment), developed by Hyndman, O'Hara, and Wang. 8 For non-seasonal data, I can describe NNAR models with the notation NNAR (p,k), where p denotes the number of non-seasonal lags used as inputs, and k means the number of nodes/neurons in the hidden layer. NNAR (p,k) is the same as AR process, but with nonlinear functions. The optimal number of non-seasonal lags is obtained by using the AICc metric, and the optimal number of neurons is identified by calculating (p+P+1)/2, where p is the nonseasonal AR order, and P is the seasonal AR order (if any). Finally, the goodness of fit is investigated using MAE, MAPE, MASE, and RMSE metrics; Hybrid models are identified by using the "hybridModel" function included in the package "forecastHybrid" (in R environment), developed by Shaub and Ellis. 9 To combine the single time series forecasting methods, I proceed as follows: i) first, I apply the Box-Cox (1964) power transformation to the inputs to make normality assumption more plausible; and ii) then, I implement the cross validation errors ("cv.errors"), which allow to give more weight to the models that perform relatively better, by producing the best forecast; and iii) finally, MAE, MAPE, RMSE, and Theil's U are used for validating the weighting procedure.
The estimated basic equation for ARIMA model is the following (Davidson, 2000) :
Where the second difference operator, p is the lag order of the AR process, ϕ is the coefficient of each parameter p, q is the order of the MA process, γ is the coefficient of each parameter q, and denotes the residuals of errors in time t. The estimated equations for basic ETS (A,N,N) model with additive error is the following (Hyndman and Athanasopoulos, 2018, section 7) 10 :
Smoothing equation:
Where is the new estimated level, ̂+ ‫|1‬ denotes each one-step-ahead prediction for time t+1, which results from the weighted average of all the observed data, 0 ≤ ≤ 1 is the smoothing parameter which controls the rate of decrease of the weights, and − −1 is the error at time t. 
Where −1 = ( −1, −2 … , − ) ′ is a vector containing the lagged values of the observed data, f is the neural network with n hidden neurons in a single layer, and is the error. Finally, the hybrid models are simple a combination of the described single models. 
In tables 2, 3, 4 and 5, I reported the main forecast accuracy measures, and the best selected parameters for the single and hybrid models. For the patients hospitalized with mild symptoms (Tables 2 and 3) , the optimal models are: ARIMA (4,2,4), ETS (A,Ad,N), 11 NNAR (4,2), hybrid ARIMA (4,2,4)-NNAR (7,4), hybrid ARIMA ( According to Lewis' interpretation (Lewis 1982, p. 40) , since MAPE is always lower than 10, all the predictive models are highly accurate. In particular, the forecasts show a forecast accuracy (100-MAPE) ranging from 95.68% to 97.84%. Since MASE and Theil's U are always and 11 The selected ETS model is also known as damped trend method with additive errors. 12 The selected ETS model is also known as Holt's linear trend method. 13 The coefficients of the parameters estimated for the single and hybrid ARIMA and ETS models are provided in Tables A1, A2 , A3, & A4 (Appendix A).
significantly lower than 1, all the proposed forecasting models performed better than the forecasts from the (no-change) "naïve" methods, i.e. the forecasts with no adjustments for casual factors (Hyndman and Koehler, 2006) . This allows to justify the use of more complex and sophisticated models, such as ARIMA, ETS, NNAR, and their hybrid combination. In tables 6 and 7, I compared the single and hybrid models by considering the minimization of MAE, MAPE, and RMSE. For the patients hospitalized with mild symptoms, hybrid ARIMA-ETS is better than the respective single models in 2 measures out of 6; the hybrid ARIMA-NNAR and ETS-NNAR are better than the respective single models in all the accuracy measures; and the hybrid ARIMA-ETS-NNAR is better than respective single models in 8 measures out of 9. For the patients hospitalized in ICU, hybrid ARIMA-ETS is better than the single models in 3 measures out of 6; and hybrid ARIMA-NNAR is better than the respective single models in all the accuracy measures; ETS-NNAR is better than the respective single models in 5 measures out of 6; and ARIMA-ETS-NNAR is better than respective single models in 8 measures out of 9. The most reliable and accurate single and hybrid models are ARIMA and ARIMA-NNAR, respectively. The results show that hybrid models always outperformed the single models, with the only exception of ARIMA-ETS. The best models are ARIMA-NNAR, ETS-NNAR, and ARIMA-ETS-NNAR.
In table 8, I also compute the percentage efficiency gainsin terms of MAE and RMSE minimizationfrom using the three best hybrid models. Specifically, hybrid ARIMA-NNAR outperforms single ARIMA and NNAR models from 1.87% to 8.91% on MAE, and from 4.19% to 9.47% on RMSE. Hybrid ETS-NNAR outperforms single ETS and NNAR models from 3.01% to 9.02% on MAE, and from 2.68% to 12.4% on RMSE. Finally, hybrid ARIMA-ETS-NNAR outperforms single ARIMA, ETS, and NNAR models from 0.69% to 9.61% on MAE, and from 0.1% to 12.4% on RMSE.
In figures 2 to 8, I fit all the seven models for both the time series. The light blue area shows the prediction interval at 80%, and the dark blue area shows the prediction interval at 95%. The forecasts of the best models show that the number of patients hospitalized with mild symptoms and in ICU will significantly increase in the next 45 days, i.e. from October 14, 2020 to November 27, 2020. However, if the number of patients hospitalized with mild symptom will approximately stabilize between the end of November 2020 and the beginning of December 2020, the number of patients in ICU will require more time to reach the plateau. This is also confirmed by the other forecasting methods, except of ETS. ARIMA-NNAR, ETS-NNAR, and ARIMA-ETS-NNAR show that: i) after 10 days (October 23), the number of patients hospitalized with mild symptoms will be 8, 829, 8,032, or 8,125, respectively; ii) after 20 days (November 2), they will be 14,038, 12,058, or 11,877; and iii) after 45 days (November 27), they will be 20,795, 14,823 or 16,052. About the number of patients hospitalized in ICU, ARIMA-NNAR, ETS-NNAR, and ARIMA-ETS-NNAR show that: i) after 10 days, the required intensive care beds will be 1,177, 1,040 or 1,100; ii) after 20 days, they will be 2,065, 1,446, or 1,829; and iii) after 45 days, they will be 3,320, 1,907 or 3,030. 14 Thus, a second wave of COVID-19 is expected in the next two months, with a peak in mid-December 2020. This has several policy implications both for the national health care system and economic activities. In particular, the predictions seem to stress the importance of implementing adequate containment measures and an increasing number of ordinary and intensive care beds, of hiring further healthcare personnel, and of buying care facilities, protective equipment, and ventilators to fight the infection and reduce deaths. 15 However, the opportunity of implementing more or less restrictive non-pharmaceutical interventions (NIP) to tackle the epidemicsuch as social distancing, travel ban, the use of face mask, hand hygiene, and bar and restaurant restrictions (ECDC, 2020)should be carefully evaluated because of the negative impact on the overall economic activity. In fact, according to Fitch Rating's (2020) previsions, the first wave of COVID-19 and the consequent massive lockdown measures could have already caused up to 9.5% contraction for the Italian 2020 GDP. Thus, to avoid new strict lockdown measures and further economic loss, it seems important to ensure the strict compliance with the basic COVID-19 control measures, such as social distancing, hand hygiene, and the use of protective equipment, rather than restricting or even closing public and private economic activities, such as provided by the recent Italian Prime Minister's Decrees of October 13, 2020. In fact, if the extension of the national state of emergency is well justified, the hours restriction for bakery, bar, ice cream shops, pubs, restaurant, and retail trade activities, and the total closure of all indoor and outdoor dance halls, discos, and similar spaces may seriously worsen the current Italy's economic downturn. 16 By the contrary, it may be more useful to expand public transport, adopt the double-shifts schooling system, and isolate older and vulnerable people. Finally, a better balance between public health, public utility, and freedom of economic initiative is recommended. 
In this paper, I attempted to forecast the short and mid-term dynamics of the real-time patients hospitalized from COVID-19 in Italy. In particular, I used both single time series forecast methods and hybrid combinations of them. The results show that: i) the best single and hybrid models are ARIMA and ARIMA-NNAR, respectively; ii) and hybrid ARIMA-NNAR, ETS-NNAR, and ARIMA-ETS-NNAR models outperformed the respective single models, by leading to more accurate and reliable predictions. Thus, hybrid models seem to enhance the chances of capturing a greater number of combinations of the linear and nonlinear epidemic patterns, compared with the use of single time series forecasting methods. Predictions seem also to give useful policy indications. In fact, they show that the number of patients hospitalized with mild symptoms and in ICU will significantly grow until mid-December 2020, when the second epidemic peak is expected. I predict that the necessary ordinary and intensive care beds will double in 10 days, and triple in about 20 days. Thus, it is necessary to strengthen the national health care system by buying protective equipment and hospital beds, managing health care facilities, and hiring and training healthcare workers. Finally, the hybrid combination of ARIMA, ETS, and NNAR have proven to be sufficiently accurate in the short and mid-term. However, the results in the mid-term should be taken with more caution because of the inevitable uncertainty and bias, which tend to grow over time. Notes: ***p-value < 0.01; **p-value < 0.05; *p-value < 0.1. Table A2 . Estimated parameters for ETS models (Fig.3) 