A regime switch analysis on Covid-19 in Romania

In this paper we propose a three stages analysis of the evolution of Covid19 in Romania. There are two main issues when it comes to pandemic prediction. The first one is the fact that the numbers reported of infected and recovered are unreliable, however the number of deaths is more accurate. The second issue is that there were many factors which affected the evolution of the pandemic. In this paper we propose an analysis in three stages. The first stage is based on the classical SIR model which we do using a neural network. This provides a first set of daily parameters. In the second stage we propose a refinement of the SIR model in which we separate the deceased into a distinct category. By using the first estimate and a grid search, we give a daily estimation of the parameters. The third stage is used to define a notion of turning points (local extremes) for the parameters. We call a regime the time between these points. We outline a general way based on time varying parameters of SIRD to make predictions.


A regime switch analysis on Covid-19 in Romania
Marian Petrica 1,2 , Radu D. Stochitoiu 3 , Marius Leordeanu 4,5 & Ionel Popescu 1,5* In this paper we propose a three stages analysis of the evolution of Covid19 in Romania. There are two main issues when it comes to pandemic prediction. The first one is the fact that the numbers reported of infected and recovered are unreliable, however the number of deaths is more accurate. The second issue is that there were many factors which affected the evolution of the pandemic. In this paper we propose an analysis in three stages. The first stage is based on the classical SIR model which we do using a neural network. This provides a first set of daily parameters. In the second stage we propose a refinement of the SIR model in which we separate the deceased into a distinct category. By using the first estimate and a grid search, we give a daily estimation of the parameters. The third stage is used to define a notion of turning points (local extremes) for the parameters. We call a regime the time between these points. We outline a general way based on time varying parameters of SIRD to make predictions.
In 2019, Covid19, a virus from the coronavirus family appeared and spread around the world very quickly. This changed dramatically our world as we knew it.
On 31st of December 2019, the first cases of infection with an unknown virus causing symptoms similar to those of pneumonia were reported in China, to the World Health Organization.
Within less than 3 months COVID-19 outbreak has become a global pandemic, spreading across almost all countries all over the world.
The fast-evolving spread of the new coronavirus, which has been officially declared a pandemic, is represented below. The charts in Fig. 1 show the countries where there have been reported at least 1000 cases of COVID-19 infection, at the mentioned date: While at the beginning of February 2020 the virus was still affecting mainly China, it has started to spread rapidly to other countries, and by the end of March 2020, the outbreak was present on all continents, affecting most of the countries in the world, which led the World Health Organization to officially name it a pandemic.
The main ideas of this paper. A basic tool in analyzing the spread of the virus is the mathematical modeling. There is a growing body of mathematical models used at the moment as for a small sample by no means exhaustive see [1][2][3][4][5][6][7][8] .
One of the characteristics of this pandemic is that there were many changes in the evolution. On one hand, the political decisions changed the course of the spread at the beginning with various measures, like quarantine, isolation, work from home and so on. Later on, other measures like relaxation, mask mandates, summer versus winter times, school openings and closing, election times, vaccination campaign, new variants of the virus and the travel ban lifting led to many changes in the status of the pandemic. Any reasonable parametric model of the epidemics has a major difficulty, namely the assumption that the parameters are constant over time is unrealistic. However, what one can still do is to use the models on short periods of time and then reassemble the local behavior to get a more general picture. This is our guiding principle in this paper.
Our approach to modelling the Covid19 evolution in Romania is in three stages. The starting point is the standard SIR model initiated in 9 and later investigated in depth in [10][11][12] . The basic SIR model uses the two basic parameters, the infection rate β and the recovery rate γ . They are assumed constant over time, however as we pointed out, the parameters of the model vary over time. However, we exploit this model on short periods of time where the assumption of constancy of parameters is still reasonable and thus we get daily estimates of the parameters based on 14 days of data. We should point out that 14 days seems to be a natural choice in the analysis for many reasons. For instance, the relative average period of recovery is 14 days. On the other hand, the infection takes some time to fully manifest. Also the global impact of a new variant of the virus takes a number of weeks till it is observed at large scale. Particularly to Romania, the health units have the obligation to report the cases involving Covid19 with an accepted delay of 14 days.
The mathematical basis for our first stage of the estimation is Proposition 1 which states that given a SIR model with constant parameters and data for two distinct days, we can completely determine the parameters. We exploit this results in combination with a neural network to do the estimation inspired by [13][14][15] . More details are outlined in "Stage I-SIR model and the first round of daily estimates" section. We only point out that the construction of this neural network is driven by the SIR model alone. We generate data and then train a neural network to learn the parameters. Using then the data, namely the number of infected and recovered, we estimate the main two parameters β and γ of the SIR model for short periods of time. This neural network construction could potentially be used in a more general framework of dynamical systems.
The second stage is driven by the idea that the number of deaths is more accurate and more reliable. Thus we propose a change of the SIR model to account for the dead as a separate category and create a differential equation associated to deaths. The idea is that the infected people are evolving into either recovered, as in the standard SIR model, or die. The parameters are now, β , the infection rate, γ 1 , the recovery rate and γ 2 the death rate, which models the rate at which the infected pass away.
The basic idea here is to take the outcome of the SIR estimates of the parameters β and γ as a first round of approximations and then proceed to a grid search of β, γ 1 , γ 2 around the suggested values from the first stage so that the model matches the observed number of deaths. We do this again, daily, utilizing the previous 14 days to have a more realistic estimates of the parameters.
The third stage is the definition of a regime. We look at the parameters β, γ 1 , γ 2 and identify the turning points (maximum/minimum) of the parameters. We pay more attention to γ 2 as predicting the number of deaths is presumably more important for the preparation of the medical units involved with fighting the pandemic. Though we payed more attention to γ 2 , the other parameters have extreme points approximately around the same values.
Having completed the three stages, we can actually use the analysis for predictions. We outline this by using a time varying SIRD model and natural estimates of the parameters using regression lines constructed in terms of 7 previous days. Based on this we show how one can make predictions. It turns out that for the prediction of deaths, this works well with two weeks of prediction and still reasonably well for three weeks forward. The  www.nature.com/scientificreports/ predictions cease in the proximity of the turning point justifying again the nomenclature of turning points. We insist on the methodology of the approach and complement this with numerical calculations and it can be extended to evolution equations for other diseases. The organization of the paper is as follows. In "The data" section we show the anomalies in data and then how we cleaned and adjusted it.
The main method is outlined in "Methodology" section. This is composed of "Stage I-SIR model and the first round of daily estimates" where we introduce the SIR model. We provide here the main mathematical result, namely Proposition 1 whose proof we postponed in the Appendix. Next, in "The neural network" and "Daily estimates of the parameters" sections we present the construction of the neural network and the first estimates of the parameters, that we get using the neural network. We continue then with "Stage II-SIRD model and the parameter estimates" section where we introduce the SIRD model and in "Parameter estimates" section we show the numerical scheme for the estimations. Furthermore in "Stage III-regimes and regime switch" section we provide the definitions of turning point and the regime, which is completed by the introduction of time varying SIRD model.
In "Predictions" section we provide an approach for the predictions using the estimates already done and we illustrate this with the case of prediction of deaths.
The last part, "Conclusions" section, is for concluding remarks. Finally, the Appendix provides the proof of Proposition 1.

