1. 3. calculate forecast of covid-19
   4. using sarima
   6. 14.4.2022
   7. v 0000.0002

1. install.packages("astsa", repos ="https://ftp.acc.umu.se/mirror/CRAN/")
2. install.packages("MLmetrics", repos ="https://ftp.acc.umu.se/mirror/CRAN/")

2. install.packages("rvest")
3. install.packages("readtext")
4. install.packages("stringi")
5. install.packages("datamart")
6. install.packages("XML")
7. install.packages("svglite")
8. install.packages("ggplot2")
9. install.packages("tidyr")
10. install.packages("stringr")
11. install.packages("stringi")
12. install.packages("tibble")

Sys.setlocale("LC_ALL","Finnish") options(encoding = "UTF-8")

library(ggplot2) library(svglite)

library(rvest) library(readtext) library(stringi) library(stringr) library(datamart) library(XML) library(jsonlite) library(rjstat)

library(tibble) library(caTools) library(mgcv) library(repmis) library(lubridate) library(tidyverse) library(tidyr) library(dplyr)

1. library(covid19.analytics)

library(R0) library(EpiEstim)

1. library(prophet)

library(forecast) library(astsa) library(MLmetrics)

1. 1. choices
   2. 1 finnish wiki data, 2 aggregated cases data
   3. 3 solanpaa finnish data 4 thl cube json data

load_data_from=3

yala=0 yyla=10000

beginday1='02/11/2021'

forecastendday1<-"2022/07/01"

1. 1. 1 from finnish wiki, 2 cases from net, 3 from net 2

plottaa=1 ## must be 1 tulosta_svg=1 # plot to out svg 0, 1 of 2

tulosfilee1="/Users/himot/akor1/sarimaa1.svg"

1. 1. ggolot smooth curves pars
2. spanni=0.2

spanni=0.5

metodi="loess"

1. 1. NOTE date limits change this

1. datelimits1=c('1/3/2020', '9/11/2020')
   1. display date limits

1. 1. this

today=Sys.Date()-15

1. or this

1. today=Sys.Date()

1. print(today)

today1=format(today, "%d/%m/%Y") today2=format(today, "%Y/%m/%d")

1. print(today1)

print(today2)

1. stop(-1)

datelimits1=c(beginday1, today1)

1. 1. dates of dataset

paivat1=seq(as.Date("2020/4/1"), as.Date(today2), "days")

calculate_r0 <- function(time1, time2, val1, val2) { td=time2-time1 gr0<-log(val2/val1) gr=gr0/td td = log(2)/gr tau<-5.0 k<-log(2.0)/td r0<-exp(k*tau) return(r0) }

moving_average <- function(x, w, FUN, ...) {

 if (w < 1) {
   stop("Window length: mustbe greater than 0")
 }
 output <- x
 for (i in 1:length(x)) {
   lower_bound <- i - w + 1
   if (lower_bound < 1) {
     output[i] <- NA_real_
     ## !!! assume NA 0
     output[i] <- 0
   } else {
     output[i] <- FUN(x[lower_bound:i, ...])
   }
 }
 return (output)

}

