library(forecast)
library(tseries)
library(readxl) 
library(lubridate)
library(bsts)
library(dplyr)
library(ggplot2)
library(zoo)
library(reshape2)
library(CausalImpact)
library(robCompositions)

#######################################################################################
setwd("C:/Research papers/COVID 2019/COVID Pakistan")
x=read.csv(file = "COVID192.csv", header = T)
#x=ts(x,frequency = 1,start = c(26/02/2020,1))
x
#x[, drop=TRUE]
#x2 <- x2[, drop=TRUE];x2
#x2 <- as.numeric(as.character(x2));x2

x1=x[,1]      #Date
x2=x[,9];x2      #New Cases
x3=x[,10];x3     #Cum. Cases
x4=x[,11];x4     #New Deaths
x5=x[,12];x5     #Cum. Deaths
x7=x[,14];x7     #Cum. Recoveries
x8=x[,15];x8     #New Tests
x9=x[,16];x9     #Cum. Tests
x10=x3-x7-x5;x10
x11=x3/x9;x11    #Ratio C.Cases to C.Tests
x12=x5/x3;x12    #Ratio C.Deaths to C.Cases
x13=x7/x3;x13    #Ratio C.Recoveries to C.Cases
x14=x10/x3;x14   #Ratio C.Active Cases to C.Cases
x15=x11[1:73]
x16=x12[22:73]
x17=x13[10:73]
x18=x11[74:137]
x19=x12[74:137]
x20=x13[74:137]
gm(x15)
gm(x16)
gm(x17)
gm(x18)
gm(x19)
gm(x20)

############################Plot of Ratios ##################################

plot(x11,type = "l",lwd=2.5,col= "green",bg="red", pch=21,cex=1.6,lty="22", xlab="Days",
     ylab="Ratios" , main = "", col.main="orange4",cex.main=1.2, cex.lab=1.2, cex.axis=1.2,
     ylim=range(x11,x12,x13,x14))
points(x12,  type = "l",lwd=2.5,col= "blue",bg="blue",pch=21,cex=1.6, lty=4)
points(x13, type = "l", lwd=2.5,col=  "red",bg="darkmagenta",pch=21,cex=1.6, lty=6)
points(x14, type = "l", lwd=2.5,col= "yellow",bg="green3",pch=21,cex=1.6, lty=5)
legend("topright",legend=c("RCT", "RDC","RRC", "RAC"), 
       lwd=3,
       col=c("green","pink3", "darkseagreen3","orange"), 
       pt.bg=c("red","blue","darkmagenta","green3"),
       pch=c(21,21,21,21),lty=c(1,1,1,1), pt.cex=1.8, cex=1.1)
mtext("(A)", side=3, col="green4",cex=1.4,line=-1.4,lwd=5.5)

##############################################################################
######################## No. of Cases ######################################

dt1 <- seq(as.Date("2020-02-26"), as.Date("2020-07-11"), by = "days")
ts <- zoo(x3, dt1)

### Run the bsts model
ss <- AddLocalLevel(list(), ts)
ss <- AddSeasonal(ss, ts, nseasons = 7)
bsts.model <- bsts(ts, state.specification = ss, niter = 20000, ping=0, seed=2016)
### Get a suggested number of burn-ins
burn <- SuggestBurn(0.1, bsts.model)
### Predict
p <- predict.bsts(bsts.model, horizon = 30, burn = burn, quantiles = c(.025, .975))
w1=c(as.numeric(p$mean))
BC12=c(x3,w1)
w10=c(as.numeric(-colMeans(bsts.model$one.step.prediction.errors[-(1:burn),])+x3))
y1=ts(x3,frequency = 1,start = c(26/02/2020,1))
F1=auto.arima(y1)
F12=forecast(F1,h=30)
w2=as.data.frame(F12)
w12=w2[,1]
AC12=c(x3,w12)
w=c(fitted(F1))
#MAPE=mean(abs(y1-w)/y1)
RMSE=(mean((y1-w)^2))^.5
MAE=mean(abs(y1-w))
RMSPE=(mean(((y1-w)/y1)^2))^.5
RMdSPE=(median(((y1-w)/y1)^2))^.5