The data
We import the data from https:// datel azi. ro which keeps a record of all the data during the pandemic in Romania starting with 17th of March 2020. We limit the data till 1st of February 2022.
One of the issues is that a first look at the raw data reveals an extreme spike in the number of recovered. This is due to the fact that the definition of recovered patient changed in October 2020. This added in one single day a very large number of recovered, approximately 44,000 cases, which is the same to the cumulative number of cases till that day. We redistributed these extra cases proportionally to the previous days.
The second correction is due to the fact that there are periods of zero reported numbers. This can range from a few days to almost three weeks. In addition, the reports during the weekends differ substantially from the reports during the weekdays. Also, by law, the medical centers have a flexibility of reporting the data with a delay of up to two weeks. Therefore, in order to alleviate these irregularities we use a moving average of two weeks. The results are presented in Fig. 2 with the data before and after the cleaning.

Methodology
In this section we outline the principles which guided us in this paper.
There are three main stages. The first stage consists in fitting a neural network on a SIR model to get a first round of parameters. Due to the uncertainties in the number of infected and recovered, we can not fully trust these results.
The second stage is to fit the parameters based on the more reliable data, namely the number of deceased people. To achieve this we propose a SIRD model which combined with the grid search provides a final estimate of the parameters.
The third and final stage introduces the definitions of turning point and regime, notions that we use in order to make predictions.
Stage I-SIR model and the first round of daily estimates. The SIR model. The first attempts of developing a mathematical model of the infectious diseases spreading were made at the beginning of the twentieth century. One of the most important models that can describe infectious diseases is the SIR model. The first ones that developed SIR epidemic models were Bernoulli, Ross, Kermack-McKendrick and Macdonald.
The SIR model is a mathematical model that can be used in epidemiology in order to analyze, at a given time for a specific population, the interactions and dependencies between the number of individuals who are susceptible to get an infectious disease, the number of people who are currently infected and those who have already been recovered or have died as cause of the infection. This model can be used to describe diseases that can be contracted just one time, meaning that a susceptible individual gets a disease by contracting an infectious agent, which is afterwards removed (death or recovery).
It is assumed that an individual can be in either one of the following three states: susceptible ( S ), infected ( Ī ) and removed/recovered ( R ). This can be represented in the following mathematical schema:

