Forecasting the number of hospitalized - for all T

Modules

Hospitalized

Hospitalized[DataT_,Ind_:0]:=Module[{},
par=0;
danIni=First[DataT][[{1,2}]]; danEnd=Last[DataT][[{1,2}]];
delU=15;  (* Data delay for deceased and hospitalized *)
axesstyle={Thickness[.004],Thickness[.004]};
delta=10;step=10;
ticksx=Table[{i,DataT[[i,{1,2}]]},{i,delta,Length[DataT],step}];
ticksdel=Table[{tabx[[i]],Data[[10i,{1,2}]]},{i,Length[Data]/10}];
If[Ind≠0,
novoz=ListPlot[DataT[[All,4]],Joined->True, Filling->{1->{1/2,LightRed}} ,Ticks->{ticksx,Automatic},GridLines->{tabx,Automatic},PlotLabel->"New cases from "<>ToString[danIni]<>" to "<>ToString[danEnd]];
hosp=ListPlot[DataT[[All,6]],Joined->True, Filling->{1->{1/2,LightBlue}} ,Ticks->{ticksx,Automatic},GridLines->{tabx,Automatic},PlotLabel->"Hospitalized from "<>ToString[danIni]<>" until "<>ToString[danEnd]];
Print[GraphicsGrid[{{novoz,hosp}},ImageSize->600]];
 ];
pod=Table[{DataT[[i,3]],DataT[[i,4]]},{i,Length[DataT]}] ;  
				(* Daily data of new cases *)
slpod=ListPlot[pod,PlotStyle->{Gray,PointSize[.015]},AxesStyle->axesstyle];   
ill=Table[{DataT[[i,3]],DataT[[i,6]]},{i,delU,Length[DataT]}] ;       
                                       (* Daily data of hospitalized with delay *)slill=ListPlot[ill,PlotStyle->{Gray,PointSize[.015]},AxesStyle->axesstyle];
				(* Cumulative data of new cases *)
POD=pod;
Do[
POD[[i]]={pod[[i,1]],pod[[i,2]]+POD[[i-1,2]]}
,{i,2,Length[pod]}];
sltocN=ListPlot[POD,PlotRange->All,PlotStyle->{Gray,PointSize[.015]}];
			(* Logistic model function of new cases:  DIRECT  *)
f[t_,a_,b_,c_]=a/(1+b Exp[-c t]);
FCompile=Compile[{{a,_Real},{b,_Real},{c,_Real}},
Evaluate[ Sum[(f[POD[[i,1]],a,b,c] -POD[[i,2]])^2,{i,Length[POD]}]]];	
F[a_,b_,c_]:=Apply[FCompile,{a,b,c}];	
dom={{1,2Last[POD][[2]]},{0,5Last[POD][[2]]/Mean[Table[POD[[s,2]],{s,5}]]},{0,10}};
Tmin=Timing[min=DIRECT[F,dom,DMin->10.^(-3),MaxIter->100, Eps->10^(-4)]];
If[Ind≠0,Print[Style["\nParameter estimation of the logistic model function:",Blue,16],"\nDIRECT: ",min]];
a0=min[[1,1]];b0=min[[1,2]]; c0=min[[1,3]];
			     (* Logistic mode function:  NonlinearModelFit *)
nlm=NonlinearModelFit[POD,f[t,a,b,c],{{a,a0},{b,b0},{c,c0}},t];aN1=a/.nlm[[1,2]]; bN1=b/.nlm[[1,2]]; cN1=c/.nlm[[1,2]];
I1={Log[bN1]/cN1,aN1/2};  
tB=Log[bN1/(2+Sqrt[3])]/cN1;
tC=Log[bN1/(2-Sqrt[3])]/cN1;
 B={tB,f[tB,aN1,bN1,cN1]};  
CC={tC,f[tC,aN1,bN1,cN1]};
Print["Parameters: (SuperStar[a]) = ",aN1,"; (SuperStar[b]) = ",bN1,"; (SuperStar[c]) = ",cN1];
If[Ind≠0,Print[Show[sltocN,ListPlot[{I1,B,CC},PlotStyle->{Blue,PointSize[.03],Opacity[.8]}],Plot[{aN1,nlm[t]},{t,0,2I1[[1]]},PlotStyle->Gray],PlotRange->{0,aN1+1},AxesStyle->axesstyle,Ticks->{ticksx,Automatic},GridLines->{tabx,Automatic},PlotLabel->"Logistic model function of new cases ",ImageSize->Medium]];
Print["I = ",Round[I1,-1],";  B = ",Round[B,-1],";  C = ",Round[CC,-1],";   Saturation level: ",Round[aN1,-1]];
];
(* The function of hospitalized phi according to Logistic model-function and prediction; 
Initial appr. by using DIRECT, and then FindMinimum (Newton)*)
phi[t_]:=aN1/(1+bN1 Exp[-cN1 t]) ;   
                            (* Logistic model-function of new cases: DIRECT *)
FCompile=Compile[{{u,_Real},{v,_Real},{del,_Real}},
Evaluate[ Sum[(u+v phi'[ill[[i,1]]-del] -ill[[i,2]])^2,{i,Length[ill]}]]];	
F[u_,v_,del_]:=Apply[FCompile,{u,v,del}];	
dom={{1,10000},{0,200},{0,50}};
Tmin=Timing[min=DIRECT[F,dom,DMin->10.^(-3),MaxIter->100, Eps->10^(-4)]];
If[Ind≠0,Print[Style["\nParameter estimation of the relating model function ψ:",Blue,16],"\nDIRECT: ",min]];
u1=min[[1,1]]; v1=min[[1,2]];del1=min[[1,3]];
							(* NonlinearModelFit *)
nlm=NonlinearModelFit[ill,u+v phi'[t-del],{{u,u1},{v,v1},{del,del1}},t];
u1=u/.nlm[[1,2]];v1=v/.nlm[[1,2]];del1=del/.nlm[[1,2]];Rdel1=Round[del1];
Print["Parameters: (SuperStar[u]) = ",u1,"; (SuperStar[v]) = ",v1,"; (SuperStar[δ]) = ",del1];
If[Ind≠0,Print[Plot[u1+v1 phi'[t-del1],{t,0,Length[DataT]+10},Ticks->{ticksdel,Automatic},GridLines->Automatic,PlotLabel->"Model function ψ",ImageSize->300]]];
If[u1>0 && v1>0 &&del1≥1,
If[Rdel1<22,
PP=Table[Round[ u1+v1 phi'[s-del1],-1],{s,Length[pod]+1,Length[pod]+Rdel1}];
If[Ind≠0,slphid=Plot[u1+v1 phi'[t-del1],{t,0,Last[ill][[1]]+del1},PlotStyle->{Red}];
				(* Characteristic points of the function ψ *)
Print[Style["Inflection point and maximum of the model function ψ given by (8)",Blue,15]];
psi[t_]:=u1+v1 phi'[t-del1];
max=NMaximize[psi[t],t]; t0=t/.max[[2]];T0=Round[{t0,max[[1]]}];
tI=del1+Log[bN1(2-Sqrt[3])]/cN1;
TI=Round[{tI,psi'[tI]}];
LD=Last[Data][[3]];
Print["Inflection point = ",TI,Style[";   date: ",Blue,14],
If[TI[[1]]≤LD,Data[[TI[[1]],{1,2}]],"after last date"]];
Print["maximum = ",T0,Style[";   date: ",Blue,14],If[T0[[1]]≤LD,Data[[T0[[1]],{1,2}]],"after last date"]  ];
						(* Confidence Bands *)
{bands90[ t_],bands95[ t_],bands99[ t_]}=Table[nlm["SinglePredictionBands",ConfidenceLevel->cl],{cl,{.90,.95,.99}}];
Print[Show[ Plot[{nlm[t],bands90[t],bands95[t],bands99[t]},{ t,First[ill][[1]]-5, T},Filling->{2->{1},3->{2},4->{3}},PlotLabel->"Confidence bands",PlotLegends->{"model function ψ", "90% confidence bands","95% confidence bands","99% confidence bands"}],ListPlot[ill],
Ticks->{ticksdel,Automatic},PlotRange->{0,Max[ill[[All,2]]]+15},GridLines->{tabx,Automatic},AxesOrigin->{delU-5,0}]];
Print["Parameters: (SuperStar[u]) = ",u1,"; (SuperStar[v]) = ",v1,"; (SuperStar[δ]) = ",del1];
];
Print[Style["Prognosis of hospitalized for the next ",Blue,14],Style[Round[del1,-1],Red,14],Style[" days, starting on ",Blue,14],Style[danEnd,Red,14],Style[":\n",Blue,14],Style[SL=Table[PP[[s]],{s,Rdel1}],Red,14]];
If[Length[pod]+Rdel1≤ Length[Data],
Print["Actually hospitalized: ",SP=Table[Data[[s,6]],{s,Length[pod]+1,Length[pod]+Rdel1}]];
Print["RMSE = ",rmse=Round[Sqrt[Mean[(PP-SP)^2]],-.1],";   MAPE = ",mape=Round[Mean[Abs[100.(PP-SP)/SP]],-.1],"%"]
],
par=13;Print[Style["GOP (9) does not give a sensible solution (predicting for more than 3 weeks makes no sense):\n(OverHat[δ]) = ",Red,12],Rdel1," > 21"]
],
par=13;Print[Style["GOP (9) does not give a sensible solution: one or several parameters are not positive ",Red,12]]
]
]

DIRECT Modul

Options[DIRECT]={Eps->0,MaxIter->20,DMin->10.^(-2)};
DIRECT[f_,dom_,OptionsPattern[]]:=Module[{n,ff,c,d,dmax,fc,dmaxrazni,pot0,pot1,pot2,temp,ind,novic,novid,novifc,red,TMin,iter=0,dtest,anim={}},
dtest=N[OptionValue[DMin]/Max[(dom[[All,2]]-dom[[All,1]])]];
DIRECT::maxiter="The maximal number of iterations (`1`) has been reached.";
n=Length[dom];
(* Scaling the domain to the unit hypercub *)
ff=Function[N[f@@({##}(dom[[All,2]]-dom[[All,1]])+dom[[All,1]]),50]];
(* Initial center and diameters *)
c={Table[1/2,{i,n}]};
d={Table[1,{i,n}]};
dmax={Max[d[[1]]]};
dmaxrazni=Union[dmax];
fc={ff@@c[[1]]};
TMin={c[[1]],fc[[1]]};
(* MAIN LOOP *) 
Monitor[
While[(dmaxrazni[[1]]>dtest)&&(iter<OptionValue[MaxIter]),
iter++;
(* STEP 1 :searching for potentially optimal hyper-rectangles *)
dmaxrazni=Union[dmax];
(* All data grouped according to d *)
pot0=Table[Select[Range[Length[fc]],dmax[[#]]==dmaxrazni[[i]]&],{i,Length[dmaxrazni]}];
(* Data with minimal f grouped according to d *)
pot1=Table[pot0[[i,Position[fc[[pot0[[i]]]],Min[fc[[pot0[[i]]]]]][[All,1]]]],{i,Length[dmaxrazni]}];
(* Checking other criteria *)
pot2={pot1[[-1]]};
Do[temp=Min[Table[(fc[[pot1[[i,1]]]]-fc[[pot2[[j,1]]]])/(dmax[[pot1[[i,1]]]]-dmax[[pot2[[j,1]]]]),{j,Length[pot2]}]];
If[temp>OptionValue[Eps]&&(i==1||Max[Table[(fc[[pot1[[i,1]]]]-fc[[pot1[[j,1]]]])/(dmax[[pot1[[i,1]]]]-dmax[[pot1[[j,1]]]]),{j,i-1}]]<temp)
,AppendTo[pot2,pot1[[i]]]];
,{i,Length[dmaxrazni]-1,1,-1}];
pot2=Flatten[pot2,1];
(* STEP 2: Dividing potentially optimal hyper-rectangles *)
Do[
(* Dividing directions *)
ind=Flatten[Position[d[[pot]],Max[d[[pot]]]]];
novic=Table[c[[pot]],{2*Length[ind]}];
novid=Table[d[[pot]],{2*Length[ind]}];
d[[pot,ind]]/=3;
dmax[[pot]]/=3;
Do[
novic[[2i-1,ind[[i]]]]-=d[[pot,ind[[i]]]];
novic[[2i,ind[[i]]]]+=d[[pot,ind[[i]]]];
,{i,Length[ind]}];
novifc=Table[ff@@novic[[i]],{i,Length[novic]}];
red=Ordering[Table[Min[{novifc[[2i-1]],novifc[[2i]]}],{i,Length[ind]}]];
Do[novid[[2red[[j]]-1,ind[[red[[i]]]]]]=novid[[2red[[j]],ind[[red[[i]]]]]]/=3,{i,1,Length[red]},{j,i,Length[red]}];
c=Join[c,novic];
d=Join[d,novid];
dmax=Join[dmax,Map[Max,novid]];
fc=Join[fc,novifc];
If[novifc[[2red[[1]]-1]]<TMin[[2]],TMin=N[{novic[[2red[[1]]-1]],novifc[[2red[[1]]-1]]},10]];
If[novifc[[2red[[1]]]]<TMin[[2]],TMin=N[{novic[[2red[[1]]]],novifc[[2red[[1]]]]},10]];
,{pot,pot2}];
(*Print[iter,": ",{TMin[[1]](dom[[All,2]]-dom[[All,1]])+dom[[All,1]],TMin[[2]]}];*)
If[iter==OptionValue[MaxIter],Message[DIRECT::maxiter,OptionValue[MaxIter]];];
(* Figure for the 2D case *)
(*If[n==2,
AppendTo[anim,
Graphics[{{Green,Rectangle[]},Table[{Blue,Rectangle[c[[i]]-d[[i]]/2+d[[i]]/10,c[[i]]+d[[i]]/2-d[[i]]/10]},{i,Length[d]}],Point[c]}]
];
];*)
];
,{iter,Length[pot2]}];
(* Back to the original domain *)
TMin[[1]]=TMin[[1]](dom[[All,2]]-dom[[All,1]])+dom[[All,1]];
(* Printing figures for the 2D case *)
(*If[n==2,Print[ListAnimate[anim,AnimationRunning->False]]];*)
Return[TMin];
]

Data input

SetDirectory[NotebookDirectory[]]; 
Data=ReadList["Croatia.txt"];
(*Print["|Data| = ",Length[Data]]*)
(*Print["First datum: ",First[Data]]*)
(*Print["Last datum: ",Last[Data]]*)

(*       To change the initial datum set  k>1 (default is k=1)              *)
k=1;
Data=Table[Data[[i]],{i,k,Length[Data]}];
Print["|Data| = ",Length[Data]]
Print["First datum: ",First[Data]]
Print["Last datum: ",Last[Data]]

(*      Graphical representation of all new cases             *)
dayIni=First[Data][[{1,2}]]; dayEnd=Last[Data][[{1,2}]];
tabx=Table[10i,{i,150}]; 
ticks=Table[i,{i,0,Length[Data],10}];fig1=ListPlot[Data[[All,4]],Joined->True, Filling->{1->{1/2,LightRed}} ,Ticks->{ticks,Automatic},GridLines->{ticks,Automatic},PlotLabel->"New cases from "<>ToString[dayIni]<>" to "<>ToString[dayEnd],ImageSize->Medium];
(* 7-day moving average of new cases *)
DataAv=Table[N[Sum[Data[[All,4]][[n]],{n,n-3,n+3}]/7],{n,0,Length[Data]-3}];
fig2=ListPlot[DataAv,Joined->True,Ticks->{ticks,Automatic},GridLines->{ticks,Automatic},PlotLabel->"7-day moving average of new cases from "<>ToString[dayIni+{4,0}]<>" to "<>ToString[dayEnd-{2,0}],ImageSize->Medium];
Print[GraphicsGrid[{{fig1,fig2}},ImageSize->800]]

Exe with details

(* INPUT:25 ≤  T ≤ Length[Data], the point which divides training and testing data in Data *)
T=50;
If[25≤ T≤  Length[Data],
DataT=Table[Data[[s]],{s,T}];
Hospitalized[DataT,1],
Print[Style["T must be between 25 and ",Red,12],Length[Data]]];

Exe without details

(* INPUT:25 ≤  T ≤ Length[Data], the point which divides training and testing data in Data *)
T=70;
If[25≤ T≤  Length[Data],
DataT=Table[Data[[s]],{s,T}];
Hospitalized[DataT,0],
Print[Style["T must be between 25 and ",Red,12],Length[Data]]];

All training data sets

(* Input: jIni ≥ 25 and jEnd ≤ Length[Data] *)
jIni=33; jEnd=112;Hosp={};
If[jIni≥  25,
If[jEnd≤   Length[Data]&& jIni<jEnd,
Do[
T=j;
Print[Style["\nT = ",Blue,18],Style[j,Blue,18]];
DataT=Table[Data[[s]],{s,j}];
Hospitalized[DataT,0];
If[par== 13,
rmse=mape=Rdel1=tsp=tpp=-Infinity,
If[j+Rdel1> Length[Data],
rmse=mape=tsp=tpp=-Infinity;
Print[Style["RMSE and MAPE cannot be determined",Red,12]],
tsp=Total[SP]; tpp=Total[PP]],
tsp=Total[SP]; tpp=Total[PP]
  ];
Hosp=Append[Hosp,{rmse,mape,Rdel1,tsp,tpp}];
,{j,jIni,jEnd}];
tab=Table[i,{i,100}];
ticksX=Table[{s,jIni+s-1},{s,1,jEnd-jIni+1,5}];
Print[ListPlot[Hosp[[All,1]],PlotStyle->{Blue,PointSize[.02],Opacity[.5]},PlotRange->{0,1.05Max[Hosp[[All,1]]]},Ticks->{ticksX,Automatic},GridLines->{tab,Automatic},AxesLabel->{"T","RMSE"},PlotLabel->"RMSE"<>ToString[ {jIni,jEnd}],ImageSize->Medium]];
Print[ListPlot[Hosp[[All,2]],PlotStyle->{Magenta,PointSize[.02],Opacity[.5]},PlotRange->{0,1.05Max[Hosp[[All,2]]]},Ticks->{ticksX,Automatic},GridLines->{tab,Automatic},AxesLabel->{"T","%"},PlotLabel->"MAPE"<>ToString[ {jIni,jEnd}],ImageSize->Medium]];
Print[ListPlot[Hosp[[All,3]],PlotStyle->{Blue,PointSize[.02],Opacity[.5]},PlotRange->{0,1.05Max[Hosp[[All,3]]]},Ticks->{ticksX,Automatic},GridLines->{tab,Automatic},AxesLabel->{"T","Days"},PlotLabel->"Delay"<>ToString[ {jIni,jEnd}],ImageSize->Medium]],
Print[Style["Upper bound of jEnd must be ≤ ",Red,12],Length[Data],Style[" and jIni < jEnd",Red,12]]],
Print[Style["Lower bound of jIni must be ≥ 25",Red,12]]];
Print["Gaps in displayed figures occur for those T's for which the GOP (9) does not give a correct solution."]

Some auxiliary figures
To evaluate cells bellow first evaluate “All training data sets”.

(* Estimated and actual numbers of hospitalized in [T,T+(OverHat[δ])], for every T≥ 33 *)
ticksX=Table[{5i-2,30+5i},{i,15}]
ListPlot[{Hosp[[All,4]],Hosp[[All,5]]},PlotStyle->{{Blue,PointSize[.02]},{Red,PointSize[.02]}},AxesLabel->{"T",None},Ticks->{ticksX,Automatic},GridLines->{Table[5i+3,{i,0,15,1}],Automatic},PlotRange->{0,1.05Max[Hosp[[All,4]]]},ImageSize->300]

(* Estimated and actual numbers of hospitalized in [T,T+(OverHat[δ])], for every T≥45 *)
ticksx=Table[{5i+1,45+5i},{i,0,12,1}]
sl1=ListPlot[Table[Hosp[[All,4]][[i]],{i,13,Length[Hosp[[All,2]]]}],PlotStyle->{Blue,PointSize[.02]}];
sl2=ListPlot[Table[Hosp[[All,5]][[i]],{i,13,Length[Hosp[[All,2]]]}],PlotStyle->{Red,PointSize[.02]}];
Print[Show[sl1,sl2,Ticks->{ticksx,Automatic},GridLines->{Table[5i+1,{i,0,15,1}],Automatic},PlotRange->{0,1.05Max[Hosp[[All,5]]]},AxesLabel->{"T",None},AxesLabel->{"T",None},ImageSize->300]]

(* MAPE for T>50 *)
ListPlot[Table[Hosp[[All,2]][[i]],{i,13,Length[Hosp[[All,2]]]}],PlotStyle->{Gray,PointSize[.02]},PlotRange->{0,33},AxesLabel->{"T","MAPE"},Ticks->{ticksx,Automatic},GridLines->{Table[5i+1,{i,0,15,1}],Automatic},PlotRange->{0,1.05Max[Hosp[[All,2]]]},ImageSize->300]

(* Delay for T>50 *)
ListPlot[Table[Hosp[[All,3]][[i]],{i,13,Length[Hosp[[All,2]]]}],PlotStyle->{Gray,PointSize[.02]},PlotRange->{0,30},AxesLabel->{"T","Delay"},Ticks->{ticksx,Automatic},GridLines->{Table[5i+1,{i,0,15,1}],Automatic},PlotRange->{0,1.05Max[Hosp[[All,3]]]},ImageSize->300]

(* Estimated and actual numbers of hospitalized in [T,T+(OverHat[δ])], for every T≥ 33 *)
sl1=Show[Table[ListPlot[{{i,Hosp[[All,4]][[i]]}},PlotStyle->{Blue,PointSize[Hosp[[All,3]][[i]]/900],Opacity[.5]}],{i,Length[Hosp[[All,4]]]}]];
sl2=Show[Table[ListPlot[{{i,Hosp[[All,5]][[i]]}},PlotStyle->{Red,PointSize[Hosp[[All,3]][[i]]/900],Opacity[.5]}],{i,Length[Hosp[[All,5]]]}]];
Show[sl1,sl2,GridLines->{Table[5i+3,{i,0,15,1}],Automatic},AxesLabel->{"T",None},Ticks->{ticksX,Automatic},PlotRange->{{0,Length[Hosp[[All,4]]]},{0,Max[Drop[Hosp[[All,5]],1]]+10000}},ImageSize->500]