calculate_multiple_r0 <- function(daata1) {

 lenu1<-length(daata1)
  
 daata2<-1:lenu1
 

for (n in 2:lenu1){ valju1=daata1[n-1] valju2=daata1[n] timex1=0 timex2=1

r0<-calculate_r0(0, 1, valju1, valju2) daata2[n]<-r0

#print (r0) } return(daata2)

}

load_data_from_finnish_wiki<-function() {

url1="https://fi.wikipedia.org/wiki/Suomen_koronaviruspandemian_aikajana" destfile1="./ward0.txt"

download.file(url1, destfile1) texti000<-readtext(destfile1) texti0<-texti000$text

etsittava1="1. huhtikuuta 2020 alkaen" len1=nchar(texti0) k1=regexpr(pattern=etsittava1, texti0) k1b=len1-k1 texti1=strtail(texti0,k1b) sink("out1.txt") print (texti1) sink()

k2=regexpr(pattern=etsittava2, texti1) texti2=strhead(texti1,k2)

sample1<-minimal_html(texti2) tabu1 <- html_table(sample1, fill=TRUE)1 colnames(tabu1) <- c("V1","V2", "V3","V4", "V5","V6", "V7","V8" )

1. print(tabu1)

sairaalassa00<-tabu1$V4 sairaalassa=as.integer(sairaalassa00)

teholla00<-tabu1$V5 teholla=as.integer(teholla00)

uusiatapauksia00<-tabu1$V3 uusiatapauksia0<-gsub(" ", "", uusiatapauksia00) uusia_tapauksia=as.integer(uusiatapauksia0)

uusiakuolleita00<-tabu1$V7 uusiakuolleita1=as.integer(uusiakuolleita00)

uusiakuolleita2<-uusiakuolleita1 uusiakuolleita2[uusiakuolleita2<0]<-0 uusia_kuolleita<-uusiakuolleita2

toipuneita00<-tabu1$V8 toipuneita01<-gsub(" ", "", toipuneita00) toipuneita0<-gsub("[^0-9.-]", "", toipuneita01) toipuneita=as.integer(toipuneita0)

tapauksia00<-tabu1$V2 tapauksia01<-gsub(" ", "", tapauksia00) tapauksia0<-gsub("[^0-9.-]", "", tapauksia01)

tapauksia=as.integer(tapauksia0)

kuolleita00<-tabu1$V6 kuolleita=as.integer(kuolleita00)

aktiivisia_tapauksia=tapauksia-kuolleita-toipuneita

1. print (paivat1)
2. print (teholla)
3. print (sairaalassa)
4. print (tapauksia)
5. print (kuolleita)
6. print (toipuneita)
7. print (uusia_tapauksia)
8. print (uusia_kuolleita)
9. plot(paivat1,aktiivisia_tapauksia)

1. xy<-data.frame(paivat1, sairaalassa)

xy<-data.frame(paivat1, uusia_tapauksia)

names(xy)<-c("daate3", "dailycases3")

xyz<-data.frame(paivat1, sairaalassa, teholla) dfout1<-data.frame(paivat1, aktiivisia_tapauksia, uusia_tapauksia, sairaalassa, teholla, uusia_kuolleita )

names(dfout1)<-c("Pvm", "Aktiivisia_tapauksia","Uusia_tapauksia", "Sairaalassa", "Teholla", "Uusia_kuolleita")

write.csv2(dfout1, "./sairaalassa.csv",row.names=FALSE )

return(xy) }

load_data_from_aggregated<-function() {

1. 2. fetch the data

dfine <- read.csv(file = 'https://datahub.io/core/covid-19/r/countries-aggregated.csv')

1. head(dfine)
2. class(dfine)

1. tail(dfine, 5)

dfinland <- dfine[ which(dfine$Country=='Finland'), ]

1. head(dfinland)

kols <- c("Date", "Confirmed","Recovered","Deaths")

tapaukset <- dfinland[kols]

1. head(tapaukset)

len1=nrow(tapaukset)

1. len1

len2=len1-1

len3=len2

confirmed<-tapaukset$Confirmed deaths<-tapaukset$Deaths

dailycases <- vector() dailycases <- c(dailycases, 0:(len2)) dailydeaths <- vector() dailydeaths <- c(dailydeaths, 0:(len2))

m=0 dailycases[1]<-tapaukset$Confirmed[1]

1. dailydeaths[1]<-tapaukset$Deaths[1]

dailydeaths[1]<-0

1. confirmed
2. deaths

m=1 for(n in 2:(len3+1)) {

a<-confirmed[n] b<-confirmed[m] #print (a) #print (b) cee<- (a-b) #print(cee) dailycases[n]=cee m=m+1 }

mm=1 for(nn in 2:(len3+1)) {

aa<-deaths[nn] bb<-deaths[mm] #print ("_") #print (aa) #print (bb) ceb=aa-bb #if (ceb<0) ceb=0 #print(ceb) dailydeaths[nn]=ceb mm=mm+1 }

1. deaths

1. dailycases

1. dailydeaths

dfout1<-dfinland

1. print(nrow(dfinland))
2. print(length(dailydeaths))

dfout1 <- cbind(dfout1, data.frame(dailycases)) dfout1 <- cbind(dfout1, data.frame(dailydeaths))

1. head(dfout1)

dfout2<-within(dfout1, rm(Country))

names(dfout2) <- c('Date','Confirmed','Recovered','Deaths', 'DailyConfirmed','DailyDeaths')

1. head(dfout2)

write.csv2(dfout2, "/Users/himot/akor1/finland_data1.csv");

daate1<-dfout2$Date dailydeaths1<-dfout2$DailyDeaths dailycases1<-dailycases

1. daate1

1. daate2<-gsub("2020-", "", daate1)

daate2<-daate1

leenu<-length(daate2)

1. alkupvm<-50

alkupvm<-1

daate3<-daate2[alkupvm:leenu] dailydeaths3<-dailydeaths1[alkupvm:leenu] dailycases3<-dailycases1[alkupvm:leenu]

1. daate3
2. dailydeaths3

# barplot(dailydeaths3, main="Koronaviruskuolemat päivittäin vuonna 2020",
# names.arg=daate3) 
 
     
 dataf1 <- data.frame("Date" = daate3, "Paivitt_kuolemat"=dailydeaths3)

1. str(dataf1)

  dataf2 <- data.frame("Date" = daate3, "Paivitt_tapaukset"=dailycases3)

1. str(dataf2)

write.csv(dataf1, "/Users/himot/akor1/dailydeaths1.csv", row.names=T)
write.csv(dataf2, "/Users/himot/akor1/dailycases1.csv", row.names=T)

indf1 <- read.csv(file = '/Users/himot/akor1/dailycases1.csv')
#head(indf1)

cases1<-indf1$Paivitt_tapaukset
dates1<-indf1$Date

len1=length(cases1)

dates2<-as.Date(dates1)
paivat<-1:len1	

xy<-data.frame(daate3, dailycases3)

}

calculate_r0_with_r0<-function(xy2) {

## calculate r0 w/r0 package dates<-as.Date(xy2$Dates) cases<-as.integer(xy2$Cases)

cases[is.na(cases)] <- 1 cases[(cases<0)] <- cases*-1 cases[cases==0] <- 1

nummeros<-1:length(dates) num<-cases #names<-nummeros names<-dates lenu=length(dates)

bekini=as.Date(dates[1]) enti=as.Date(dates[lenu])

#print(bekini) #print(enti)

#stop(-1)

#enti=lenu

#bekini=enti*0+1

#enti=as.integer(enti) #bekini=as.integer(bekini)

df1 <- setNames(num, names)

mGT<-generation.time("gamma", c(3, 1.5)) #TD <- est.R0.TD(df1, mGT, begin=1, end=length(dates), nsim=200) #TD <- est.R0.TD(df1, mGT, begin=bekini, end=enti, nsim=200) TD <- est.R0.TD(df1, mGT, begin=bekini, end=enti, nsim=200) TD.5D <- smooth.Rt(TD, 5) paivat1<-TD.5D$epid$t paivat2<-as.Date(paivat1) r0t1<-TD.5D$R conf1<-TD.5D$conf.int

xypaluu<-data.frame(paivat1,r0t1) names(xypaluu)<-c("paivat","r0") return(xypaluu) }

calculate_r0_with_epiestim<-function(xy2) {

## calculate r0 w/r0 package dates<-as.Date(xy2$Dates) cases<-as.integer(xy2$Cases) nummeros<-1:length(dates) num<-cases #names<-nummeros names<-dates lenu=length(dates)

cases[is.na(cases)] <- 1 cases[(cases<0)] <- cases*-1 cases[cases==0] <- 1

incid<-cases

bekini=as.Date(dates[1]) enti=as.Date(dates[lenu])

config<-make_config( list(mean_si = 2.6,std_si = 1.5) )

res<-estimate_R(incid,method="parametric_si", config = config)

#plot(res) resr<-res$R

str(resr)

meanr<-resr$Mean medianr<-resr$Median quantile95<-resr$Quantile.0.95 quantile05<-resr$Quantile.0.05 quantile75<-resr$Quantile.0.75 quantile25<-resr$Quantile.0.25 meanr

daydexes<-resr$t_start

daydexes

#plot(daydexes, meanr)

dayss<-as.Date(dates[daydexes])

print (dayss) #stop(-1) #plot(dayss, meanr)

xypaluu<-data.frame(dayss,meanr) names(xypaluu)<-c("paivat","r0") return(xypaluu) }

calculate_r0_with_simple_exponent_moving_average<-function(xy2, madays1, madays2) {

## calculate r0 w/r0 package dates<-as.Date(xy2$Dates) cases<-as.integer(xy2$Cases) nummeros<-1:length(dates) num<-cases #names<-nummeros names<-dates lenu=length(dates)

cases[is.na(cases)] <- 1 cases[(cases<0)] <- cases*-1 cases[cases==0] <- 1

# compute a MA(7) ma1<-moving_average(cases,madays1,mean) r0t1<-calculate_multiple_r0(ma1) r0avg1<-moving_average(r0t1, madays2, mean) xypaluu<-data.frame(dates,r0t1)

#plot(r0t1) #print (r0t1) #stop(-1)

names(xypaluu)<-c("paivat","r0")

return(xypaluu)

}

forecast_profet<-function(xy2, futuredays) { ## calculate r0 w/r0 package ds<-as.Date(xy2$Dates) y<-as.integer(xy2$Cases) nummeros<-1:length(ds)

lenu=length(ds)

df<-data.frame(ds,y)

m <- prophet(df)

future <- make_future_dataframe(m, periods = futuredays)

forecast <- predict(m, future)

   #str(future)
   #str(forecast)
   
   futu_days=future$ds
   futu_trendi=forecast$trend
   futu_trendi_upper=forecast$trend_upper
   futu_trendi_lower=forecast$trend_lower
   futu_yhat=forecast$yhat
   futu_yhat_upper=forecast$yhat_upper
   futu_yhat_lower=forecast$yhat_lower
   futu_weekly=forecast$weekly
   futu_weekly_upper=forecast$weekly_upper
   futu_weekly_lower=forecast$weekly_lower
   
   xypaluu<-data.frame(as.Date(futu_days),futu_yhat)

# xypaluu<-data.frame(as.Date(futu_days),futu_weekly) #plot(r0t1) #print (r0t1) #stop(-1)

names(xypaluu)<-c("paivat","r0") return(xypaluu)

}

lataa_thl_tapaukset_kuolleet<-function() {

   url1<-"https://sampo.thl.fi/pivot/prod/fi/epirapo/covid19case/fact_epirapo_covid19case.json?row=measure-492118&column=dateweek20200101-508804L"

cube1 <- fromJSONstat(url1, naming = "label", use_factors = F, silent = T) res01 <- cube11 #res00 url2<-"https://sampo.thl.fi/pivot/prod/fi/epirapo/covid19case/fact_epirapo_covid19case.json?row=measure-444833&column=dateweek20200101-508804L" cube2 <- fromJSONstat(url2, naming = "label", use_factors = F, silent = T) res02 <- cube21 #res02

#stop (-1) paiva=as.Date(res01$dateweek20200101) kuolleet=as.integer(res01$value) tapaukset=as.integer(res02$value)

kuolin_prosentit=kuolleet/tapaukset kuolin_prosentit=kuolin_prosentit*10000 kuolin_prosentit=as.integer(kuolin_prosentit) kuolin_prosentit=as.double(kuolin_prosentit) kuolin_prosentit=kuolin_prosentit/100.0

#print (paiva) #print (kuolleet) #stop(-1) #print (tapaukset) #print (kuolin_prosentit )

df1<-data.frame(paiva,tapaukset, kuolleet, kuolin_prosentit)

names(df1)<-c("Paiva", "Tapauksia", "Kuolleita", "Kuolinprosentti") #write.csv2(df1, "./kuolleet_ikaryhmittain.csv", sep = ";" )

write.csv(df1, "./thl_tapaukset_kuolleet.csv")

xy0<-data.frame(paiva, tapaukset) names(xy0)<-c("Dates", "Cases") xy<-na.omit(xy0)

#return(df1)

}

download_solanpaa_finnish_data<-function() { solanpaa_fi="https://covid19.solanpaa.fi/data/fin_cases.json" cache_file="solanpaa_fi.json"

download.file(solanpaa_fi, cache_file)

j1 <- fromJSON(cache_file)

 ## maybe errori

dates<-as.Date(j1$date)

dailycases<-j1$new_cases dailydeaths<-j1$new_deaths

   dataf1 <- data.frame("Date" = dates, "Paivitt_kuolemat"=dailydeaths) 
   dataf2 <- data.frame("Date" = dates, "Paivitt_tapaukset"=dailycases)

write.csv(dataf1, "./dailydeaths1.csv", row.names=T) write.csv(dataf2, "./dailycases1.csv", row.names=T)

xy0<-data.frame(dates, dailycases) names(xy0)<-c("Dates", "Cases") xy<-na.omit(xy0)

return(xy)

}

if(load_data_from==1) { xy<-load_data_from_finnish_wiki() print (xy) }

if(load_data_from==2) { xy<-load_data_from_aggregated() }

if(load_data_from==3) { xy<-download_solanpaa_finnish_data() }

if(load_data_from==4) { xy<-lataa_thl_tapaukset_kuolleet() }

names(xy)<-c("Dates","Cases")

#print (xy)

#print (beginday1)
select_datelimit_begin=as.Date(beginday1,format="%d/%m/%Y")
select_datelimit_end=as.Date(today1,format="%d/%m/%Y" )

#format(select_datelimit_begin, "%Y-%m-%d")

#print(select_datelimit_begin)

#print(select_datelimit_end)
#2020-12-16
 ##xy2<-xy[xy$Dates >= select_datelimit_begin,]

1. xy2<-xy[xy$Dates >= select_datelimit_begin & xy$Dates <= select_datelimit_end]

xy2a<-xy[xy$Dates >= select_datelimit_begin, ]

xy2<-xy2a[ xy2a$Dates <= select_datelimit_end,]

print("Sel")
print(today1)
print(select_datelimit_begin)
print(select_datelimit_end)
print(xy2)

1. 1. stop(-1)

cases1<-xy2$Cases
dates1<-xy2$Dates

xy3<-data.frame( as.Date(dates1),as.integer(cases1) )
names(xy3)<-c("Dates", "Cases")

len1=length(cases1)

dates2<-as.Date(dates1)
paivat<-1:len1

## test code
arrat0<-calculate_r0_with_simple_exponent_moving_average(xy2, 14,7)

arrat1<-calculate_r0_with_r0(xy2)
arrat2<-calculate_r0_with_epiestim(xy2)

 #plot(arrat$paivat, arrat$r0)

 arrat<-arrat2
 

#str(arrat)
#head(arrat)

 sarrat1<-arrat1
 sarrat2<-sarrat1
 names(sarrat1)<-c("Dates","Cases")
  

 datelimits2=c(today1, as.Date(forecastendday1,"%Y/%m/%d"))

 datelimits3=c(as.Date(beginday1, "%d/%m/%Y" ), as.Date(forecastendday1,"%Y/%m/%d"))

 daysek1<-seq(today, as.Date(forecastendday1), "days")
   forecastdaynum1<-length(daysek1)
 lendaysek1<-length(daysek1)
 daysek2<-seq(as.Date(beginday1),as.Date(forecastendday1), "days")
 daysek3<-seq(as.Date(beginday1),as.Date(forecastendday1), "days")

daysek4<-seq(as.Date(today), as.Date(forecastendday1,"%Y/%m/%d"),"days" )

# farrat1<-forecast_profet(xy2, lendaysek1)
 
 
# names(farrat1)<-c("Dates", "Cases")
 
# print(datelimits3)
 
 
 #print(farrat1)
 
  
 
 
# stop(-1)
 
 
 #plot(arrat$paivat, arrat$r0)
# plot(farrat1$paivat, farrat1$r0)
# lines(farrat1$paivat, farrat1$r0  )
# lines(sarrat2$paivat, sarrat1$r0  )

1. stop(-1)

#arrat<-farrat1



# names(arrat)<-c("Dates","Forecast")
 
# arrat$Forecast<-as.integer(arrat$Forecast)

#marrat <- left_join(farrat1, xy3, by=c("Dates"))

# names(marrat)<-c("Dates","Forecast","Cases")
 

1. print (marrat)

#stop(-1)

1. days2 <- seq(ISOdate(2020,11,22), by = "day", length.out = daynum)
2. days3 <- seq(ISOdate(2020,10,1), by = "day", length.out = daynum2)

y=xy3$Cases

y[y == 0] <- 1

yy=y

print (yy)

1. stop(-1)

logy=log(yy)

y=logy

print (length(y)) print (forecastdaynum1)

1. sarima_forecast = sarima.for(y, n.ahead=daynum,p=0,d=1,q=1,P=1,D=1,Q=0,S=7)

1. sarima_forecast = sarima.for(y, n.ahead=daynum,p=0,d=1,q=1,P=1,D=1,Q=0,S=7)

sarima_forecast = sarima.for(y, n.ahead= forecastdaynum1,p=0,d=1,q=1,P=1,D=1,Q=0, S=6)

str(sarima_forecast)

y1=sarima_forecast$pred y2=sarima_forecast$se

y1=exp(y1) y2=exp(y2)

days2<-daysek4

print (length(y2)) print (length(days2))

1. stop(-1)

1. print (length(days3))
2. print (length(y1))

ydelta1=y2*50 ydelta2=y2*25 ydelta3=y2*10 ydelta4=y2*5

ya1=y1+ydelta1 ya2=y1-ydelta1 ya3=y1+ydelta2 ya4=y1-ydelta2 ya5=y1+ydelta3 ya6=y1-ydelta3 ya7=y1+ydelta4 ya8=y1-ydelta4

lines(days2, ya1) lines(days2, ya2)

plot(days2, y1, type="l", main="Koronatapaukset ennuste",xlabel="Pvm", ylabel="Tapauksia" ) print (y2)

data<-data.frame(days2,y1,ya1,ya2,ya3,ya4)

if(tulosta_svg==1) {

svg(filename=tulosfilee1, width=12, height=6, pointsize=12)

}

if(plottaa==1) {

 metodi="loess"

1. spanni=0.1
2. metodi="loess"

ggplot(data, aes(x =days2 , y = y1) ) + ylim(yala, yyla) + geom_line()+

#ggtitle(" Ennuste - koronavirustapauksia Suomessa") + labs(title = "Ennuste: koronavirustapauksia Suomessa",

             subtitle = "jos muutos jatkuu samaa vauhtia",
             caption = "")+

xlab("Kuukausi") + ylab("Tapauksia")+ theme(title=element_text(size=15), axis.text=element_text(size=12,face="bold"),axis.title=element_text(size=14,face="bold"))+

#geom_smooth( fill="#a0a0ff",span=spanni, method=metodi, level=0.9999, size=2)+ geom_smooth( fill="#9090ff", span=spanni,method=metodi, level=0.7) + geom_smooth( fill="#8a08af", span=spanni, method=metodi,level=0.5) +

   geom_ribbon( aes(ymin=ya1,ymax=ya2), fill="blue", alpha=0.5) +
       
   geom_ribbon( aes(ymin=ya3,ymax=ya4), fill="blue", alpha=0.25)+
   geom_ribbon( aes(ymin=ya5,ymax=ya6), fill="blue", alpha=0.1) +
       
   geom_ribbon( aes(ymin=ya7,ymax=ya6), fill="blue", alpha=0.5)

}

if(tulosta_svg==1) { dev.off() }