Susceptible
Infected Removed β γ where: • β = infection rate • γ = removed rate. www.nature.com/scientificreports/ We consider N as the total population in the affected area. We assume N to be fixed, with no births or deaths by other causes, for a given period of n days. Therefore, N is the sum of the three categories previously defined: the number of susceptible people, the ones infected, and the ones removed: Therefore, we analyze the following SIR model: at time t, we consider S (t) as the number of susceptible individuals, Ī (t) as the number of infected individuals, and R (t) as the number of removed/recovered individuals. The equations of the SIR model are the following: where: • dS dt is the rate of change of the number of individuals susceptible to the infection over time; N =S +Ī +R.
(1) The rows (from top to bottom) present the number of recovered, infected and dead individuals. The left column represents the raw data and the right column represents the adjusted data as we presented. As a detail, notice the scale and the spike in the first picture which is adjusted as we pointed out. The data is scaled by 10,000,000. is the rate of change of the number of individuals infected over time; • dR dt is the rate of change of the number of individuals recovered over time.
Because there is no canonical choice of N, we will transform the system (1) by dividing it by N and considering S(t) =S(t)/N , I(t) =Ī/N and R (t) =R(t)/N . It is customary to choose N = 10 6 for convenience but this is just an arbitrary choice. For instance, analysis on smaller communities or cities involves less than 10 6 , however 10 6 is a common choice because countries number their populations in multiples of 10 6 . With these notations we translate (1) into where β =β/N and γ is the same as in (1).
Notice that now we actually have that S(t) + I(t) +R(t) = S 0 + I 0 +R 0 = 1 for all t ≥ 0 . Since we are interested in the reverse problem, namely determining the parameters β, γ from the observations, we put this as a formal mathematical result as follows.
Proposition 1 Referring to the system (2), if we know I 0 , S 0 and the values I(t 1 ), S(t 1 ) for some t 1 > 0, these determine uniquely the parameters β and γ of the system.
It is important to remark that one of the main assumption is that the parameters β, γ do not change in time.
The neural network. Our next goal is to get a rough round of estimates on the parameters β, γ of the SIR model. For this we add the deaths to the recovered. Taking into consideration the recommendations/restrictions that have been applied by the authorities, in almost all countries (school closure, the ban of public events, social distancing recommendation/constraint, self-isolation if experiencing symptoms, quarantine for people tested positive), we presume that these parameters are not constant over time. We estimate the parameters based on two weeks of (cleaned) data.
To train the neural network we first use the SIR model to simulate data. We build a dataset based on the following procedure. Proposition 1 guarantees that if we know any values I 0 , I(t), R 0 , R(t) we can uniquely determine the parameters β, γ . We use the range of t = 1, 2, . . . , 14 to make the model more robust. At the same time the two weeks period is also consistent with the average recovery time of an infected patient of Covid19 and also corresponds to the lawful period of reporting of data.
We pick a sub sample (of size 80%) from the rows of and set which are the columns of corresponding to t, I 0 , I(t), R 0 , R(t) and the output data is exactly the pair The neural network we used is of the following form Daily estimates of the parameters. To estimate the parameters from the real data, we proceed as follows. For each day k = 0, 1, 2, . . . , T − 14 (T is the data range) we use the neural network to predict the parameters β k,t , γ k,t YTrain = �(β, γ ). www.nature.com/scientificreports/ ( t = 1, 2, . . . , 14 ) based on the real data I 0 = I real (k) , I(t) = I real (k + t) , R 0 = R real (k) , R(t) = R real (k + t) . Thus for each day k we determine 14 estimates of the parameters which are plotted in Fig. 3. Here by I real and R real we refer to the cleaned real data. In Fig. 4 we take the average of the parameters. For each day k we plot β k = 1   www.nature.com/scientificreports/ We should also comment on the fact that the data that is available shows the number of individuals that have been tested positive, but it is very likely that the real number of people infected is in fact much higher, as there are also asymptomatic individuals, people that are not being tested although they present the specific symptoms, so they are not part of the official reports.
Thus the above predictions for the parameters constitutes a good starting point for an optimization procedure we describe now. In order to reduce the effects of the above deficiencies, we consider another model which accounts for the number of deceased as a separate compartment.
Stage II-SIRD model and the parameter estimates. The SIRD model. In the sequel we propose a model which refines the SIR model. The guiding line is that the number of deaths is the most reliable number we can account for, as the number of infected and recovered people could be largely unaccounted.
We have now four variables changing with time. These are S(t), I(t), R(t) and D(t) where R(t) is the proportion of recovered and alive people while the D(t) is the proportion of deceased people. We set the interaction as follows where: • γ 1 = recovery rate • γ 2 = mortality rate Notice that in this setup the removed population, from classical SIR model, R , bifurcates into recovered ones, accounted by R and the deceased ones accounted by D. We can observe that R is the sum of the two factors R + D . In this way we separate the dead people from the recovered ones which are mixed up in the classical SIR model and we are going to manipulate these equations and reduce the computations to a single equation involving only one of these quantities, the most reliable one, namely D(t). To do this we will write all the other quantities as functions of D as follows: The easiest to deal with is R because from the last two equations we get from which we deduce that R(t) = γ 1 γ 2 (D(t) − D 0 ) + R 0 . Now, we deal with the function u from S(t) = u(D(t)) . Dividing the first and the last we get which can be integrated and gives S in terms of D as On the other hand this allows us to solve for I = v(D) . First we notice that from which we deduce that therefore we obtain that (as functions of D) This last implication works in the case the parameters β, γ 1 , γ 2 are all assumed constant in time. However, if they vary with time, then, the equation is a little bit different, the main equation becomes now www.nature.com/scientificreports/ At this moment we can use the data on the death cases to estimate the parameters involved. As we pointed out already, the proportion of infected (or recovered) is grossly underestimated since there are probably more infected people than the reported cases tested.
From the technical standpoint, Eq. (5) is not easy to handle and we will use Eq. (4) instead together with the implicit assumption that the parameters are constant for short periods of time. More precisely, in our approach we take the time interval on which we assume the parameters constant to be two weeks, which is in accordance with the lawful time of reporting and also with the dynamic of the time to recovery. In other words, we fit the number of deceased on pieces of two weeks where we assume that the parameters do not change.
Parameter estimates. The estimation of the parameters is done using a grid search based on the values already found with the SIR model estimates. The main reason is that we now use the previously found set of parameters as the starting point of the grid search. We do this dynamically, starting with any given day k and use the next 14 days forward to search for the set of parameters I k , R k , β k , γ 1,k , γ 2,k to find the ones which best predict the number of deaths.
We detail here the main steps.
1. For each day k we take the next 14 days of data.  Stage III-regimes and regime switch. We define first a turning point or a regime switch as a time where the parameters attain a local extreme (maximum or minimum). A regime is a period between two consecutive turning points.
To model the dynamic inside a regime, we use the linear regression in order to make prediction for the future periods.
Given times T = {t 1 , t 2 , . . . , t k } and positions X = {x 1 , . . . , x k } we define the regression line obtained using the least square optimization. Precisely, we set: We talked about the existence of different regimes in the spreading of Covid19 because of the measures that have been taken, which had a significant impact on the evolution of the infection rate. We consider now a regime starting at time p 1 , with a set of parameters, ending at time p 3 and having the intermediary time p 2 . We adapt the SIRD model as follows where σ β , σ γ 1 , σ σ 2 represent the regression functions based on a number of consecutive days and the parameters β, γ 1 , γ 2 evaluated at these days obtained from the optimization scheme (6). More precisely, we take for a fixed day k ≥ 7 , the previous 7 days and we use the regression line based on these days and the values of the estimated parameters. The effect of using the regression line has a smoothing effect on the parameters. As we plainly see, the parameters are not constant in time, not even for relatively short periods of time (say for instance 30 days). This new model is robust to the fluctuation of the parameters.
Below, in Figs. 7, 8, 9 we have the evolution of the estimates obtained in (6) together with the regression lines constructed in terms of 7 days. It is important to remark that the predictions lose their power around the turning points. In each of the cases we have a good indication of the accuracy with which the regression predicts the behavior of the parameters. With vertical dotted line we mark the points where we have a local extreme (minimum or maximum) in the evolution of γ 2 . The turning days are those days k for which the value of γ 2 (k) is an extreme value for the period [k − 7, k + 7] days centered at that given k day.

