1. 3. Calculate forecast of Covid-19 cases from r0 in Finland
   6. 31.8.2021
   7. v 0000.0000

new_in_skript=0

if (new_in_skript==1) { #install.packages("ggplot2", "plotly", repos ="https://ftp.acc.umu.se/mirror/CRAN/") install.packages("svglite") install.packages("ggplot2")

install.packages("rvest") install.packages("readtext") install.packages("stringi") install.packages("datamart") install.packages("XML")

install.packages("tidyr") install.packages("stringr") install.packages("stringi") install.packages("tibble") install.packages("readr") install.packages("data.table")

install.packages("caTools") install.packages("mgcv") install.packages("repmis") install.packages("lubridate") install.packages("tidyverse") install.packages("R0") install.packages("EpiEstim") install.packages("jsonlite") install.packages("rjstat") install.packages("dplyr") }

library(ggplot2) library(svglite)

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

library(tibble) library(caTools) library(mgcv) library(repmis) library(lubridate) library(tidyverse) library(dplyr) library(tidyr) library(readr) library(data.table) library(rjstat)

library (R0) library(EpiEstim)

library(forecast) library(prophet)

library(astsa) library(MLmetrics) library(tseries)

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

load_data_from=3

1. beginday1='28/02/2020'

1. beginday1='1/1/2021'

beginday1='1/8/2020'

forecastendday1<-"2021/10/16"

1. today=Sys.Date()-1

today=Sys.Date()-3

spanni=0.3

1. yala=0.6
2. yyla=3.2

yala=0.7 yyla=1.4

metodi="loess"

widthi=10 heighti=5

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

tulosfilee1="./R0_Suomessa_2.svg"

beginday0=as.Date(beginday1) beginday2=format(beginday0, "%Y/%m/%d")

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

print(today1)

datelimits1=c(beginday1, today1)

paivat1=seq(as.Date(beginday2), as.Date(today2), "days")

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)

}

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)

pv0<-tabu1$V1 len1=length(pv0) daates1 <- vector(mode="character", length=len1)

1. print(pv0)

n=1

for(n in 1:len1) { it1<-pv0[n] #print(it1) qq1<-str_split(it1, "\\[")1 qq2<-qq1[1] qq3<-gsub(" ", "", qq2, fixed = TRUE) daates1[n]=qq3 }

daates2=as.Date(daates1, format="%d.%m.%Y")

print(daates2)

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(daates2, sairaalassa)

xy<-data.frame(daates2, uusia_tapauksia) names(xy)<-c("Dates", "Cases")

xyz<-data.frame(daates2, sairaalassa, teholla) dfout1<-data.frame(daates2, 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

srkurl='https://datahub.io/core/covid-19/r/countries-aggregated.csv'

dfine000 <- read.csv(file=srkurl)

dfine<-as.data.frame(dfine000)

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

1. print(dfine)

1. tail(dfine)

1. stop(-1)

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

1. head(dfinland)

1. print(dfinland)

1. quit(-2)

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

pv0<-tabu1$V1 len1=length(pv0) daates1 <- vector(mode="character", length=len1)

1. print(pv0)

n=1

for(n in 1:len1) { it1<-pv0[n] #print(it1) qq1<-str_split(it1, "\\[")1 qq2<-qq1[1] qq3<-gsub(" ", "", qq2, fixed = TRUE) daates1[n]=qq3 }

daates2=as.Date(daates1, format="%d.%m.%Y")

print(daates2)

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

1. str(dataf1)

  dataf2 <- data.frame("Date" = daates2, "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(daates2, dailycases3)


return(xy)

}

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)

}

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)

}

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)

}

nth_element <- function(vector, starting_position, n) {

 vector[seq(starting_position, length(vector), n)] 
 }   