##############################################################################
############################## No. of deaths #################################

dt1 <- seq(as.Date("2020-02-26"), as.Date("2020-07-11"), by = "days")
ts <- zoo(x5, dt1)


### Run the bsts model
ss <- AddLocalLevel(list(), ts)
ss <- AddSeasonal(ss, ts, nseasons = 7)
bsts.model <- bsts(ts, state.specification = ss, niter = 20000, ping=0, seed=2016)
### Get a suggested number of burn-ins
burn <- SuggestBurn(0.1, bsts.model)
### Predict
p <- predict.bsts(bsts.model, horizon = 30, burn = burn, quantiles = c(.025, .975))
w3=c(as.numeric(p$mean))
BC22=c(x5,w3)
w11=c(as.numeric(-colMeans(bsts.model$one.step.prediction.errors[-(1:burn),])+x5))
y2=ts(x5,frequency = 1,start = c(26/02/2020,1))
F2=auto.arima(y2)
F22=forecast(F2,h=30)
w4=as.data.frame(F22)
w22=w4[,1]
AC22=c(x5,w22)
y21=y2[23:137]
w=c(fitted(F2))
w21=w[23:137]
#MAPE=mean(abs(y1-w)/y1)
RMSE=(mean((y21-w21)^2))^.5
MAE=mean(abs(y21-w21))
RMSPE=(mean(((y21-w21)/y21)^2))^.5
RMdSPE=(median(((y21-w21)/y21)^2))^.5

##############################################################################
############################ No. of Recoveries ###############################

dt1 <- seq(as.Date("2020-02-26"), as.Date("2020-07-11"), by = "days")
ts <- zoo(x7, dt1)

### Run the bsts model
ss <- AddLocalLevel(list(), ts)
ss <- AddSeasonal(ss, ts, nseasons = 7)
bsts.model <- bsts(ts, state.specification = ss, niter =20000, ping=0, seed=2016)
summary(bsts.model)

### Get a suggested number of burn-ins
burn <- SuggestBurn(0.1, bsts.model)
### Predict
p <- predict.bsts(bsts.model, horizon = 30, burn = burn, quantiles = c(.025, .975))
w5=c(as.numeric(p$mean))
BC32=c(x7,w5)
w12=c(as.numeric(-colMeans(bsts.model$one.step.prediction.errors[-(1:burn),])+x7))
y3=ts(x7,frequency = 1,start = c(26/02/2020,1))
F3=auto.arima(y3)
F32=forecast(F1,h=30)
w6=as.data.frame(F32)
w32=w6[,1]
AC32=c(x7,w32)
y31=y3[11:137]
w=c(fitted(F3))
w31=w[11:137]
#MAPE=mean(abs(y1-w)/y1)
RMSE=(mean((y31-w31)^2))^.5
MAE=mean(abs(y31-w31))
RMSPE=(mean(((y31-w31)/y31)^2))^.5
RMdSPE=(median(((y31-w31)/y31)^2))^.5

##############################################################################
############################ No. of Active Cases #############################

dt1 <- seq(as.Date("2020-02-26"), as.Date("2020-07-11"), by = "days")
ts <- zoo(x10, dt1)

### Run the bsts model
ss <- AddLocalLevel(list(), ts)
ss <- AddSeasonal(ss, ts, nseasons = 7)
bsts.model <- bsts(ts, state.specification = ss, niter = 20000, ping=0, seed=2016)
### Get a suggested number of burn-ins
burn <- SuggestBurn(0.1, bsts.model)
### Predict
p <- predict.bsts(bsts.model, horizon = 30, burn = burn, quantiles = c(.025, .975))
w6=c(as.numeric(p$mean))
BC42=c(x10,w6)
w13=c(as.numeric(-colMeans(bsts.model$one.step.prediction.errors[-(1:burn),])+x10))
y4=ts(x10,frequency = 1,start = c(26/02/2020,1))
F2=auto.arima(y4)
F22=forecast(F2,h=30)
w4=as.data.frame(F22)
w22=w4[,1]
AC22=c(x10,w22)
y41=y4[4:137]
w=c(fitted(F2))
w41=w[4:137]
RMSE=(mean((y41-w41)^2))^.5
MAE=mean(abs(y41-w41))
RMSPE=(mean(((y41-w41)/y41)^2))^.5
RMdSPE=(median(((y41-w41)/y41)^2))^.5

