prediction
Christos Papahristodoulou
May 24, 2020
In this document I apply the standard SIRD model , using Italian COVID19 data. Up to date, May 23,
there are 89 daily observations, starting from February 24. The model is adjusted and does not
cover all possible people to be infected, but only those who been tested and were found positive
until May 23. In Italy more than 3.3 million tests were curried so far and the accumulated number of
infected people is 228,658 (almost 7%). The infected/tested ratio over the recent days is very low, at
most 1% (about 650 new infected per day), compared to more than 30% two months ago. Today is
the 89th day and I will solve the problem to predict the key variables until the middle of June (t =
115). Assuming the current number of new infected is around 600 per day, the susceptible
popula-tion is therefore set equal to about 245,000 over the period under investigapopula-tion.
After I have retrieved the data, I use 79 out of 89 observations to estimate the key SIRD parameters.
Thereafter, using these estimates I solve the system of ODE equations and plot their paths together
with the current observations.
Retrieving the data
In[ ]:= URLDownload[ "https://raw.githubusercontent.com/pcm-dpc/COVID-19/master/dati-andamento-nazionale /dpc-covid19-ita-andamento-nazionale.csv", Directory[] <> "\\pd.read_csv(url)"]; In[ ]:= SemanticImport[ "https://raw.githubusercontent.com/pcm-dpc/COVID-19/master/dati-andamento-nazionale /dpc-covid19-ita-andamento-nazionale.csv"];
Since I will apply the standard SIRD model, I will use the following columns (all are accumulated
statistics): infected (7), recovery (10), deaths (11), and total (12).
Estimates on parameters
In order to estimate the key SIRD (Susceptible, Infected, Recovery, Death) parameters, β, γ, δ, I will
use 79/89 of the observations (until May 13).
In[ ]:= infstart = Import[ "https://raw.githubusercontent.com/pcm-dpc/COVID-19/master/dati-andamento-nazionale/dpc-covid19-ita-andamento-nazionale.csv", {"Data", 2 ;; 80, 7}]; recstart = Import[ "https://raw.githubusercontent.com/pcm-dpc/COVID-19/master/dati-andamento-nazionale/dpc-covid19-ita-andamento-nazionale.csv", {"Data", 2 ;; 80, 10}]; deathstart = Import[ "https://raw.githubusercontent.com/pcm-dpc/COVID-19/master/dati-andamento-nazionale/dpc-covid19-ita-andamento-nazionale.csv", {"Data", 2 ;; 80, 11}]; totstart = Import[ "https://raw.githubusercontent.com/pcm-dpc/COVID-19/master/dati-andamento-nazionale/dpc-covid19-ita-andamento-nazionale.csv", {"Data", 2 ;; 80, 12}] // Reverse;
I need to normalize the data by dividing their values by 221,216, i.e. the total infected number until
May 13. The model is:
eqS = s '[t] ⩵ - β s[t] × m[t];
eqI = m '[t] ⩵ β s[t] × m[t] - γ m[t] - δ m[t]; eqR = r '[t] ⩵ γ m[t];
eqD = q '[t] ⩵ δ m[t];
where, s = susceptible, m = infected, r = recovery, q = deaths.
In[ ]:= time = Range[79];
In[ ]:= datastart = ThreadJoin{time}, {totstart} 221 216,
{infstart} 221 216, {recstart} 221 216, {deathstart} 221 216 // N;
In[ ]:= f[β_ ? NumericQ, γ_ ? NumericQ, δ_ ? NumericQ] :=
SumTotalyi[[i, All]] - Map[pfun[β, γ, δ][[i]], xi] ^ 2, {i, 1, 4} // Quiet data = datastart;(*take only 79 days*)
xi = data[[1, All]];(*independent variable*) yi = data[[2 ;; 5, All]];
(*four dependent variables*)pfun = ParametricNDSolveValue[{s '[t] ⩵ - β s[t] × m[t], m '[t] ⩵ β s[t] × m[t] - γ m[t] - δ m[t], r '[t] ⩵ γ m[t], q '[t] ⩵ δ m[t], s[xi[[1]]] ⩵
yi[[1, 1]], m[xi[[1]]] ⩵ yi[[2, 1]], r[xi[[1]]] ⩵ yi[[3, 1]], q[xi[[1]]] ⩵ yi[[4, 1]]}, {s, m, r, q}, {t, 0, 80}, {β, γ, δ}];(*the time period is set to 80 days*)
In[ ]:= fit = NMinimize[f[β, γ, δ], {β, γ, δ}];
params = fit // Last
Out[]= {β → 0.200637, γ → 0.0222591, δ → 0.00896591}
Solve the SIRD system using the estimated parameters
Using these estimates I solve the system of ODE equations and plot their paths together with the
current observations.
In[ ]:= Clear[t];
In[ ]:= Clear[s, m, r, q];
In[ ]:= eqS = s '[t] ⩵ - β s[t] × m[t];
eqI = m '[t] ⩵ β s[t] × m[t] - γ m[t] - δ m[t]; eqR = r '[t] ⩵ γ m[t];
eqD = q '[t] ⩵ δ m[t];
In[ ]:= β = 0.200637; γ = 0.0222591; δ = 0.00896591;
In[ ]:= solu = NDSolve[{eqS, eqI, eqR, eqD, s[0] ⩵ 1, m[0] ⩵ 0.00096, r[0] ⩵ 0, q[0] ⩵ 0},
{s, r, m, q}, {t, 115}];(*initially there were about 220 persons*)
In[ ]:= soluS = First[s /. solu];
soluI = First[m /. solu]; soluR = First[r /. solu]; soluD = First[q /. solu];
In[ ]:= Plot[{soluS[t], soluI[t], soluR[t], soluD[t]}, {t, 0, 115},
PlotRange → {0, 1}, PlotStyle → {Blue, Red, Green, Magenta}, PlotLegends → {"Susceptible", "Infected", "Recovered", "Dead"}]
Out[]= 0 20 40 60 80 100 0.2 0.4 0.6 0.8 1.0 Susceptible Infected Recovered Dead
Plot the data and the solution of the ODE
We can now plot all data as well as the ODE.
In[ ]:= infall = Import[
"https://raw.githubusercontent.com/pcm-dpc/COVID-19/master/dati-andamento-nazionale /dpc-covid19-ita-andamento-nazionale.csv", {"Data", 2 ;; 90, 7}]; recall = Import[ "https://raw.githubusercontent.com/pcm-dpc/COVID-19/master/dati-andamento-nazionale /dpc-covid19-ita-andamento-nazionale.csv", {"Data", 2 ;; 90, 10}]; deathall = Import[ "https://raw.githubusercontent.com/pcm-dpc/COVID-19/master/dati-andamento-nazionale /dpc-covid19-ita-andamento-nazionale.csv", {"Data", 2 ;; 90, 11}]; totall = Import["https://raw.githubusercontent.com/pcm-dpc/COVID-19/master/dati-andamento-nazionale/dpc-covid19-ita-andamento-nazionale.csv", {"Data", 2 ;; 90, 12}] // Reverse;
We normalize again these data by dividing their values by 245,000, i.e. the total expected number of
infected, by the middle of June.
In[ ]:= beta = ListPlot
Tooltiptotall 245 000, infall 224 500, recall 224 500, deathall 224 500, PlotStyle → {Blue, Red, Green, Magenta},
PlotLegends → {"Susc", "Inf", "Rec", "Dead"};
In[ ]:= Show[Plot[{soluS[t], soluI[t], soluR[t], soluD[t]}, {t, 1, 115}, PlotLabel →
Style["Italia:COVID; Feb 24 - May 23 ", Black, Bold, 14], PlotRange → {0, 1}, PlotStyle → {Darker[Blue], Darker[Red], Darker[Green], Darker[Magenta]}, GridLines → Automatic, Frame → True, PlotRange → All, ImageSize → 500], beta]
Out[]= 0 20 40 60 80 100 0.0 0.2 0.4 0.6 0.8
1.0 Italia:COVID; Feb 24 - May 23
Susc Inf Rec Dead
In[ ]:= Manipulate[
Show[Plot[a {soluS[t], soluI[t], soluR[t], soluD[t]}, {t, 1, T}, PlotRange → Automatic], beta, PlotStyle → {Darker[Blue], Darker[Red], Darker[Green], Darker[Magenta]},
ImageSize → 450, Evaluated → True, PlotPoints → 1000, GridLines → Automatic, Frame → True, PlotTheme → "Detailed"], {{T, 90}, 50, 115, Appearance → "Labeled"}, {{a, 0.94}, 0.9, 1.0, Appearance → "Labeled"}, ControlPlacement → Top,
Alignment → Right, SaveDefinitions → True]
Out[]=
T 55.4
In[ ]:= Manipulate
ColumnTextStyle"Susceptible =" Text[Style[N[228 658 soluS[t] /. solu ]]], 16, Blue, TextStyle"Infected =" Text[Style[N[228 658 soluI[t] /. solu]]], 16, Red,
Text[Style["Recovered =" Text[Style[N[228 658 soluR[t] /. solu]]], 16, Darker[Green]]], Text[Style["Dead =" Text[Style[N[228 658 soluD[t] /. solu]]], 16, Magenta]],
Text[Style["Ro =" Text[Style[N[(0.200637 / 0.022259) soluS[t] /. solu]]], 16, Black]], {t, 8, Text[Style[Row[{"Days"}], 16]]}, 2, 115, 1, LabelStyle → DirectiveBlack, Italic, 18,
Appearance -> "Labeled", ControlPlacement → Top,
ImageMargins → {50, 0}, Automatic, 0, SaveDefinitions → True
Out[]=