get_thl_hospital_data<-function() { url_base1="https://sampo.thl.fi/pivot/prod/fi/epirapo/covid19care/fact_epirapo_covid19care.json" request1="?row=dateweek20200101-508804L&column=measure-547523.547516.547531.456732.&fo=1"

url1 <- paste0(url_base1, request1)

cube1 <- fromJSONstat(url1, naming = "label", use_factors = F, silent = T)

res1 <- cube11

days0<-as.Date(res1$dateweek20200101) values1<-as.integer(res1$value)

dimx<-4 dimy<-length(values1)/dimx

values2<-values1 days1<-days0

dim(values2)<-c(dimx,dimy) dim(days1)<-c(dimx,dimy)

days2<-nth_element(days1, 1, 4)

esh1<-values2[1,] esh2<-values2[2,] ## esh esh3<-values2[3,] ## perusterv teho1<-values2[4,]

esh1[is.na(esh1)] <- 0 esh2[is.na(esh2)] <- 0 esh3[is.na(esh3)] <- 0 teho1[is.na(teho1)] <- 0

#print(esh1)

sairaalassa1<-esh1+esh2+esh3+teho1

#df1<-data.frame(days2, esh2, teho1) df1<-data.frame(days2, sairaalassa1, teho1)

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

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

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

}

1. 1. 1. main program

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() }

if(load_data_from==5) { xy<-get_thl_hospital_data() names(xy)<-c("Dates","Cases") }

print (xy)

1. quit(-1)

#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 & xy$Dates <= select_datelimit_end,]

#xy2<-xy[xy$Dates >= select_datelimit_begin,]

print("xy2")
print(xy2)

#stop(-1)

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

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)

#print("calcu ok")

 #plot(arrat$paivat, arrat$r0)

# arrat<-arrat2
arrat<-arrat0

 
str(arrat)
head(arrat)

 #plot(arrat$paivat, arrat$r0)
# stop(-1)

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" )

y=arrat$r0

print("R in")

y[is.na(y)] <- 1 y[!is.finite(y)] <- 1

basedatay=y basedatadays=arrat$paivat

print(basedatay)

1. stop(-1)

1. print (length(y))
2. print (forecastdaynum1)

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

print("YYYYYYYYYY")

1. print (y1)
2. print (y2)

1. plot(y1)

1. stop(-1)

print (basedatadays) print (length(basedatadays)) print (length(basedatay))

1. 1. as.Date(a$Number, origin = “1964-10-22”)

#basedata<-c( as.Date(basedatadays, origin="2021-06-01"), basedatay )
#basedata<-c( as.numeric(basedatadays), basedatay )
basedata<-data.frame(basedatadays, basedatay)

names(basedata)<-c("ds", "y")

 print(basedata)
 
# stop(-1)
 

1. farrat1<-forecast_profet(basedata, lendaysek1)

1. stop(-1)

 # profeta <- prophet(basedata,interval.width = 0.1, yearly.seasonality =0, weekly.seasonality = TRUE)
 #   profeta <- prophet(basedata,interval.width = 0.1, yearly.seasonality =0, daily.seasonality = TRUE)
#  profeta <- prophet(basedata,interval.width = 1.0, yearly.seasonality =1,weekly.seasonality = TRUE, daily.seasonality = TRUE)
  # profeta <- prophet(basedata, interval.width = 0.25, daily.seasonality = TRUE)
# profeta <- prophet(basedata, seasonality.mode = 'multiplicative')
profeta <- prophet(seasonality.mode = 'multiplicative')

#profeta <- add_seasonality(profeta, 'quarterly', period = 60, fourier.order = 8, mode = 'additive')
 

1. profeta <- add_seasonality(profeta, 'quarterly', period = 90, fourier.order = 8, mode = 'additive')

1. profeta <- add_seasonality(profeta, 'monthly', period = 60, fourier.order = 8, mode = 'additive')
2. profeta <- add_seasonality(profeta, 'weekly', period = 30, fourier.order = 8, mode = 'additive')

##  profeta <- add_seasonality(profeta, 'p0', period = 365, fourier.order = 8, mode = 'additive')
# profeta <- add_seasonality(profeta, 'p1', period = 180, fourier.order = 8, mode = 'additive')

