## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
## Registered S3 methods overwritten by 'forecast':
##   method             from    
##   fitted.fracdiff    fracdiff
##   residuals.fracdiff fracdiff

Todo modelo predictivo debe probarse con datos reales, para determinar su validez. Una forma de medir el grado de error es el error medio cuadrático, que se calcula

\[ EMC= \sum_i^k \frac{(y_i-\hat{y_i})^2}{n} \]

Donde \(y_i\) es el dato observado en el momento \(i\), entanto que \(\hat{y_i}\) es valor predicho. En nuestros ejemplos, como todo lo hacemos sobre el logaritmo de las curvas, se calcular el EMC sobre los valores al logaritmo

Probaremos la adecuación de los tres modelos, usando el modelo calculado el día 28 de Marzo. Tal como se puede ver, justo al día siguiente de la predicción se observa un cambio en la curva que dejo todos los modelos fuera del rango predictivo.

calcular.mse<-function(x,y) {
  mean((log(x$casos)-log(y$casos))^2)
}

dato.al.28.3<-datos.pais$Chile[datos.pais$Chile$fecha<=43917,]
dato.al.3.4<-datos.pais$Chile[datos.pais$Chile$fecha<=43924,]

obtener.prediccion<-function(x.parc,x.comp,min.casos) {
  dato.obs<-data.frame(
    dia=x.comp$dia,
    casos=x.comp$casos,
    li=x.comp$casos, 
    ls=x.comp$casos,
    tipo="Observado",
    pais="Chile"
)
  prediccion<-prediccion.casos(list(Chile=x.parc),min.casos = min.casos)
  out<-list()
  out$dato.obs<-dato.obs
  out$dato.obs.parc<-tail(x.comp,7)
  out$p.tar21<-tail(prediccion$tar2$Chile,7)
  out$p.arima1<-tail(prediccion$arima$Chile,7)
  out$p.expon1<-tail(prediccion$exp$Chile,7)
  out$para.ver<-rbind(out$dato.obs, out$p.tar21,out$p.arima1,out$p.expon1)
  out
}
library(ggplot2)

pred.250<-obtener.prediccion(dato.al.28.3, dato.al.3.4, 250)
ggplot(pred.250$para.ver, aes(x=dia, y=casos, color=tipo, ymin=li,ymax=ls))+geom_point()+geom_line()+geom_ribbon(alpha=0.1)+scale_y_continuous(trans="log10")+geom_hline(yintercept = 250)

Podemos ver que, a diferencia de la evaluación del día 28 que consideró el modelo del día 21 que consideraba 20 casos (revisar Github), donde el MSE de T+AR(1) era de 0.003 para el periodo, ahora ascendió a 0.21, siendo el ARIMA y el Exponencial prácticamente idénticos en su capacidad predictiva.

data.frame("T+AR(1)"=calcular.mse(pred.250$dato.obs.parc, pred.250$p.tar21),
           "ARIMA(1,0,0)"=calcular.mse(pred.250$dato.obs.parc, pred.250$p.arima1),
           "Exponencial"=calcular.mse(pred.250$dato.obs.parc, pred.250$p.expon1)
           )

Viendo la situación en contexto, podemos ver que la razón de la falla podría ser que solo consideramos la curva sobre 250 casos. Si consideramos la curva desde los 50 casos, el modelo de predicción T+AR(1) presenta un mucho mejor comportamiento.

pred.50<-obtener.prediccion(dato.al.28.3, dato.al.3.4, 50)
ggplot(pred.50$para.ver, aes(x=dia, y=casos, color=tipo, ymin=li,ymax=ls))+geom_point()+geom_line()+geom_ribbon(alpha=0.1)+scale_y_continuous(trans="log10")+geom_hline(yintercept = 50)

Claramente, al utilizar una serie más larga el modelo ARIMA y exponencial se dañan, pero el de tendencia funciona mucho mejor.

data.frame("T+AR(1)"=calcular.mse(pred.50$dato.obs.parc, pred.250$p.tar21),
           "ARIMA(1,0,0)"=calcular.mse(pred.50$dato.obs.parc, pred.250$p.arima1),
           "Exponencial"=calcular.mse(pred.50$dato.obs.parc, pred.250$p.expon1)
           )

Si utilizamos el corte del periodo anterior, con 20 casos, el modelo T+AR(1) subestima un poco, pero manteniendo un MSE similar.

pred.20<-obtener.prediccion(dato.al.28.3, dato.al.3.4, 20)
ggplot(pred.20$para.ver, aes(x=dia, y=casos, color=tipo, ymin=li,ymax=ls))+geom_point()+geom_line()+geom_ribbon(alpha=0.1)+scale_y_continuous(trans="log10")+geom_hline(yintercept = 50)

Claramente, al utilizar una serie más larga el modelo ARIMA y exponencial se dañan, pero el de tendencia funciona mucho mejor.

data.frame("T+AR(1)"=calcular.mse(pred.20$dato.obs.parc, pred.250$p.tar21),
           "ARIMA(1,0,0)"=calcular.mse(pred.20$dato.obs.parc, pred.250$p.arima1),
           "Exponencial"=calcular.mse(pred.20$dato.obs.parc, pred.250$p.expon1)
           )

##Conclusión##

Durante la próxima semana, la predicción se realizará usando un modelo que parta del día 50, no desde el caso 250.