At the beginning of 2020, there happened to be an out break of pneumonia, started from an unknown aetiology in Wuhan City, Hubei Province, China. A novel strain of coronavirus (COVID-19) was found from a cluster of patients in January 2020 [22] . Some scientists proposed a hypothesis that this strain of virus might have originated from bats that could be traced back to Huanan South China Seafood Market, a seafood market [43] . There has been evidence about person-to-person transmission since the middle of December 2019 [22] .
The number of cases has increased exponentially since then and authorities all over the world have to promulgate laws to alleviate the widespread of disease. Europe made the strict act as to close borders of all Schengen area [5] . At the end of the first quarter, USA surpassed China to become the world's most infected territory at more than 100, 000 confirmed cases [11] . Roughly a month later, USA, again, became the first country to surpass 1 million confirm cases, however, the actually cases could be 10 times higher due to the lack of test kits [20, 33] . Most lately, after half a year of outbreak, the total confirm cases all over the world has reached 10 million, with more than 500, 000 people died. No country has succeeded in producing vaccination or special drugs, while the current situation is extremely complex and unpredictable, the need of data modelling and prediction is a must.
The number of COVID-19 infected cases can be regarded as a time series problem. According to previous research, time-varying process forecasting is one of the most important fields in terms of statistics. By collecting and analyzing different data points, researchers could investigate the situation and surmount current problem. However, most of these problems require a significant amount of data, which might be difficult due to the lack of information, or the ambiguity and non-linearity of data. For example, data of exchange rate between two currencies may fluctuate throughout a day or a period of time. Essentially, there is an urge to develop a method to solve this real-world dilemma.
The remainder of this paper will be organized as follows. In Section 2, we will provide a brief description of the traditional Grey models, including GM(1, 1) and the Nonlinear Grey Bernoulli Model GM (1, 1) . Section 3 briefly introduces methods of optimizing parameters simultaneously in the formation of the optimized ONGBM(1, 1) model. In section 4, the application of the rolling mechanism is presented to propose a new RONGBM(1, 1) model. Section 5 provides some statistical and modern machine learning models as a benchmarking tool for the proposed model in the previous section. Section 6 introduces some modern machine learning models as an additional benchmarking tool for the newly proposed model in the empirical examples. Section 7 adopts the example of predicting Vietnam's GDP to demonstrate the accuracy of the new RONGBM(1, 1) model in comparison with models in Section 5, then applies this new model into predicting the total COVID-19 infected cases until 2020-02-18. Section 8 concludes the paper by drawing some conclusions and future work.
2 Description of the traditional GM(1, 1) and NGBM(1, 1) models In this section, a brief introduction of the mathematical modelling of two traditional models: GM(1, 1) and Nonlinear Grey Bernoulli Model NGBM(1, 1) is introduced.
Assume that x (0) is the non -negative original historical time series of data with m entries as:
Next, we define x (1) using one-time accumulated generating operation (1 -AGO), which is as
As x (1) is a monotonic increasing sequence which is similar to the solution of a first-order linear differential equation, one can assume that the solution of the following differential equation
represents the Grey predicted value complement to the initial conditionx (1) (1) = x (0) (1) and parameters a and b.
By definition, dx dt can be written as:
However, to discretize the differential equation, ∆t can be set to be equal to 1, which means that (5) can be rewritten as
The Grey predicted value is now defined aŝ
with P usually set as 1 2 in the traditional models. The differential equation can now be discretized as:
By the least-squared method, the coefficients a and b can be determined by
...
The particular solution of equation (4) with the initial condition is:
The prediction of the historical time series of data at point k + 1 can now be deduced by:
with
In order to obtain higher accuracy in predicting comparing to the original GM(1, 1) model, Professor Chen [4] proposed the Nonlinear Bernoulli Grey Model NBGM(1, 1) as follows:
Similarly to the traditional Grey model GM(1, 1), assume that x (0) is the non -negative original historical time series of data with m entries as:
Next, we define x (1) using one-time accumulated generating operation (1 -AGO), which is as
where
x (0) (k), k = 2, 3, ..., m
As indicated previously, equation (4) is a linear differential equation. A similar form of this equation, which is nonlinear and has the form of:
where n ∈ R any real number is called a Bernoulli equation, or the with then differential equation of the NGBM(1, 1) model.
The background value is now also defined as:
with P = 1 2 for the traditional model. Then, discretizing the ODE, one obtains
which is called the basic Grey differential equation of the NGBM(1, 1) model. It can easily be recognized that for n = 0, this equation turns to equation (4) , which is the traditional GM(1, 1) model; for n = 2, the equation turns to the Grey-Verhulst equation.
By the least square method, the parameters a and b can be determined by:
...
The particular solution of equation (17) , or the discrete time function, with the initial condition is:
The prediction of the historical time series of data at point k can now be deduced by:
3 The Optimized Nonlinear Grey Bernoulli Model ONGBM(1, 1)
Based on the basic principle of the information overlapping of the Grey system, Wang et al. (2009) [37] proposed a formula determining the value of the exponential parameter by:
with γ(k) =
However, with the existence of modern programming languages, this formula yields the problem of inflexibility. Moreover, as further investigated in Section 7, this formula yield a result which is significantly less effective, comparing to the method of iterating on a pre-defined range of the exponential parameter of [−1, 1).
According to Wang et al. (2011) [38] , the precision of any Grey model is highly dependent on the value of P in the background value z (1) (k). Change et al. (2005) [7] produced an even stronger statement by stating that forecasting accuracy of a Grey prediction model can be improved by optimizing parameter P .
The background value z (1) (k) of any (derived) Grey model is considered to be an approximation of the integral region generated by the curve of the function x (1) and the abscissa axis within the interval [k − 1, k]. Mathematically, this value is expressed as:
Applying the mean value theorem, it can be derived that:
which are exactly equations (8) and (18). Mentioned beforehand, P = 1 2 is the usually chosen values for traditional Grey models. However, it is trivial that the accuracy of the model depends heavily on the parameter P , which means that fixing this parameter at 1 2 as usual will reduce massively the effectiveness when there are values fluctuating and deviating from the overall trend.
Zhang (1993) [45] has proven that the parameter P determining the background values and parameter a satisfy the following condition:
From this formula, we derive:
= lim a→0 e a ae a + e a + e a L'Hôpital rule on
which is a valid explanation for the setting of the value of P to 1 2 for traditional Grey prediction models. However, when the value of a is significant (the absolute value of a is large enough), customarily setting the coefficient of P to 1 2 is incorrect, reducing the efficiency of the model.
As a result, a new calculating formula was proposed by Tan (2000) [34] as:
with q being given empirically by [26] :
This means that, we have defined the value of P by:
Once again, this formula yields the problem of inflexibility. As 1-AGO sequence is a strictly increasing sequence,
which means that
which has already omitted the range of 0, 1 2 for consideration. Moreover, as applied later, upon applying this formula, the algorithm yields higher error than the empirically derived value of the same parameter. The optimal solution for this problem would also be iterating on a well-defined range of values of P ([0, 1]) with a small enough step size.
Within the previous session, the parameters P and n are not optimized simultaneously.
Considering the particular solution of equation (17), the initial condition is set to be x (0) (1), the first point of the original time series of data. Using the principle of new information prior choosing, the last item of the 1 -AGO series, x (1) (m), can be used to enhance the accuracy of the model. Dang (2004) [12] has confirmed this theory by testing on two different Grey models: GM(1, 1) and the Grey-Verhulst model. Replacing x (1) (m) for (0) (1) as the initial condition, the time response function for the whitening equation of the NGBM(1, 1) model can be rewritten as:
Furthermore, since the model parameters a and b are obtained with least square method, there will be high chances that the model will not pass through the point x (0) (n), As a result, the new initial condition x (1) (m) will require correction in order to maximize optimization. This means that the new initial term, instead, would be x (1) 
The correction parameter c would be chosen in order to minimize the following function ( proposed by Lu et al. 2016 [26] ):
Substituting the discrete time response function deduced from Theorem 1, we obtain:
Let
Then, f (c) will be rewritten as:
To find the minimum value of f (c), we set the first -order condition (FOC) to 0, which means that ∂f (c) ∂c = 0, or:
We deduce
which means that the discrete time response function of the NGBM(1, 1) model can be now rewritten as:
Based on the simultaneous optimization of three parameters: the exponential parameter, the initial condition and the background value, we propose a further optimized algorithm of the Optimized NGBM(1, 1) given by Wu et al. (2019) [41] . This algorithm is more practical and less time-consuming upon implementing in the currently available computer languages (R, Python or Matlab):
1. Given the original time series of the data
2. Indicate the range of values for P and n. According to the mean value theorem, P must belong to the range [0, 1], while n is usually considered within the range of [−1, 1). Considering currently available computational tools, the step size is usually set at 0.001 − 0.01. Remark: Empirically, for n / ∈ [−1, 1), the model has been proven to yield higher error on average, reducing the efficiency of the model. As a result, the values of n outside the interval is not within our interest.
3. For the pair of values of P and n, calculate the matrices B and Y , which yields parameters a and b.
4. Calculate RPE and ARPE for each a and b. Then, choose the parameters which yield the minimum ARPE.
Calculate the values of E(k) and A(k) to deduce the error correction first initial condition x (1) (m) + c. Then, computex (1) based on the proposed discrete time response function (44).
6. Fromx (1) , using the inverse AGO, computex (0) .
As mentioned by Akay and Atak (2007) [3] , the rolling mechanism is an efficient technique to increase forecasting accuracy of Grey prediction when there are exponential increasing/decreasing or chaotic data. This mechanism provides the ability to update input data by removing the oldest data points for each loop before applying the prediction algorithm. Its purpose is that, for each rolling step, the set of data points used for the upcoming forecasting is the most recent.
Wu et al. (2019) [41] is the first paper to ever mention the application of rolling mechanism to the NGBM(1, 1) to enhance accuracy. However, the proposed model did not reach the maximum efficiency due to the following reasons:
1. In the proposed algorithm, the optimization of the parameters P and n are determined by running a loop through a set of pre-defined value by the user. In this case, the author defined the set of values of P to be between 0 and 1, with the step of 0.01; while n is considered between −1 and 0.99, with the step of 0.01. This step size for both parameters is still large enough to miss the optimal value.
2. In Algorithm 1, the background value is set as
which is not the optimal approximation and not compatible with the mean value theorem.
3. As proposed, re-initiating the initial value to x (1) (m) will enhance accuracy of the algorithm. However, this can further be improved by correcting this initial value itself with the error correction parameter c proposed in the previous section.
4. The rolling mechanism implemented in Algorithm 1 of [41] is incorrect. The rolling mechanism requires the insertion of the predicted value into the new array, while removing the oldest data point. In the algorithm, however, the author inserted the actual value instead of the predicted value, which develops biases for this model.
Assume that x (0) is a non -negative input of time series of data with m entries as:
and the output sequence of predicted data would be:
We propose a rolling Optimized NGBM(1, 1) model with the following algorithm:
1. Set i = 1.
2. Set the training data to
Use this training data to predict the following data points: 5. Our output sequence of predicted data would now be:
In order to examine the precision of each model, it is necessary to evaluate the difference between the fitted values and the actual values. For each specific time stamp k, the relative percentage error (RPE) calculates the relative difference of the forecast and recorded value by:
where x (0) (k) is the actual value andx((0)(k) is the fitted value.
The overall precision of the model can be measured by the average relative percentage error (ARPE), which is determined by:
The classification of the model precision based on ARPE is described in Table 1 below [41] . Besides ARPE, the root-mean-square error (RMSE) or root-mean-square deviation (RMSD) is also used to evaluate the precision of the model:
Based on the concept of RMSE, the posterior error ratio (c) is established and is also to evaluate the efficiency of a model:
The precision classification with respect to different levels of posterior error ratio is described in the following table 2 [6] . Within the scope of this paper, RPE and ARPE will be used as the main tool of benchmarking and comparing efficiency between different models. This approach is the most popular and widely-used upon considering empirical results of any newly-proposed Grey model.
In this section, we will indicate some alternative machine learning algorithms based on neural network. For this paper, we used 2 Neural Network models: Artificial Neural Network (ANN) and Long Short-term Memory (LSTM). Historically, recurrent Neural Network (RNN) and Convolutional Neural Network (CNN) were first introduced to handle such heavy-load tasks, like Image Processing or Natural Language Processing, and thus they were too complex for typical time series problems. However, researchers have started to adopt these deep learning methods, and turned them into perfect fits for a time-varying process. There are 2 main properties that make NN a competitive candidate for time series. NNs have general nonlinear function mapping capability that can approximate any continuous function. This gives NNs a capability of handling any complex problems with limited data. Moreover, NN is a non-parametric data-driven model, hence it does not require (mostly) any constraints. This is the key point for solving nonlinear, non-pattern problems that may have unique characteristics that cannot be captured by parametric model [36] .
Artificial Neural Network (ANN) is a brain-inspired system, which simulates the way human learns. A single layer neural network is called a Perceptron, which can give a single output.
To use ANN, we have to decide an Activation function, which converts an input signal of a node in an ANN to an output signal. This output signal is used as input to the next layer in the stack [8] . There are 4 Activation functions that can be taken into consideration.
• Threshold Activation function (Binary Step function) A threshold-based binary activation function. For classification only (because of activated cases only 1 or 0).
• Sigmoid Activation function (Logistic function) A S-shaped curve or sigmoid, ranges between 0 and 1. For predicting the probability
• Hyperbolic Tangent function (tanh) A similarity to sigmoid, but better in performance, ranges between -1 and 1. Perfect fit for layer stacking.
• Rectified Linear Units (ReLU) The most common Activation function in both CNN and ANN, ranges from zero to infinity [0, ∞).
Based on the difference between the actual value and the predicted value, an error value called Cost Function is computed and sent back through the system. Between an input layer and an output layer is at least a hidden layer. The hidden layers perform computations on the weighted inputs and produce net input which is then applied with activation functions to produce the actual output. ANN uses a procedure called Back-propagation cyclically until the error value is kept at minimum [8] . 
Long Short-term Memory (LSTM) is a RNN-based algorithm with the capability of learning long-term dependencies, which was proposed by Sepp Hochreiter and Jürgen Schmidhuber in 1997. Vanilla RNN method suffers from short-term memory during feed-forward. This means that, if processing data which has long time piece, it will be hard for RNN to carry (or remember) data from earlier time steps to later ones. Thus, RNN may leave out important information from the beginning. It also suffers from vanishing gradient problem in back-propagation. If a gradient value is too smal, it does not contribute much in learning process.
LSTM was created as a solution to this problem. LSTM has feedback connection, which allows it to process the whole data sequence. A general LSTM consists of multiple cells, with forget gates, an input and an output. They have the ability to choose which information to carry on or leave behind. [32, 15, 30] . [16] . Since then, the economic status has revived significantly, with the annual growth rate of 6-8% [17] .
As mentioned above, we used GDP data between 2005 and 2018, as from 2005 and after, Vietnam's GDP has some certain growth in economic status and annual growth rate. Moreover, Vietnam is on the verge of development and integration into the world, using this range of data will enhance the forecasting accuracy, as it can avoid the inaccurate that those previous world events bring.
In this section, we will regard the most important macroeconomics variable as our testing data: the annual Gross Domestic Product (GDP) of Vietnam, between 2005 and 2018. These data can be retrieved from World Bank Open Data (https://data.worldbank.org/)
In order to demonstrate the accuracy and effectiveness of the newly proposed Rolling Optimized NGBM(1, 1) (abberivated as RONGBM(1, 1)), we will proceed to compare this model with the existing Grey models, including:
• The traditional Grey model GM(1, 1)
• The traditional Nonlinear Grey Bernoulli Grey Model NGBM(1, 1)
• The Optimized Nonlinear Grey Bernoulli Model ONGBM (1, 1) . This model has adapted all possible optimization mentioned previously, including optimization of the background value, initial condition and exponential parameter.
The data from 2004 to 2013 will be used for training, and data from 2014 to 2018 will be used for forecasting and testing models. 
From Table 3 and Table 4 , it can easily be observed that the ARPE generated by the newly proposed Rolling Optimized NGBM(1, 1) model is significantly smaller than that of comparing models. This is due to the fact that the initial condition x (1) (m) + c and the original training data have been updated, which yields closer result to the actual data.
This empirical result also explains the conclusion driven in the previous sections that the estimating formula for the exponential parameter n and the background value P is not optimal, namely (25) and (31) . Upon using the equations above, the values of P and n would be: [31, 43, 22] . Researchers ared scared that this could be a new epidemic similar to the SARS-CoV and MERS-CoV epidemic in 2003 and 2012 separately that terrified worldwide.
The total infected cases were retrieved free of charge from Worldometer (https://www.worldometers.info/ coronavirus/coronavirus-cases/). We decided to use 2020-01-28 as the starting data point, as it would be more precise in terms of medical data collecting and analysis. Apart from comparing the Grey models, alternative modern machine learning models are also taken into account to compare accuracy and efficiency with Grey models.
In this section, the Rolling Optimized NGBM(1, 1) will not be taken into account considering the ARPE. As the training data sample is relatively small, this model will only be used to predict from the available data, based on the significant efficiency it has proven comparing to the Optimized NGBM in the previous example.
Considering the press release given be Reuters [39] stating that the peak of the COVID-19 outbreak will happen at midor late-February, we will proceed to predict the total infected cases everyday globally, starting from 2020-02-09 to 2020-02-18. 
Considering 5 different models predicting within the time range 2020-01-28 to 2020-02-08, the Optimized NGBM(1, 1) and even the traditional NGBM(1, 1) perform significantly better than two modern machine learning model (ANN and LSTM). This is due to the fact illustrated previously that NGBM(1, 1) yields relatively high accuracy upon handling short one -dimensional time series data.
Within the next 10 days, from 2000-02-09, similar to the example of predicting Vietnam's GDP 2013 − 2018, the predictions of the rolling optimized NGBM(1, 1) are lower than that of ONGBM(1, 1). This is due to the update of new data, reducing the increasing rate. However, both of the two predictions confirm that the peak of this COVID-19 outbreak will still not be happening within the next 10 days. From plot 9, it can be easily derived that the new infected cases by day is increasing, and there has been no signs of stopping. Figure 6 : Total COVID-19 infected cases prediction by ANN and LSTM (from 2020-01-28 to 2020-02-08) 
In this paper, a new rolling optimized NGBM(1, 1) model has been proposed. There are two following major advantages of this model in comparison with the existing ones:
• Within the proposed model, all of the parameters are optimized simultaneously, including the exponential parameter, the background value and the initial condition.
• The first application of the rolling mechanism upon the combination of all currently available methods of parameter optimization.
The efficiency of this newly proposed model has also been tested with the other Grey models (GM(1, 1), NGBM(1, 1) and ONGBM(1, 1)), which showed significant improvement. Later, while this model is applied to predict the total COVID-19 infected cases, it has also been put into comparison with the latest and most effective machine learning models, including Artificial Neural Network (ANN) and Long Short-Term Memory (LSTM). Within the context of extremely limit amount of data and time constrain, RONGBM(1, 1) outperforms both of the two other models, which additionally confirms the effectiveness of Grey models upon dealing with small amount of data.
From the results obtained by applying RONGBM(1, 1) on available data until the final day of retrieval (2020-02-10), there are now signs that the outbreak of COVID-19 will reach its peak within the next 10 days as expected. The number of newly infected patients reported everyday will even be higher than that of the present, increasing day by day. However, as Grey models are only applied on one-dimensional time series data, it is impossible to detect or predict changes due to any other related or outside effects (medical advancement, natural conditions, etc.). As a result, with enough amount of background knowledge into the field, applying a multivariate Grey model would yield more accurate and medically meaningful results.
Considering the proposed model RONGBM(1, 1) itself, there are two major problems that should be discussed and investigated in the future:
• Training data size: As the training data size has a significant impact on any Grey model [40] , the amount of data used in the COVID-19 infected case estimation has been manually chosen. In the future, for Grey models to reach maximum optimality, there should be a formal, mathematically derived procedure to define the optimal size of training data. • Parameter optimization: Within the proposed model, in order to obtain the optimization of the two parameters P and n, two loops have been run through. Due to the resource constrain of the current computational tools, the step size of these loops are still large enough to omit actual optimal values, leading to the model not perfectly optimized. Moreover, using the proposed formula, as demonstrated within the practical section, does not yield the optimal error. As a result, in the future, there should be a mathematically exact way to derive the parameter values optimally.