2 Iterative estimation of parameters and outof sample prediction.

This is interesting to monitor the evolution of covid-19.

If model is reliable, then also parameter evolution is reliable.

Further we check its out-of-sample prediction of cases since beginning of march. This is more reliable than the above goodness of fit.

ndays=sum(cvirus$fecha>"2020-03-10")
evpars=data.frame(NULL)
windows=list(1:(n-ndays+1))
for(i in (n-ndays+2):n) windows=c(windows,list(1:i))
nw=length(windows)
for(i in 1:nw){
  fit = sampling(m1,
           data=list(y=cvirus$casos[windows[[i]]],
                     n=length(windows[[i]]),np=2),
           iter=niter,chains = 4,
            init = init_f(),
#           control = list(adapt_delta = 0.99,max_treedepth=50),
           seed = 11,refresh=0)
  post.par=extract(fit,pars=c(int.par,"ypf"))
  for(j in 1:length(int.par)){
    evpars=rbind(evpars,data.frame(day=cvirus$fecha[max(windows[[i]])],
                              value=post.par[[j]],
                              param=int.par[j]))
  }
  evpars=rbind(evpars,data.frame(day=cvirus$fecha[max(windows[[i]])],
                              value=post.par[[j+1]][,1],
                              param="outpred"))

  cat("i=",i,"/",nw," - ")
}
save(evpars,file="evpars.RData")

2.1 Evolution of parameters (using data since the beginning)

load(file="evpars.RData")
ggplotly(ggplot(evpars[evpars$param!="outpred",], aes(x=day, y=value, colour=param)) + geom_smooth()+
  geom_vline(xintercept=as.Date("2020-03-10"),linetype="dashed",color = "red", size=2)+
  geom_hline(yintercept=0)+
  xlab("Days")+  ylab("Posterior mean and 95% C.I."))

Reaching a peak means \(\alpha_2<0\) and \(\beta_2<0\), while disappear of the covid means all parameters less than 0.

2.2 Out of sample prediction

outpred=evpars[evpars$param=="outpred",]
outpred=aggregate(value~day,outpred,
                             function(xx) c(quantile(xx,0.025),mean(xx),quantile(xx,0.975)))
outpred=data.frame(do.call(cbind,outpred))
colnames(outpred)=c("day","inf","mean","sup")
outpred$predday=sort(unique(evpars$day))+1
ggplot(outpred, aes(x=predday, y=mean)) + geom_line() +
geom_point(data=cvirus[cvirus$fecha%in%outpred$day,],
           aes(x=fecha,y=casos),color="red")+
  geom_errorbar(aes(ymin=inf, ymax=sup), width=.2,position=position_dodge(0.05))+
  xlab("Days")+ ylab("Posterior outofsample mean and 95% C.I.")