infection. In fact this parameter m represents the mean value of the infected individuals by one individual per day. Using the observed statistics some methods for estimation are proposed and corresponding graphics are presented. Two models with and without immigration are compared. In this way we are able to give a prediction of the possible development of the mean value of the infected individuals.
Notice that both type processes with or without immigration have an exponential growth in the supercritical case m > 1 but in the critical case m = 1 and in the subcritical case m < 1 the asymptotic behaviour is essentially different. In the critical case the mean value of the process with immigration grows linearly while for the process without immigration the mean value is constant. In the subcritical case the mean value of the immigration process converges to a positive constant but for the process without immigration the mean value goes to zero. The estimation of the immigration mean a is a serious problem because this parameter cannot be estimated by the observed data and one needs additional information. Hence the processes with immigration need more careful investigation.
As it is given in the paper the proposed estimators can be applied also in the case when the processes are inhomogeneous in time. The behaviour of the estimators shows that the observed processes are able to change the criticality during the development of the epidemics. The estimated values for the main epidemical parameter m vary greater than 1, equal (or very closed) to 1 and even less than 1. Moreover it seems that the real epidemic process develops like a mixture of both type of the models with and without immigration. In this case four stages of epidemic development are available: exponential growth (m > 1), linear growth (m = 1 and an immigration component), non-increasing and almost stable population (due to m = 1 or m < 1 with an immigration component), convergence to zero (extinction of the epidemics due to m < 1 without immigration). It is obvious that the restriction of the immigration component is very important to the limitation of the epidemic process. It seems that in some countries exist some regions which can be considered as immigration sources for the other regions and in general for the whole country. The abroad immigration plays also an important role.
The paper continues the investigations started in [16] where some results for the model without immigration were presented only for Bulgaria, Italy and globally.
In the present paper as an illustration of the models with and without immigration the obtained results are presented for several countries: USA, Italy, France, Germany, Spain and Bulgaria. Additional information, reports and plots, related to this research for all countries all over the world can be found on the site http://ir-statistics.net/covid-19. The data used for the estimation of the parameters of the model are taken from European Centre for Disease Prevention and Control [21] , similar to the data provided by World Health Organization [20] . Since these databases are updated daily, the proposed here model is applied regularly on each new data set. Using these results one can compare the infection rate on different countries and regions on the basis of es-timated growth rate. For example, on Table 3 , the 10 countries with lowest and highest growth rate are shown. Even in the cases where the infection growth is less than 1, the confidence interval goes above 1, which states that there is a possibility for increasing of the infection growth in the future.
The theoretical model based on two type branching process is described in detail in Section 2. The p.g.f.'s and the mathematical expectations are obtained. Regardless of its simplicity the model has a great advantage using only the observed official data for the lab-confirmed cases. The two-type branching process assuming an additional immigration component is considered in Section 3. The estimation problems are presented in Section 4. Some conclusive remarks are given in Section 5.
Finally the estimation of the mean parameter of infection can be considered as a fast test to estimate the rate of Covid-19 epidemic in a country or a region. It can be used also as a first stage of a construction to some more complicated epidemiological models where it is not possible to estimate directly this parameter. Obviously the solution of this problem requires the collaboration of specialists in various fields as epidemiology, mathematics, medicine, microbiology, molecular biology and informatics among others.
2 Two-type branching process as model of Covid-19 population dynamics.
Assume that the epidemic process of infection begins with some finite number of immigrants and then the process of immigration is isolated under the quarantine.
To describe this situation we can consider a two type branching process {Z 1 (n), Z 2 (n)} where type T 1 are infected (but still healthy) individuals who don't know that they are Covid-19 infected and type T 2 of discovered with Covid-19 virus individuals (and this is the data we use). Every individual of type T 1 (infected) produces per day a random number of new individuals of type T 1 (infected) or only one individual of type T 2 (more precisely, in this case the individual type T 1 is transformed into an individual type T 2 ). Note that T 2 is a final type, i.e. the individuals of this type don't take part in the further evolution of the process because they are isolated under the quarantine.
Let
2 ) be the offspring vector of type T 1 . Then the offspring joint probability generating function (p.g.f.) of type T 1 can be defined as follows:
Obviously h 2 (s 1 , s 2 ) ≡ 1 because the type T 2 has (0, 0) offspring. Note that p 0 is the probability that type T 1 goes out of the reproduction process (the individual becomes healthy or goes out of the country, i.e. emigrates), p j is the probability to produce new j infected individuals of type T 1 and q is the probability that the individual type T 1 is confirmed ill (or dead). In other words, q = P{T 1 → T 2 }, i.e. with probability q an individual of type T 1 is transformed into an individual of type T 2 . Then from (1) we can obtain also that the marginal p.g.f. are
If we assume that Z 1 (0) > 0 and Z 2 (0) = 0 then for n = 1, 2, ...
2 (n; j)} are independent and identically distributed (iid) as (ξ
2 ). The recurrent formula (4) defines two-type branching process {(Z 1 (n), Z 2 (n)), n = 0, 1, 2, ...}. Notice that Z 1 (n) is the total number of individuals (type T 1 ) in the n-th day infected by the individuals of the (n − 1)-th day; Z 2 (n) is the total number of the registered Covid-19 individuals (type T 2 ) in the n-th day. The process starts with Z 1 (0) infected individuals, where Z 1 (0) can be an integervalued random variable with a p.g.f. h 0 (s) = Es Z0 = K k=1 p 0k s k , |s| ≤ 1, or Z 0 = N for some integer value, N = 1, 2, ... . The random variable ξ (1) 1 (n; j) is the number of individuals of type T 1 in the n-th day infected by the j-th individual of type T 1 from the (n − 1)-th day, j = 1, 2, ..., Z 1 (n − 1). Similarly the random variable ξ (1) 2 (n; j) is the number of the confirmed infected individuals (type T 2 ) in the n-th day transformed by the j-th infected individual type T 1 from the (n − 1)-th day, j = 1, 2, ..., Z 1 (n − 1).
Note that P{ξ
.., l; l = 0, 1, 2, ... In other words the probability q can be interpreted as a proportion of the confirmed individuals in the day n among all infected individuals in the day n − 1.
Let h 0 (s) = Es Z1(0) , F 1 (n; s) = E(s Z1(n) ), F 2 (n; s) = E(s Z2(n) ). Introduce the following p.g.f.
Then it is not difficult to check that for n = 0, 1, 2, ..., we are able to obtain the p.g.f. of the process: 
Notice that the asymptotic behaviour of the process depends essentially of parameter m. Especially, if m > 1 (supercritical case) then the mean value of the infected individuals M 1 (n) grows exponentially, in the critical case m = 1 it is a constant and for m < 1 (subcritical case) M 1 (n) → 0 as n → ∞.
We will use these results to present in the next section a more complicated model with immigration.
Note that we can observe only Z 2 (1), Z 2 (2), ..., Z 2 (n). What can be estimated with these observations?
Note first that EZ2(n+1) EZ2(n) = m. Hence we can consider (11) m n = Z2(n+1) Z2(n) , n = 1, 2, ...., as an estimator of the parameter m (similar to Lotka-Nagaev estimator for the classical BGW branching process). It is possible to use also the following Harris type estimator (12) m n = n+1 i=2 Z 2 (i)/ n j=1 Z 2 (j), n = 1, 2, ...., or Crump and Hove type estimators (13) m [12] for more details. Estimating m we are able to predict the mean value of the infected (non observed) individuals in the population. In the case when we assume that Z 1 (0) = 1 then M 1 (n) = EZ 1 (n) can be approximated respectively by m n n , or m n n , or m n n,N . In fact it means that we can obtain three types of estimators (14) M 1 (n) = m n n , M 1 (n) = m n n and M 1 (n) = m n n,N . In other words we could say that we have at least three prognostic lines. Therefore if we have the observations (Z 2 (1), Z 2 (2), ..., Z 2 (n)) over the first n days, we are able to predict the mean value of the infected individuals for the next k days by the relations: (15) M 1 (n+k) = m n+k n , M 1 (n+k) = m n+k n and M 1 (n+k) = m n+k n,N , k = 1, 2, ... We are able to estimate also the proportion α(n) of the registered infected individuals among the population in the n-th day. Then we can obtain the following three types of estimators:
All obtained estimators will be presented by the observed registered labconfirmed cases. The quality of the estimation, however, depends on the representativeness of the sample due to the specifics of the data collection in each country.
Remark 1. In fact our model (4) can be generalized as non-homogeneous in time. In this case
jp j (l) will be the mean value of the new infected individuals by one infected individual in the day l = 1, 2, .... Therefore instead (9) we obtain (16)
Notice that in this case EZ2(n+1) EZ2(n) = m(n). Hence we can use (11) − (13) to estimate M 1 (n) from (16). Therefore (17) M 1 (n) = Π n l=1 m l , M 1 (n) = Π n l=1 m l , n = 1, 2, .. 3. The two-type branching process with immigration as model of Covid-19 epidemic development.
The model considered in Section 2 assumes that the process of infection begins with some random number of infected immigrants and then the immigration process is bounded and it is not essential for the Covid-19 population dynamics. But in some cases the role of the immigration process cannot be negligible. That is why we will introduce random variables {I n }, where I n gives the number of infected immigrants in the n − th day which take part in the process of infection. We will assume first that {I n } are iid r.v. with a p.g.f. (18) g(s) = Es In = l k=0 g k s k , |s| ≤ 1. Then instead of (4) we will consider the following branching process with immigration
2 (n; j)} are independent and identically distributed as (ξ
2 ) with p.g.f. (1) − (3) and they are also independent of {I n }. We can assume that Y 1 (0) > 0 is some random variable independent of
2 (n; j)} and {I n }, and also the case Y 2 (0) = 0. Another possible assumption is Y 1 (0) = Y 2 (0) = 0, which means that in fact the process starts with the first real immigrants.
Interpretation: Y 1 (n) is the total number of individuals (type T 1 ) in the n-th day infected by the individuals of the (n − 1)-th day plus the new infected immigrants I n ; Y 2 (n) is the total number of the officially registered infected individuals (type T 2 ) in the n-th day.
Then it is not difficult to check that for n = 0, 1, 2, ..., we are able to obtain from (18) and (19) the p.g.f.'s of the process: G 1 (n; s) = E(s Y1(n) ) = g(s)G 1 (n − 1; h * (s)) = f (h * n (s))Π n−1 k=0 g(h * k (s)), G 2 (n; s) = E(s Y2(n) ) = G 1 (n − 1; h(s)) = g( h(s))G 1 (n − 2; h * ( h(s))) = f (h * n−2 ( h(s)))Π n−1 k=0 g(h * k ( h(s))), where the p.g.f. h * k (s) and h k (s) are obtained after k iterations of the p.g.f. h * (s) and h(s) as it is given in (8) and
Notice that if we assume that Y 1 (0) = 0 then f (s) ≡ 1 and (20) (18) we can introduce the immigration mean a = EI n = d ds g(s)| s=1 = l k=0 kg k . Then from (18)−(20) it is not difficult to obtain that for n = 1, 2, ...
where it is assumed that A 1 (0) = EY 1 (0) = 0 and the parameters m and q are well defined in Section 2 by (5).
Hence from (21) and (22) one has (23) A 1 (n) = a(m n − 1)/(m − 1), m = 1, and A 1 (n) = an, m = 1, (24) A 2 (n) = qa(m n−1 − 1)/(m − 1), m = 1, and A 2 (n) = qa(n − 1), m = 1 Therefore by (23) and (24) one obtains as n → ∞ A 1 (n) ∼ am n /(m − 1), m > 1; A 1 (n) = an, m = 1; A 1 (n) → a/(1 − m), m < 1, A 2 (n) ∼ qam n−1 /(m − 1), m > 1; A 2 (n) ∼ qan, m = 1; A 1 (n) → qa 1−m , m < 1. In the general case A 1 (0) = EY 1 (0) = M 0 > 0 and instead of (19) and (20) one has (25) A 1 (n) = M 0 m n + a n−1 k=0 m k , (26) A 2 (n) = q(M 0 m n−1 + a n−2 k=0 m k ). We would like to point out once again that we can observe only the statistics Y 2 (1), Y 2 (2), ..., Y 2 (n) and we have to use for estimation only these observations.
Notice first that for m ≥ 1 we obtain lim n→∞ EY2(n+1) EY2(n) = m. Hence for large enough n we can consider (27) m n = Y 2 (n + 1)/Y 2 (n) as an estimator of the parameter m (similar to Lotka-Nagaev estimator for the classical BGW branching process). It is possible to use also for m > 1 and large enough n the following Harris type estimator (28) m n = n+1 i=2 Y 2 (i)/ n j=1 Y 2 (j) or Crump and Hove type estimators
. See [12] for more details. Estimating m > 1 we are able to predict the mean value of the infected (non observed) individuals in the population. In the case when we assume that the process begins with the first immigrants then A 1 (n) = EY 1 (n) can be approximated using the estimators (27) − (29).
The problem is how to estimate the immigration mean a. First of all there is an special case when a = m. Then using (23) and (24) with the Harris estimator we have (30) A 1 (n) = m n ( m n n − 1)/( m n − 1), m > 1; A 1 (n) = m n n, m = 1. In general we have to use some additional information. For example, if we can observe {I k } then we can apply the estimator a * n = n −1 n k=1 I k . Hence (31) A 1 (n) = a * n ( m n n − 1)/( m n − 1), m > 1; A 1 (n) = a * n n, m = 1. One can proceed similarly for the other estimators m n n and m n n,N . Remark 2. Similarly as it is shown in Remark 1 from Section 2 the model (19) can be generalized in the case with non-homogeneous in time offspring distributions. In this case m(l) = d ds h * (l; s)| s=1 = Eξ 
Recall that both type processes with or without immigration have exponential growth in the supercritical case m > 1. In the critical case m = 1 and in the subcritical case m < 1 the asymptotic behaviour is essentially different. In the critical case the mean value of the process with immigration grows linearly while for the process without immigration the mean value is constant. In the subcritical case the mean value of the immigration process converges to a positive constant but for the process without immigration the mean value limit is equal to zero. The estimation of the immigration mean a cannot be estimated by the observed statistics and we need some additional information.
We would like to point out once again that the considered in Section 2 model is versatile but the application in each country is specific because it depends essentially on the official data from the country. The plots and tables below illustrate well some specific details for different countries as well as the common trend.
The data used for the estimation of the parameters of the model come from European Centre for Disease Prevention and Control [21] .
We will consider first the process without immigration. Note that the observed data is the number of the newly (daily) registered individuals denoted by Z 2 (n). The data about the new number of infected individuals (denoted by Z 1 (n)) is unobservable. The initial number m 0 = EZ 1 (0) is also unknown. Here n is the corresponding day from the beginning of the infection.
The estimation of the parameters of the defined model can be summarized in the following steps. Firstly, we will demonstrate the approach described above by the data of the reported laboratory-confirmed COVID-19 daily cases for USA provided by the European Centre for Disease Prevention and Control [21] (the data are retrieved on 27.04.2020). Table 1 represents the estimated model parameters for the last 5 days of the available data set. Every row in the table represents the Harris estimatem n , as well as it's 95 % confidence interval ( CI l -CI r ), the proportion of the registered infected individuals α(n) and the expected values of the non-confirmed cases M 1 ( or A 1 for the process with immigration) , based on n − k observations, i.e. Z 2 (1), ..., Z 2 (n − k); k = 0, ..., 4. Table 1 : Estimation of the model parameters
The following figures Figure 1 -6 represent the parameters of the process for the USA data. Figure 3 . After the initially large estimated values it stabilizes below 1.1, which is determined by the branching processes theory as a slightly supercritical process. This corresponds to the exponential growth shown above. The next results shows that the Harris type estimator has more stable behaviour than the Lotka-Nagaev type estimator. The last 5 days results for Italy, Germany, France, Spain and Bulgaria can be compared on Table 2 (the data are retrieved on 27.04.2020).
The value of the proportion of the registered infected individuals is considerably higher in USA and France than in Germany and Italy, while countries with a longer infection period observe relatively small values of the Harris estimator of the mean value of the number of the confirmed infected individuals by one infected individual. Even more, the lower boundary of the confidence interval for the Harris estimator falls beneath the value of 1.
Using the same data set one can compare the infection rate for different countries and regions on the basis of the estimated growth rate. For example, on Table 3 , the 10 countries with lowest and highest growth rate are shown. Even in the cases where the infection growth is less than 1, the upper bound of the confidence interval goes above 1, which states that there is a possibility that the infection growth will increase in the future.
In the countries, where the infection is growing, the growth rate is isoclinically above 1, even though the lower boundary of the confidence interval is less First of all the estimation of the mean value of reproduction m allows us to classify the contamination process as supercritical (m > 1), critical (m = 1) and subcritical (m < 1). Recall that both type processes with or without immigration have exponential growth in the supercritical case m > 1. In the critical case m = 1 and in the subcritical case m < 1 the asymptotic behaviour is essentially different. In the critical case the mean value of the process with immigration grows linearly while for the process without immigration the mean value is constant. In the subcritical case the mean value of the immigration process converges to a positive constant but for the process without immigration the mean value limit is equal to zero.
Finally the estimating of the mean parameter of infection can be considered as a first stage to construction of a more complicated epidemiological model. As an example, one can use a branching process with random migration considered in [17 − 19] or some other model of controlled branching processes (see [8] ). But for a general pandemic model the collaboration with specialists of epidemiology, mathematics, medicine, microbiology, molecular biology and informatics is absolutely necessary for the application of all information about Covid-19 phenomenon.
Remark 3. For more detailed investigation and simulation the following models can be applied in the considered situation: (i) h * (s) = q + p 0 + p 1 s + p 2 s 2 + ... + p k s k , where q = 1 − k j=0 p j and p j , j = 0, 1, ..., k, can be specially chosen for k = 2, 3, 4, 5, 6, 7, 8.
(ii) h * (s) = q + p 0 + ∞ k=1 (1 − p)p k s k = q + p 0 + (1 − p)ps)/(1 − ps), where q + p 0 = 1 − p. It is possible to consider also the restricted geometrical distribution up to some k = 2, 3, 4, 5, 6, 7, 8.
(iii) h * (s) = q + p 0 + ∞ k=1 e −λ λ k k! s k = q + p 0 + e −λ(1−s) − e −λ , where q + p 0 = e −λ . Similarly it is possible to consider also the restricted Poisson distribution up to some k = 2, 3, 4, 5, 6, 7, 8.
Note that the parameters of these distributions can be set in the manner that d ds h * (s)| s=1 is equal to m n , or m n , or m n,N . Then with this individual distributions it is possible to simulate the trajectories of the non-observed process of infection for further studies. Additional information, reports and plots, related to this research can be found on http://ir-statistics.net/covid-19. The presented results are updated every day following the new data which are provided day by day from European Centre for Disease Prevention and Control [21] .