Predictions
Based on the method outlined above we predict the proportion of deaths into the future. The idea is to use now the system (8). More precisely, for any given day k we use the following recipe: 1. Use 7 previous days of the parameters produced by the estimation procedure from (6) to get the regression functions σ β , σ γ 1 , σ γ 2 . 2. Based on these data, we solve the system (8) for the next 14 or 21 days to see the fit with the real data of deceased.
The results are illustrated below in Fig. 10. Figure 6. This is the plot of the parameter γ 2 by day.

Conclusions
One of the characteristic of the Covid19 pandemic is that there were continuous changes of the conditions due to many factors. For instance, we can mention the quarantine periods, relaxation, mask enforcement, school openings and closings, vaccination campaign, major events, new variants of the virus, time of the year. If we take any parametric model, it is natural to assume that the parameters of the model change over time.
On the other hand, the reported numbers of infected and recovered during the pandemic is not reliable. Because of this reason we base our estimates on the reported number of deceased people which is, in our opinion, more realistic. However, discarding the data of infected and recovered entirely is not reasonable.
In this light, we proposed a three layers approach to the analysis of the pandemic of Covid19. The first layer consists in using the SIR model to estimate the first round of its parameters, β and γ assumed constant for periods of up to 14 days. Technically, we do this using a neural network trained on data simulated with the SIR model. Using this neural network and real data for any period of 14 days we generate a pair (β, γ ) which best describes the period from the point of view of the SIR model.
The second layer is based on a modification of the SIR model to account for the dead individuals. We call this SIRD and to our knowledge it is not used in this framework in the literature (though it is also recently used in cite). The SIRD model is depending on three parameters β, γ 1 , γ 2 . Based on the result of the previous layer, we use an optimization problem to fit the number of reported deaths. We do this procedure assuming the parameters are constant on pieces of 14 days. The result provides, for any given day, a prediction of the parameters of the model using the previous 14 days. www.nature.com/scientificreports/ The last layer consists in using the parameters to determine the turning points, the extreme points of the parameters (more precisely we do this for γ 2 ). We define a regime as the time period between two turning points. This perspective is consistent with our assumption that the parameters are not constant for long periods of time.
This last layer, combined with an adaptation of the SIRD model to account for the change in parameters, leads to a way to predict the evolution of the pandemic, at least for short periods of time.
We believe that this proposed methodology is a general one and can be extended to the analysis of spread of Covid19 in any country provided that we have relevant data. This methodology can also be extended to many other diseases which can be modeled by SIR and SIRD. www.nature.com/scientificreports/ Figure 9. The first plot shows the estimated parameter γ 2 together with the regression lines started at each day, using the last 7 days data and plotted for 7 days. In the second one, we plotted the regression line for 14 days, based on 7 days from the data. Typically the initial value of S 0 is close to 1 and I 0 is relatively small. In particular, if we assume that the epidemic ends somewhere then we definitely have I(t) = 0 and thus S(t) solves the equation In particular if we assume that I(t ∞ ) = 0 and S(t) converges as t → t ∞ , then we get in the limit that S(t ∞ ) solves (11). One consequence of this argument is that for all time 0 ≤ t ≤ t ∞ , we have that S(t) − ρ log(S(t)) ≤ α := I 0 + S 0 − ρ log(S 0 ).
Another important consequence of this model is that if we assume S 0 and I 0 fixed (obviously R 0 will also be determined) but, for a given time t = t 1 > 0 , knowing S(t 1 ) and I(t 1 ) (therefore R(t 1 ) as well), we can determine uniquely the parameters β and γ . Indeed this is clearly seen from (10) which gives On the other hand, from (10) in the first line of (2), then we obtain that I(t) + S(t) − ρ log(S(t)) constant in t.