###########################################################################
###################### Plot for combine forecasts #########################

BC112=BC12[24:167]
w110=w10[24:137]
BC122=BC22[24:167]
w111=w11[24:137]
BC132=BC32[24:167]
w112=w12[24:137]
BC142=BC112-BC132
w113=w110-w112

dat <- list(BC112,w110,BC122,w111,BC132,w112,BC142,w113)
plot(unlist(dat),type="n",log="y",xlim=c(1,max(sapply(dat,length))),xlab="Days",
     ylab="No. of Cases" , main = "", col.main="orange4",cex.main=1.2, cex.lab=1.2, cex.axis=1.2)
mapply(lines,dat,col=c("cyan","magenta","green","red","blue","yellow","forestgreen","yellowgreen"),pch=8,lty=1,lwd=5.5)
legend("bottomright",legend=c("Forecasted Cases", "Fitted Cases","Forecasted Deaths", "Fitted Deaths", "Forecasted Recoveries","Fitted Recoveries","Forecasted Active Cases", "Fitted Active Cases"), 
       lwd=2.5,
       col=c("cyan","magenta","green","red","blue","yellow","forestgreen","yellowgreen"), 
       pt.bg=c("cyan","magenta","green","red","blue","yellow","forestgreen","yellowgreen"),
       pch=c(21,21,21,21,21,21),lty=c(5,5,5,5,5,5), pt.cex=1.2, cex=1)
mtext("(B)", side=3, col="blue",cex=1.4,line=-1.4,lwd=5.5)


######################################################################################
################### Seasonality Trend and Regression##################################

########################### Total No of Cases ########################################

dt1 <- seq(as.Date("2020-02-26"), as.Date("2020-07-11"), by = "days")
ts <- zoo(x3, dt1)
ts
ss <- AddLocalLevel(list(), ts)
ss <- AddSeasonal(ss, ts, nseasons = 7)
bsts.reg <- bsts(x3 ~x9 , state.specification = ss, data =
                   ts, niter = 500, ping=0, seed=2016)
burn <- SuggestBurn(0.1, bsts.reg)
components.withreg <- cbind.data.frame(
  colMeans(bsts.reg$state.contributions[-(1:burn),"trend",]),
  colMeans(bsts.reg$state.contributions[-(1:burn),"seasonal.7.1",]),
  colMeans(bsts.reg$state.contributions[-(1:burn),"regression",]),
  as.Date(dt1))  
names(components.withreg) <- c("Trend", "Seasonality", "Regression", "Date")
components.withreg <- melt(components.withreg, id.vars="Date")
names(components.withreg) <- c("Date", "Component", "Value")

ggplot(data=components.withreg, aes(x=Date, y=Value)) + geom_line() + 
  theme_bw() + theme(legend.title = element_blank()) + ylab("") + xlab("") + 
  facet_grid(Component ~ ., scales="free") + guides(colour=FALSE) + 
  theme(axis.text.x=element_text(angle = 0, hjust = 1,size=13,color="black"),
        axis.text.y=element_text(angle = 0, hjust = 1,size=13,color="black"),
        axis.line = element_line(colour = "green3"),
        panel.border =element_rect(fill=NA, color='green3', size=1.3, linetype="solid") )+
  geom_point(color='green',size=2.2)+geom_line(color='red',size=.74)


########################################################################################
############################ Causal Impact of Lockdown #################################

set.seed(1)
time.points <- seq.Date(as.Date("2020-02-26"), by = 1, length.out = 137)
data <- zoo(cbind(x3,x9), time.points)
head(data)
pre.period <- as.Date(c("2020-02-26", "2020-05-08"))
post.period <- as.Date(c("2020-05-09", "2020-07-11"))
impact <- CausalImpact(data, pre.period, post.period)
plot(impact,"original")
summary(impact)