1. profeta <- add_seasonality(profeta, 'p2', period = 90, fourier.order = 8, mode = 'additive')
   1. profeta <- add_seasonality(profeta, 'p3', period = 30, fourier.order = 8, mode = 'additive')

  profeta <- add_seasonality(profeta, 'p1', period = 120, fourier.order = 8, mode ='multiplicative')
profeta <- add_seasonality(profeta, 'p2', period = 60, fourier.order = 8, mode = 'multiplicative')
 profeta <- add_seasonality(profeta, 'p3', period = 30, fourier.order = 8, mode = 'multiplicative')

#profeta <- add_regressor(profeta, 'regressor', mode = 'additive')
profeta<-fit.prophet(profeta,basedata)

  future <- make_future_dataframe(profeta, periods = lendaysek1 )
  forecast1 <- predict(profeta, future)
  
  
 	svg(filename="prophet_r10.svg", width=widthi, height=heighti, pointsize=12) 
  
  plot(profeta, forecast1,  main="R0 Suomessa", xlab="pvm", ylab="R0",font.main=4, font.lab=4, font.sub=4, cex.main=1.25, cex.lab=3, cex.axis=3)
  dev.off() 
    
        str(forecast1)
      
   yhat1<-forecast1$yhat
   ds<-forecast1$ds
   
   rohve<-data.frame(ds, yhat1)
  
   

 #  stop(-1)
    
    
       

2. forecasti

#md = rwf(y,h=forecastlen1,drift=T,level=c(90,95),fan=FALSE,lambda=NULL) md = rwf(y,h=lendaysek1,drift=T,level=c(30,40),fan=FALSE,lambda=NULL)

str(md)

fitted1<-md$fitted mean1<-md$mean lower1<-md$lower[,1] upper1<-md$upper[,1]

#plot(fitted1) plot(mean1) lines(lower1)

lines(upper1)

#stop(-1)

plot(md)

dif_data <- diff(y)

arima1<-auto.arima(dif_data) foca1<-forecast(arima1)

str(foca1)

plot(foca1)

foca1_fitted1<-foca1$fitted foca1_mean1<-foca1$mean foca1_lower1<-foca1$lower[,1] foca1_upper1<-foca1$upper[,1]

plot(foca1_mean1) lines(foca1_lower1) lines(foca1_upper1)

#mdarima1_fitted1<-focal_fitted1 mdarima1_mean1<-mean1 mdarima1_lower1<-lower1 mdarima1_upper1<-upper1

#fit1 <- stl(y, s.window="periodic",t.window=length(y) ) #fit1 <- stl(y ) #arima <- forecast(fit1,h=lendaysek1,method ='arima')

#fit = arima(y, c(0, 1, 1), seasonal = list(order = c(0, 1, 1), period = 7))

#fit = arima(y, c(0, 1,1), seasonal = list(order = c(0, 1,1), period = 7))

fit = arima(y, c(0, 1,1), seasonal = list(order = c(0, 1,1), period = lendaysek1))

pred <- predict(fit, n.ahead = lendaysek1)

mdarima1_pred<-pred$pred mdarima1_se<-pred$pred #mdarima1_pred<-pred$se #mdarima_mean1<-arima$mean

bat1 <- tbats(y) fc2 <- forecast(bat1, h=lendaysek1) plot(fc2, ylab="Tapauksia")

sensor <- ts(y,frequency=7) # consider adding a start so you get nicer labelling on your chart. fit <- auto.arima(sensor) fcast <- forecast(fit, h=lendaysek1) plot(fcast) grid() print(" Fcast 1 ...")

str(fcast)

mdarima_pred<-fcast$mean

#stop(-1)

#plot(fc) str(pred)

#plot(pred)

print(mdarima1_pred)

#quit(-1)

# stop(-1)

1. r0=1.2

r0=arrat$r0

tau=5.5

log(r0)

k=log(r0)/tau

1. print(k)

print(arrat) print print("Koo ....")

1. sel1=as.Date(today1, format="%d/%m/%Y")

sel2=as.Date(forecastendday1,format="%Y/%m/%d") sel1=as.Date("2021/08/28", format="%Y/%m/%d")

print (sel1) print (sel2)

errat<-rohve

1. print (errat)

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

1. print (errat)

1. str(arrat)

1. errat5<-arrat[errat$paivat >= sel1 & errat$paivat <= sel2,]

errat5<-errat[errat$paivat >= sel1,]

1. arrat5 <- arrat[arrat$paivat > sel1, ]

1. arrat5<-subset(arrat, paivat> sel1 & paivat < sel2 )

1. arrat5<-filter(arrat, filter(paivat>sel1 & paivat<sel2))

1. arrat5 <- subset(arrat, paivat > "2021-08-28" & paivat < "2021-09-15")

print("ERRAT5")

errat5

ennustepaivat=errat5$paivat

1. 1. print (xy2)

print (str(xy2)) xytoday<-xy2[xy2$Dates == sel1,] print (xytoday) nykymaara=xytoday$Cases print (nykymaara) nykypaiva=xytoday$Dates

1. 1. print (length(errat5))

r0=errat5$r0 k=log(r0)/tau

maara=nykymaara edmaara=maara

print ("FORECSZT")

mamma<-c(nykymaara) pappa<-c(nykypaiva)

for (i in 1:length(k)) { kaija=exp(k[i]) maara=edmaara*kaija paiva=ennustepaivat[i] print(paiva) print (maara)

if(i>1) { print("*") mamma=append(mamma,maara) pappa=append(pappa,paiva) }

edmaara=maara }

print("Mamma") print (mamma)

mennuste<-data.frame(pappa, mamma) names(mennuste)<-c("paivat","maarat")

print (mennuste)

 	svg(filename="covid_cases_forecast_from_r0.svg", width=widthi, height=heighti, pointsize=12) 
	
	ggplot( mennuste, aes(x =paivat , y = maarat) ) +

ggtitle("Koronavirustapauksia - ennuste R0+Prophet") + xlab("Pvm") + ylab("Uusia tapauksia")+ theme(title=element_text(size=15), axis.text=element_text(size=12,face="bold"),axis.title=element_text(size=14,face="bold"))+ geom_point() + geom_smooth( fill="#a0a0ff",span=spanni, method=metodi, level=0.99, size=3)+ geom_smooth( fill="#9090ff", span=spanni,method=metodi, level=0.7) + geom_smooth( fill="#8a08af", span=spanni, method=metodi,level=0.5)

  dev.off() 
    

1. r0

1. k

stop(-1)

if(tulosta_svg==1) {

#svg(filename=tulosfilee1, width=8, height=3, pointsize=12) svg(filename=tulosfilee1, width=widthi, height=heighti, pointsize=12)

}

if(plottaa==1) {

 metodi="loess"

print ("ggplot")

#ggplot(arrat, aes(x =paivat , y = r0)) +ylim(0.6, 1.8)+xlim(as.Date(datelimits1, format="%d/%m/%Y") )+ ggplot(arrat, aes(x =paivat , y = r0)) +ylim(yala, yyla)+xlim(as.Date(datelimits1, format="%d/%m/%Y") )+ ggtitle("Arvioitu koronaviruksen perusuusiutumisluku R0") + xlab("Kuukausi") + ylab("R0")+ theme(title=element_text(size=15), axis.text=element_text(size=12,face="bold"),axis.title=element_text(size=14,face="bold"))+ geom_point() + geom_smooth( fill="#a0a0ff",span=spanni, method=metodi, level=0.99, size=3)+ geom_smooth( fill="#9090ff", span=spanni,method=metodi, level=0.7) + geom_smooth( fill="#8a08af", span=spanni, method=metodi,level=0.5) + geom_hline(yintercept=1.0, linetype="dashed", color = "red", size=1)

}

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