Epidemic.TA System for Forecasting Covid-19 Cases Using Time Series and Neural Networks Models

Cumulative Covid 19 Infection cases in Siberian Federal
Makarovskikh Tatyana Anatolyevna “Макаровских Татьяна Анатольевна”
Abotaleb mostafa“Аботалеб Мостафа”
Department of Electrical Engineering and Computer Science
South ural state university, Chelyabinsk, Russian federation
#Import
library(fpp2)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
## -- Attaching packages ---------------------------------------------- fpp2 2.4 --
## v ggplot2   3.3.2     v fma       2.4  
## v forecast  8.13      v expsmooth 2.3
## 
library(forecast)
library(ggplot2)
library("readxl")
library(moments)
library(forecast)
require(forecast)  
require(tseries)
## Loading required package: tseries
require(markovchain)
## Loading required package: markovchain
## Package:  markovchain
## Version:  0.8.5-3
## Date:     2020-12-03
## BugReport: https://github.com/spedygiorgio/markovchain/issues
require(data.table)
## Loading required package: data.table
library(Hmisc)
## Loading required package: lattice
## Loading required package: survival
## Loading required package: Formula
## 
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:base':
## 
##     format.pval, units
library(ascii)
library(pander)
## 
## Attaching package: 'pander'
## The following object is masked from 'package:ascii':
## 
##     Pandoc
##Global vriable##
Full_original_data <- read_excel("data2.xlsx", sheet = "Siberian _Federal") # path of your data ( time series data)
## New names:
## * region -> region...2
## * infection -> infection...3
## * `daily infection` -> `daily infection...4`
## * region -> region...5
## * infection -> infection...6
## * ...
original_data<-Full_original_data$`cumulative total`
y_lab <- "Cumulative Covid 19 Infection cases in Siberian Federal"   # input name of data
Actual_date_interval <- c("2020/03/12","2021/03/22")
Forecast_date_interval <- c("2021/03/23","2021/03/29")
validation_data_days <-8
Number_Neural<-5       # Number of Neural For model NNAR Model
NNAR_Model<- TRUE     #create new model (TRUE/FALSE)
frequency<-"days"
country.name <- "Siberian Federal"
# Data Preparation & calculate some of statistics measures
summary(original_data) # Summary your time series
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       0   30142  110674  163972  297772  446774
# calculate standard deviation 
data.frame(kurtosis=kurtosis(original_data))   # calculate Cofficient of kurtosis
##   kurtosis
## 1 1.896722
data.frame(skewness=skewness(original_data))  # calculate Cofficient of skewness
##    skewness
## 1 0.5821243
data.frame(Standard.deviation =sd(original_data))
##   Standard.deviation
## 1           148957.9
#processing on data (input data)
rows <- NROW(original_data) # calculate number of rows in time series (number of days)
training_data<-original_data[1:(rows-validation_data_days)] # Training data
testing_data<-original_data[(rows-validation_data_days+1):rows] #testing data
AD<-fulldate<-seq(as.Date(Actual_date_interval[1]),as.Date(Actual_date_interval[2]), frequency)  #input range for actual date
FD<-seq(as.Date(Forecast_date_interval[1]),as.Date(Forecast_date_interval[2]), frequency)  #input range forecasting date
N_forecasting_days<-nrow(data.frame(FD))  #calculate number of days that you want to forecasting
validation_dates<-tail(AD,validation_data_days) # select validation_dates
validation_data_by_name<-weekdays(validation_dates) # put names of validation dates
forecasting_data_by_name<-weekdays(FD)  # put names of Forecasting dates
#NNAR Model 
if(NNAR_Model==TRUE){
  data_series<-ts(training_data)
  model_NNAR<-nnetar(data_series, size = Number_Neural)
  saveRDS(model_NNAR, file = "model_NNAR.RDS")
  my_model <- readRDS("model_NNAR.RDS")
  accuracy(model_NNAR)  # accuracy on training data #Print Model Parameters
  model_NNAR
}
## Series: data_series 
## Model:  NNAR(1,5) 
## Call:   nnetar(y = data_series, size = Number_Neural)
## 
## Average of 20 networks, each of which is
## a 1-5-1 network with 16 weights
## options were - linear output units 
## 
## sigma^2 estimated as 49832
if(NNAR_Model==FALSE){
  data_series<-ts(training_data)
  #model_NNAR<-nnetar(data_series, size = Number_Numeral)
  model_NNAR <- readRDS("model_NNAR.RDS")
  accuracy(model_NNAR)  # accuracy on training data #Print Model Parameters
  model_NNAR
}

# Testing Data Evaluation
forecasting_NNAR <- forecast(model_NNAR, h=N_forecasting_days+validation_data_days)
validation_forecast<-head(forecasting_NNAR$mean,validation_data_days)
MAPE_Per_Day<-round(  abs(((testing_data-validation_forecast)/testing_data)*100)  ,3)
paste ("MAPE % For ",validation_data_days,frequency,"by using NNAR Model for  ==> ",y_lab, sep=" ")
## [1] "MAPE % For  8 days by using NNAR Model for  ==>  Cumulative Covid 19 Infection cases in Siberian Federal"
MAPE_Mean_All<-paste(round(mean(MAPE_Per_Day),3),"% MAPE ",validation_data_days,frequency,y_lab,sep=" ")
MAPE_Mean_All_NNAR<-round(mean(MAPE_Per_Day),3)
MAPE_NNAR<-paste(round(MAPE_Per_Day,3),"%")
MAPE_NNAR_Model<-paste(MAPE_Per_Day ,"%")
paste (" MAPE that's Error of Forecasting for ",validation_data_days," days in NNAR Model for  ==> ",y_lab, sep=" ")
## [1] " MAPE that's Error of Forecasting for  8  days in NNAR Model for  ==>  Cumulative Covid 19 Infection cases in Siberian Federal"
paste(MAPE_Mean_All,"%")
## [1] "0.193 % MAPE  8 days Cumulative Covid 19 Infection cases in Siberian Federal %"
paste ("MAPE that's Error of Forecasting day by day for ",validation_data_days," days in NNAR Model for  ==> ",y_lab, sep=" ")
## [1] "MAPE that's Error of Forecasting day by day for  8  days in NNAR Model for  ==>  Cumulative Covid 19 Infection cases in Siberian Federal"
print(ascii(data.frame(date_NNAR=validation_dates,validation_data_by_name,actual_data=testing_data,forecasting_NNAR=validation_forecast,MAPE_NNAR_Model)), type = "rest")
## 
## +---+------------+-------------------------+-------------+------------------+-----------------+
## |   | date_NNAR  | validation_data_by_name | actual_data | forecasting_NNAR | MAPE_NNAR_Model |
## +===+============+=========================+=============+==================+=================+
## | 1 | 2021-03-15 | Monday                  | 440602.00   | 440437.61        | 0.037 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 2 | 2021-03-16 | Tuesday                 | 441488.00   | 441166.91        | 0.073 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 3 | 2021-03-17 | Wednesday               | 442399.00   | 441877.19        | 0.118 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 4 | 2021-03-18 | Thursday                | 443296.00   | 442568.76        | 0.164 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 5 | 2021-03-19 | Friday                  | 444193.00   | 443241.90        | 0.214 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 6 | 2021-03-20 | Saturday                | 445078.00   | 443896.93        | 0.265 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 7 | 2021-03-21 | Sunday                  | 445932.00   | 444534.15        | 0.313 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 8 | 2021-03-22 | Monday                  | 446774.00   | 445153.88        | 0.363 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
print(ascii(data.frame(FD,forecating_date=forecasting_data_by_name,forecasting_by_NNAR=tail(forecasting_NNAR$mean,N_forecasting_days))), type = "rest")
## 
## +---+------------+-----------------+---------------------+
## |   | FD         | forecating_date | forecasting_by_NNAR |
## +===+============+=================+=====================+
## | 1 | 2021-03-23 | Tuesday         | 445756.45           |
## +---+------------+-----------------+---------------------+
## | 2 | 2021-03-24 | Wednesday       | 446342.16           |
## +---+------------+-----------------+---------------------+
## | 3 | 2021-03-25 | Thursday        | 446911.36           |
## +---+------------+-----------------+---------------------+
## | 4 | 2021-03-26 | Friday          | 447464.36           |
## +---+------------+-----------------+---------------------+
## | 5 | 2021-03-27 | Saturday        | 448001.50           |
## +---+------------+-----------------+---------------------+
## | 6 | 2021-03-28 | Sunday          | 448523.10           |
## +---+------------+-----------------+---------------------+
## | 7 | 2021-03-29 | Monday          | 449029.50           |
## +---+------------+-----------------+---------------------+
plot(forecasting_NNAR,xlab = paste ("Time in", frequency ,y_lab , sep=" "), ylab=y_lab)
x1_test <- ts(testing_data, start =(rows-validation_data_days+1) )
lines(x1_test, col='red',lwd=2)

graph1<-autoplot(forecasting_NNAR,xlab = paste ("Time in", frequency ,y_lab , sep=" "), ylab=y_lab)
graph1

##bats model
# Data Modeling
data_series<-ts(training_data) # make your data to time series
autoplot(data_series ,xlab=paste ("Time in", frequency, sep=" "), ylab = y_lab, main=paste ("Actual Data :", y_lab, sep=" "))

model_bats<-bats(data_series)
accuracy(model_bats)  # accuracy on training data
##                    ME     RMSE      MAE MPE MAPE       MASE        ACF1
## Training set 2.872936 48.47552 32.83925 NaN  Inf 0.02741029 0.005111758
# Print Model Parameters
model_bats
## BATS(1, {3,3}, 0.985, -)
## 
## Call: bats(y = data_series)
## 
## Parameters
##   Alpha: 0.8365729
##   Beta: 0.1800211
##   Damping Parameter: 0.985003
##   AR coefficients: 1.211127 -0.232627 0.013333
##   MA coefficients: -0.545915 -0.011527 -0.12166
## 
## Seed States:
##           [,1]
## [1,] -82.31101
## [2,]   5.69624
## [3,]   0.00000
## [4,]   0.00000
## [5,]   0.00000
## [6,]   0.00000
## [7,]   0.00000
## [8,]   0.00000
## 
## Sigma: 48.47552
## AIC: 5064.634
#ploting BATS Model
plot(model_bats,xlab = paste ("Time in", frequency ,y_lab , sep=" "))

# Testing Data Evaluation
forecasting_bats <- predict(model_bats, h=N_forecasting_days+validation_data_days)
validation_forecast<-head(forecasting_bats$mean,validation_data_days)
MAPE_Per_Day<-round(  abs(((testing_data-validation_forecast)/testing_data)*100)  ,3)
paste ("MAPE % For ",validation_data_days,frequency,"by using bats Model for  ==> ",y_lab, sep=" ")
## [1] "MAPE % For  8 days by using bats Model for  ==>  Cumulative Covid 19 Infection cases in Siberian Federal"
MAPE_Mean_All.bats_Model<-round(mean(MAPE_Per_Day),3)
MAPE_Mean_All.bats<-paste(round(mean(MAPE_Per_Day),3),"% MAPE ",validation_data_days,frequency,y_lab,sep=" ")
MAPE_bats<-paste(round(MAPE_Per_Day,3),"%")
MAPE_bats_Model<-paste(MAPE_Per_Day ,"%")
paste (" MAPE that's Error of Forecasting for ",validation_data_days," days in bats Model for  ==> ",y_lab, sep=" ")
## [1] " MAPE that's Error of Forecasting for  8  days in bats Model for  ==>  Cumulative Covid 19 Infection cases in Siberian Federal"
paste(MAPE_Mean_All.bats,"%")
## [1] "0.004 % MAPE  8 days Cumulative Covid 19 Infection cases in Siberian Federal %"
paste ("MAPE that's Error of Forecasting day by day for ",validation_data_days," days in bats Model for  ==> ",y_lab, sep=" ")
## [1] "MAPE that's Error of Forecasting day by day for  8  days in bats Model for  ==>  Cumulative Covid 19 Infection cases in Siberian Federal"
print(ascii(data.frame(date_bats=validation_dates,validation_data_by_name,actual_data=testing_data,forecasting_bats=validation_forecast,MAPE_bats_Model)), type = "rest")
## 
## +---+------------+-------------------------+-------------+------------------+-----------------+
## |   | date_bats  | validation_data_by_name | actual_data | forecasting_bats | MAPE_bats_Model |
## +===+============+=========================+=============+==================+=================+
## | 1 | 2021-03-15 | Monday                  | 440602.00   | 440608.23        | 0.001 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 2 | 2021-03-16 | Tuesday                 | 441488.00   | 441514.30        | 0.006 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 3 | 2021-03-17 | Wednesday               | 442399.00   | 442410.26        | 0.003 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 4 | 2021-03-18 | Thursday                | 443296.00   | 443296.68        | 0 %             |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 5 | 2021-03-19 | Friday                  | 444193.00   | 444173.63        | 0.004 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 6 | 2021-03-20 | Saturday                | 445078.00   | 445041.21        | 0.008 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 7 | 2021-03-21 | Sunday                  | 445932.00   | 445899.51        | 0.007 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 8 | 2021-03-22 | Monday                  | 446774.00   | 446748.63        | 0.006 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
print(ascii(data.frame(FD,forecating_date=forecasting_data_by_name,forecasting_by_bats=tail(forecasting_bats$mean,N_forecasting_days),lower=tail(forecasting_bats$lower,N_forecasting_days),Upper=tail(forecasting_bats$lower,N_forecasting_days))), type = "rest")
## 
## +---+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## |   | FD         | forecating_date | forecasting_by_bats | lower.80. | lower.95. | Upper.80. | Upper.95. |
## +===+============+=================+=====================+===========+===========+===========+===========+
## | 1 | 2021-03-23 | Tuesday         | 447588.68           | 446734.74 | 446282.70 | 446734.74 | 446282.70 |
## +---+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 2 | 2021-03-24 | Wednesday       | 448419.75           | 447412.68 | 446879.56 | 447412.68 | 446879.56 |
## +---+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 3 | 2021-03-25 | Thursday        | 449241.94           | 448068.85 | 447447.85 | 448068.85 | 447447.85 |
## +---+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 4 | 2021-03-26 | Friday          | 450055.35           | 448703.22 | 447987.44 | 448703.22 | 447987.44 |
## +---+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 5 | 2021-03-27 | Saturday        | 450860.06           | 449315.77 | 448498.27 | 449315.77 | 448498.27 |
## +---+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 6 | 2021-03-28 | Sunday          | 451656.18           | 449906.51 | 448980.29 | 449906.51 | 448980.29 |
## +---+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 7 | 2021-03-29 | Monday          | 452443.80           | 450475.45 | 449433.46 | 450475.45 | 449433.46 |
## +---+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
plot(forecasting_bats)
x1_test <- ts(testing_data, start =(rows-validation_data_days+1) )
lines(x1_test, col='red',lwd=2)

graph1<-autoplot(forecasting_bats,xlab = paste ("Time in", frequency ,y_lab , sep=" "), ylab=y_lab)
graph1

## TBATS Model
# Data Modeling
data_series<-ts(training_data)
model_TBATS<-forecast:::fitSpecificTBATS(data_series,use.box.cox=FALSE, use.beta=TRUE,  seasonal.periods=c(6),use.damping=FALSE,k.vector=c(2))
accuracy(model_TBATS)  # accuracy on training data
##                    ME     RMSE      MAE MPE MAPE       MASE        ACF1
## Training set 2.645299 51.29567 35.25475 NaN  Inf 0.02942646 0.004142224
# Print Model Parameters
model_TBATS
## TBATS(1, {0,0}, 1, {<6,2>})
## 
## Call: NULL
## 
## Parameters
##   Alpha: 0.8647007
##   Beta: 0.9528763
##   Damping Parameter: 1
##   Gamma-1 Values: -0.001926545
##   Gamma-2 Values: -4.569249e-05
## 
## Seed States:
##             [,1]
## [1,] -84.7403429
## [2,]   5.9287757
## [3,]   1.3172002
## [4,]  -0.2169835
## [5,]   0.3413989
## [6,]   1.8851393
## 
## Sigma: 51.29567
## AIC: 5092.253
plot(model_TBATS,xlab = paste ("Time in", frequency ,y_lab , sep=" "), ylab=y_lab)

# Testing Data Evaluation
forecasting_tbats <- predict(model_TBATS, h=N_forecasting_days+validation_data_days)
validation_forecast<-head(forecasting_tbats$mean,validation_data_days)
MAPE_Per_Day<-round(  abs(((testing_data-validation_forecast)/testing_data)*100)  ,3)
paste ("MAPE % For ",validation_data_days,frequency,"by using TBATS Model for  ==> ",y_lab, sep=" ")
## [1] "MAPE % For  8 days by using TBATS Model for  ==>  Cumulative Covid 19 Infection cases in Siberian Federal"
MAPE_Mean_All.TBATS_Model<-round(mean(MAPE_Per_Day),3)
MAPE_Mean_All.TBATS<-paste(round(mean(MAPE_Per_Day),3),"% MAPE ",validation_data_days,frequency,y_lab,sep=" ")
MAPE_TBATS<-paste(round(MAPE_Per_Day,3),"%")
MAPE_TBATS_Model<-paste(MAPE_Per_Day ,"%")
paste (" MAPE that's Error of Forecasting for ",validation_data_days," days in TBATS Model for  ==> ",y_lab, sep=" ")
## [1] " MAPE that's Error of Forecasting for  8  days in TBATS Model for  ==>  Cumulative Covid 19 Infection cases in Siberian Federal"
paste(MAPE_Mean_All.TBATS,"%")
## [1] "0.038 % MAPE  8 days Cumulative Covid 19 Infection cases in Siberian Federal %"
paste ("MAPE that's Error of Forecasting day by day for ",validation_data_days," days in TBATS Model for  ==> ",y_lab, sep=" ")
## [1] "MAPE that's Error of Forecasting day by day for  8  days in TBATS Model for  ==>  Cumulative Covid 19 Infection cases in Siberian Federal"
print(ascii(data.frame(date_TBATS=validation_dates,validation_data_by_name,actual_data=testing_data,forecasting_TBATS=validation_forecast,MAPE_TBATS_Model)), type = "rest")
## 
## +---+------------+-------------------------+-------------+-------------------+------------------+
## |   | date_TBATS | validation_data_by_name | actual_data | forecasting_TBATS | MAPE_TBATS_Model |
## +===+============+=========================+=============+===================+==================+
## | 1 | 2021-03-15 | Monday                  | 440602.00   | 440619.18         | 0.004 %          |
## +---+------------+-------------------------+-------------+-------------------+------------------+
## | 2 | 2021-03-16 | Tuesday                 | 441488.00   | 441554.70         | 0.015 %          |
## +---+------------+-------------------------+-------------+-------------------+------------------+
## | 3 | 2021-03-17 | Wednesday               | 442399.00   | 442489.38         | 0.02 %           |
## +---+------------+-------------------------+-------------+-------------------+------------------+
## | 4 | 2021-03-18 | Thursday                | 443296.00   | 443419.93         | 0.028 %          |
## +---+------------+-------------------------+-------------+-------------------+------------------+
## | 5 | 2021-03-19 | Friday                  | 444193.00   | 444357.79         | 0.037 %          |
## +---+------------+-------------------------+-------------+-------------------+------------------+
## | 6 | 2021-03-20 | Saturday                | 445078.00   | 445292.30         | 0.048 %          |
## +---+------------+-------------------------+-------------+-------------------+------------------+
## | 7 | 2021-03-21 | Sunday                  | 445932.00   | 446220.33         | 0.065 %          |
## +---+------------+-------------------------+-------------+-------------------+------------------+
## | 8 | 2021-03-22 | Monday                  | 446774.00   | 447155.85         | 0.085 %          |
## +---+------------+-------------------------+-------------+-------------------+------------------+
print(ascii(data.frame(FD,forecating_date=forecasting_data_by_name,forecasting_by_TBATS=tail(forecasting_tbats$mean,N_forecasting_days),Lower=tail(forecasting_tbats$lower,N_forecasting_days),Upper=tail(forecasting_tbats$upper,N_forecasting_days))), type = "rest")
## 
## +---+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## |   | FD         | forecating_date | forecasting_by_TBATS | Lower.80. | Lower.95. | Upper.80. | Upper.95. |
## +===+============+=================+======================+===========+===========+===========+===========+
## | 1 | 2021-03-23 | Tuesday         | 448090.53            | 447918.39 | 447827.26 | 448262.68 | 448353.81 |
## +---+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 2 | 2021-03-24 | Wednesday       | 449021.08            | 448840.13 | 448744.34 | 449202.04 | 449297.83 |
## +---+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 3 | 2021-03-25 | Thursday        | 449958.95            | 449769.62 | 449669.39 | 450148.28 | 450248.51 |
## +---+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 4 | 2021-03-26 | Friday          | 450893.46            | 450696.15 | 450591.71 | 451090.76 | 451195.21 |
## +---+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 5 | 2021-03-27 | Saturday        | 451821.49            | 451616.59 | 451508.12 | 452026.38 | 452134.85 |
## +---+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 6 | 2021-03-28 | Sunday          | 452757.00            | 452544.83 | 452432.52 | 452969.17 | 453081.49 |
## +---+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 7 | 2021-03-29 | Monday          | 453691.68            | 453472.50 | 453356.47 | 453910.87 | 454026.90 |
## +---+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
plot(forecasting_tbats)
x1_test <- ts(testing_data, start =(rows-validation_data_days+1) )
lines(x1_test, col='red',lwd=2)

graph2<-autoplot(forecasting_tbats,xlab = paste ("Time in", frequency ,y_lab , sep=" "), ylab=y_lab)
graph2

## Holt's linear trend
# Data Modeling
data_series<-ts(training_data)
model_holt<-holt(data_series,h=N_forecasting_days+validation_data_days,lambda = "auto")
accuracy(model_holt)  # accuracy on training data
##                     ME     RMSE      MAE MPE MAPE       MASE      ACF1
## Training set -10.36445 60.39636 43.58967 NaN  Inf 0.03638346 0.4907488
# Print Model Parameters
summary(model_holt$model)
## Holt's method 
## 
## Call:
##  holt(y = data_series, h = N_forecasting_days + validation_data_days,  
## 
##  Call:
##      lambda = "auto") 
## 
##   Box-Cox transformation: lambda= 0.5277 
## 
##   Smoothing parameters:
##     alpha = 0.9364 
##     beta  = 0.4524 
## 
##   Initial states:
##     l = -2.8431 
##     b = 0.8594 
## 
##   sigma:  0.6259
## 
##      AIC     AICc      BIC 
## 1835.240 1835.406 1854.781 
## 
## Training set error measures:
##                     ME     RMSE      MAE MPE MAPE       MASE      ACF1
## Training set -10.36445 60.39636 43.58967 NaN  Inf 0.03638346 0.4907488
# Testing Data Evaluation
forecasting_holt <- predict(model_holt, h=N_forecasting_days+validation_data_days,lambda = "auto")
validation_forecast<-head(forecasting_holt$mean,validation_data_days)
MAPE_Per_Day<-round(  abs(((testing_data-validation_forecast)/testing_data)*100)  ,3)
paste ("MAPE % For ",validation_data_days,frequency,"by using holt Model for  ==> ",y_lab, sep=" ")
## [1] "MAPE % For  8 days by using holt Model for  ==>  Cumulative Covid 19 Infection cases in Siberian Federal"
MAPE_Mean_All.Holt_Model<-round(mean(MAPE_Per_Day),3)
MAPE_Mean_All.Holt<-paste(round(mean(MAPE_Per_Day),3),"% MAPE ",validation_data_days,frequency,y_lab,sep=" ")
MAPE_holt<-paste(round(MAPE_Per_Day,3),"%")
MAPE_holt_Model<-paste(MAPE_Per_Day ,"%")
paste (" MAPE that's Error of Forecasting for ",validation_data_days," days in holt Model for  ==> ",y_lab, sep=" ")
## [1] " MAPE that's Error of Forecasting for  8  days in holt Model for  ==>  Cumulative Covid 19 Infection cases in Siberian Federal"
paste(MAPE_Mean_All.Holt,"%")
## [1] "0.049 % MAPE  8 days Cumulative Covid 19 Infection cases in Siberian Federal %"
paste ("MAPE that's Error of Forecasting day by day for ",validation_data_days," days in holt Model for  ==> ",y_lab, sep=" ")
## [1] "MAPE that's Error of Forecasting day by day for  8  days in holt Model for  ==>  Cumulative Covid 19 Infection cases in Siberian Federal"
print(ascii(data.frame(date_holt=validation_dates,validation_data_by_name,actual_data=testing_data,forecasting_holt=validation_forecast,MAPE_holt_Model)), type = "rest")
## 
## +---+------------+-------------------------+-------------+------------------+-----------------+
## |   | date_holt  | validation_data_by_name | actual_data | forecasting_holt | MAPE_holt_Model |
## +===+============+=========================+=============+==================+=================+
## | 1 | 2021-03-15 | Monday                  | 440602.00   | 440631.48        | 0.007 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 2 | 2021-03-16 | Tuesday                 | 441488.00   | 441573.98        | 0.019 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 3 | 2021-03-17 | Wednesday               | 442399.00   | 442517.44        | 0.027 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 4 | 2021-03-18 | Thursday                | 443296.00   | 443461.85        | 0.037 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 5 | 2021-03-19 | Friday                  | 444193.00   | 444407.21        | 0.048 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 6 | 2021-03-20 | Saturday                | 445078.00   | 445353.52        | 0.062 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 7 | 2021-03-21 | Sunday                  | 445932.00   | 446300.78        | 0.083 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 8 | 2021-03-22 | Monday                  | 446774.00   | 447248.99        | 0.106 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
print(ascii(data.frame(FD,forecating_date=forecasting_data_by_name,forecasting_by_holt=tail(forecasting_holt$mean,N_forecasting_days),Lower=tail(forecasting_holt$lower,N_forecasting_days),Upper=tail(forecasting_holt$upper,N_forecasting_days))), type = "rest")
## 
## +---+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## |   | FD         | forecating_date | forecasting_by_holt | Lower.80. | Lower.95. | Upper.80. | Upper.95. |
## +===+============+=================+=====================+===========+===========+===========+===========+
## | 1 | 2021-03-23 | Tuesday         | 448198.15           | 444848.09 | 443079.48 | 451560.08 | 453344.58 |
## +---+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 2 | 2021-03-24 | Wednesday       | 449148.26           | 445307.49 | 443280.62 | 453004.60 | 455052.33 |
## +---+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 3 | 2021-03-25 | Thursday        | 450099.32           | 445745.68 | 443449.09 | 454472.95 | 456796.30 |
## +---+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 4 | 2021-03-26 | Friday          | 451051.34           | 446163.46 | 443586.16 | 455964.35 | 458575.32 |
## +---+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 5 | 2021-03-27 | Saturday        | 452004.30           | 446561.59 | 443693.00 | 457478.14 | 460388.41 |
## +---+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 6 | 2021-03-28 | Sunday          | 452958.21           | 446940.72 | 443770.64 | 459013.70 | 462234.65 |
## +---+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 7 | 2021-03-29 | Monday          | 453913.07           | 447301.44 | 443820.01 | 460570.50 | 464113.26 |
## +---+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
plot(forecasting_holt)
x1_test <- ts(testing_data, start =(rows-validation_data_days+1) )
lines(x1_test, col='red',lwd=2)

graph3<-autoplot(forecasting_holt,xlab = paste ("Time in", frequency ,y_lab , sep=" "),  ylab=y_lab)
graph3

#Auto arima model
##################
require(tseries) # need to install tseries tj test Stationarity in time series 
paste ("tests For Check Stationarity in series  ==> ",y_lab, sep=" ")
## [1] "tests For Check Stationarity in series  ==>  Cumulative Covid 19 Infection cases in Siberian Federal"
kpss.test(data_series) # applay kpss test
## Warning in kpss.test(data_series): p-value smaller than printed p-value
## 
##  KPSS Test for Level Stationarity
## 
## data:  data_series
## KPSS Level = 5.912, Truncation lag parameter = 5, p-value = 0.01
pp.test(data_series)   # applay pp test
## 
##  Phillips-Perron Unit Root Test
## 
## data:  data_series
## Dickey-Fuller Z(alpha) = -1.985, Truncation lag parameter = 5, p-value
## = 0.9691
## alternative hypothesis: stationary
adf.test(data_series)  # applay adf test
## 
##  Augmented Dickey-Fuller Test
## 
## data:  data_series
## Dickey-Fuller = -3.3457, Lag order = 7, p-value = 0.06348
## alternative hypothesis: stationary
ndiffs(data_series)    # Doing first diffrencing on data
## [1] 2
#Taking the first difference
diff1_x1<-diff(data_series)
autoplot(diff1_x1, xlab = paste ("Time in", frequency ,y_lab , sep=" "), ylab=y_lab,main = "1nd differenced series")

##Testing the stationary of the first differenced series
paste ("tests For Check Stationarity in series after taking first differences in  ==> ",y_lab, sep=" ")
## [1] "tests For Check Stationarity in series after taking first differences in  ==>  Cumulative Covid 19 Infection cases in Siberian Federal"
kpss.test(diff1_x1)   # applay kpss test after taking first differences
## Warning in kpss.test(diff1_x1): p-value smaller than printed p-value
## 
##  KPSS Test for Level Stationarity
## 
## data:  diff1_x1
## KPSS Level = 3.9886, Truncation lag parameter = 5, p-value = 0.01
pp.test(diff1_x1)     # applay pp test after taking first differences
## Warning in pp.test(diff1_x1): p-value greater than printed p-value
## 
##  Phillips-Perron Unit Root Test
## 
## data:  diff1_x1
## Dickey-Fuller Z(alpha) = 1.3463, Truncation lag parameter = 5, p-value
## = 0.99
## alternative hypothesis: stationary
adf.test(diff1_x1)    # applay adf test after taking first differences
## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff1_x1
## Dickey-Fuller = -0.42186, Lag order = 7, p-value = 0.9851
## alternative hypothesis: stationary
#Taking the second difference
diff2_x1=diff(diff1_x1)
autoplot(diff2_x1, xlab = paste ("Time in", frequency ,y_lab , sep=" "), ylab=y_lab ,main = "2nd differenced series")

##Testing the stationary of the first differenced series
paste ("tests For Check Stationarity in series after taking Second differences in",y_lab, sep=" ")
## [1] "tests For Check Stationarity in series after taking Second differences in Cumulative Covid 19 Infection cases in Siberian Federal"
kpss.test(diff2_x1)   # applay kpss test after taking Second differences
## Warning in kpss.test(diff2_x1): p-value smaller than printed p-value
## 
##  KPSS Test for Level Stationarity
## 
## data:  diff2_x1
## KPSS Level = 0.92814, Truncation lag parameter = 5, p-value = 0.01
pp.test(diff2_x1)     # applay pp test after taking Second differences
## Warning in pp.test(diff2_x1): p-value smaller than printed p-value
## 
##  Phillips-Perron Unit Root Test
## 
## data:  diff2_x1
## Dickey-Fuller Z(alpha) = -489.65, Truncation lag parameter = 5, p-value
## = 0.01
## alternative hypothesis: stationary
adf.test(diff2_x1)    # applay adf test after taking Second differences
## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff2_x1
## Dickey-Fuller = -3.8363, Lag order = 7, p-value = 0.01745
## alternative hypothesis: stationary
####Fitting an ARIMA Model
#1. Using auto arima function
model1 <- auto.arima(data_series,stepwise=FALSE, approximation=FALSE, trace=T, test = c("kpss", "adf", "pp"))  #applaying auto arima
## 
##  ARIMA(0,2,0)                    : 3943.233
##  ARIMA(0,2,1)                    : 3935.131
##  ARIMA(0,2,2)                    : 3929.594
##  ARIMA(0,2,3)                    : 3930.803
##  ARIMA(0,2,4)                    : 3932.742
##  ARIMA(0,2,5)                    : 3912.492
##  ARIMA(1,2,0)                    : 3932.868
##  ARIMA(1,2,1)                    : 3934.146
##  ARIMA(1,2,2)                    : 3892.752
##  ARIMA(1,2,3)                    : Inf
##  ARIMA(1,2,4)                    : 3930.628
##  ARIMA(2,2,0)                    : 3932.783
##  ARIMA(2,2,1)                    : 3899.371
##  ARIMA(2,2,2)                    : 3926.514
##  ARIMA(2,2,3)                    : Inf
##  ARIMA(3,2,0)                    : 3925.848
##  ARIMA(3,2,1)                    : 3925.008
##  ARIMA(3,2,2)                    : 3896.302
##  ARIMA(4,2,0)                    : 3927.486
##  ARIMA(4,2,1)                    : 3927.075
##  ARIMA(5,2,0)                    : 3913.38
## 
## 
## 
##  Best model: ARIMA(1,2,2)
model1 # show the result of autoarima 
## Series: data_series 
## ARIMA(1,2,2) 
## 
## Coefficients:
##         ar1      ma1     ma2
##       0.978  -1.2887  0.3633
## s.e.  0.013   0.0567  0.0547
## 
## sigma^2 estimated as 2397:  log likelihood=-1942.32
## AIC=3892.64   AICc=3892.75   BIC=3908.25
#Make changes in the source of auto arima to run the best model
arima.string <- function (object, padding = FALSE) 
{
  order <- object$arma[c(1, 6, 2, 3, 7, 4, 5)]
  m <- order[7]
  result <- paste("ARIMA(", order[1], ",", order[2], ",", 
                  order[3], ")", sep = "")
  if (m > 1 && sum(order[4:6]) > 0) {
    result <- paste(result, "(", order[4], ",", order[5], 
                    ",", order[6], ")[", m, "]", sep = "")
  }
  if (padding && m > 1 && sum(order[4:6]) == 0) {
    result <- paste(result, "         ", sep = "")
    if (m <= 9) {
      result <- paste(result, " ", sep = "")
    }
    else if (m <= 99) {
      result <- paste(result, "  ", sep = "")
    }
    else {
      result <- paste(result, "   ", sep = "")
    }
  }
  if (!is.null(object$xreg)) {
    if (NCOL(object$xreg) == 1 && is.element("drift", names(object$coef))) {
      result <- paste(result, "with drift        ")
    }
    else {
      result <- paste("Regression with", result, "errors")
    }
  }
  else {
    if (is.element("constant", names(object$coef)) || is.element("intercept", 
                                                                 names(object$coef))) {
      result <- paste(result, "with non-zero mean")
    }
    else if (order[2] == 0 && order[5] == 0) {
      result <- paste(result, "with zero mean    ")
    }
    else {
      result <- paste(result, "                  ")
    }
  }
  if (!padding) {
    result <- gsub("[ ]*$", "", result)
  }
  return(result)
}


bestmodel <- arima.string(model1, padding = TRUE)
bestmodel <- substring(bestmodel,7,11)
bestmodel <- gsub(" ", "", bestmodel)
bestmodel <- gsub(")", "", bestmodel)
bestmodel <- strsplit(bestmodel, ",")[[1]]
bestmodel <- c(strtoi(bestmodel[1]),strtoi(bestmodel[2]),strtoi(bestmodel[3]))
bestmodel
## [1] 1 2 2
strtoi(bestmodel[3])
## [1] 2
#2. Using ACF and PACF Function
#par(mfrow=c(1,2))  # Code for making two plot in one graph 
acf(diff2_x1,xlab = paste ("Time in", frequency ,y_lab , sep=" ") , ylab=y_lab, main=paste("ACF-2nd differenced series ",y_lab, sep=" ",lag.max=20))    # plot ACF "auto correlation function after taking second diffrences

pacf(diff2_x1,xlab = paste ("Time in", frequency ,y_lab , sep=" "), ylab=y_lab,main=paste("PACF-2nd differenced series ",y_lab, sep=" ",lag.max=20))   # plot PACF " Partial auto correlation function after taking second diffrences

library(forecast)   # install library forecast             
x1_model1= arima(data_series, order=c(bestmodel)) # Run Best model of auto arima  for forecasting
x1_model1  # Show result of best model of auto arima 
## 
## Call:
## arima(x = data_series, order = c(bestmodel))
## 
## Coefficients:
##         ar1      ma1     ma2
##       0.978  -1.2887  0.3633
## s.e.  0.013   0.0567  0.0547
## 
## sigma^2 estimated as 2377:  log likelihood = -1942.32,  aic = 3892.64
paste ("accuracy of autoarima Model For  ==> ",y_lab, sep=" ")
## [1] "accuracy of autoarima Model For  ==>  Cumulative Covid 19 Infection cases in Siberian Federal"
accuracy(x1_model1)  # aacuracy of best model from auto arima
##                     ME     RMSE      MAE       MPE     MAPE       MASE
## Training set 0.4105546 48.62617 32.85027 0.2766923 1.504977 0.02741949
##                    ACF1
## Training set 0.02186529
x1_model1$x          # show result of best model from auto arima 
## NULL
checkresiduals(x1_model1,xlab = paste ("Time in", frequency ,y_lab , sep=" "), ylab=y_lab)  # checkresiduals from best model from using auto arima 

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,2,2)
## Q* = 23.794, df = 7, p-value = 0.001238
## 
## Model df: 3.   Total lags used: 10
paste("Box-Ljung test , Ljung-Box test For Modelling for   ==> ",y_lab, sep=" ")
## [1] "Box-Ljung test , Ljung-Box test For Modelling for   ==>  Cumulative Covid 19 Infection cases in Siberian Federal"
Box.test(x1_model1$residuals^2, lag=20, type="Ljung-Box")   # Do test for resdulas by using Box-Ljung test , Ljung-Box test For Modelling
## 
##  Box-Ljung test
## 
## data:  x1_model1$residuals^2
## X-squared = 35.591, df = 20, p-value = 0.01717
library(tseries)
jarque.bera.test(x1_model1$residuals)  # Do test jarque.bera.test 
## 
##  Jarque Bera Test
## 
## data:  x1_model1$residuals
## X-squared = 2025.4, df = 2, p-value < 2.2e-16
#Actual Vs Fitted
plot(data_series, col='red',lwd=2, main="Actual vs Fitted Plot", xlab='Time in (days)', ylab=y_lab) # plot actual and Fitted model 
lines(fitted(x1_model1), col='black')

#Test data
x1_test <- ts(testing_data, start =(rows-validation_data_days+1) ) # make testing data in time series and start from rows-6
forecasting_auto_arima <- forecast(x1_model1, h=N_forecasting_days+validation_data_days)
validation_forecast<-head(forecasting_auto_arima$mean,validation_data_days)
MAPE_Per_Day<-round(abs(((testing_data-validation_forecast)/testing_data)*100)  ,3)
paste ("MAPE % For ",validation_data_days,frequency,"by using bats Model for  ==> ",y_lab, sep=" ")
## [1] "MAPE % For  8 days by using bats Model for  ==>  Cumulative Covid 19 Infection cases in Siberian Federal"
MAPE_Mean_All.ARIMA_Model<-round(mean(MAPE_Per_Day),3)
MAPE_Mean_All.ARIMA<-paste(round(mean(MAPE_Per_Day),3),"% MAPE ",validation_data_days,frequency,y_lab,sep=" ")
MAPE_auto_arima<-paste(round(MAPE_Per_Day,3),"%")
MAPE_auto.arima_Model<-paste(MAPE_Per_Day ,"%")
paste (" MAPE that's Error of Forecasting for ",validation_data_days," days in bats Model for  ==> ",y_lab, sep=" ")
## [1] " MAPE that's Error of Forecasting for  8  days in bats Model for  ==>  Cumulative Covid 19 Infection cases in Siberian Federal"
paste(MAPE_Mean_All.ARIMA,"%")
## [1] "0.006 % MAPE  8 days Cumulative Covid 19 Infection cases in Siberian Federal %"
paste ("MAPE that's Error of Forecasting day by day for ",validation_data_days," days in bats Model for  ==> ",y_lab, sep=" ")
## [1] "MAPE that's Error of Forecasting day by day for  8  days in bats Model for  ==>  Cumulative Covid 19 Infection cases in Siberian Federal"
print(ascii(data.frame(date_auto.arima=validation_dates,validation_data_by_name,actual_data=testing_data,forecasting_auto.arima=validation_forecast,MAPE_auto.arima_Model)), type = "rest")
## 
## +---+-----------------+-------------------------+-------------+------------------------+-----------------------+
## |   | date_auto.arima | validation_data_by_name | actual_data | forecasting_auto.arima | MAPE_auto.arima_Model |
## +===+=================+=========================+=============+========================+=======================+
## | 1 | 2021-03-15      | Monday                  | 440602.00   | 440610.99              | 0.002 %               |
## +---+-----------------+-------------------------+-------------+------------------------+-----------------------+
## | 2 | 2021-03-16      | Tuesday                 | 441488.00   | 441524.05              | 0.008 %               |
## +---+-----------------+-------------------------+-------------+------------------------+-----------------------+
## | 3 | 2021-03-17      | Wednesday               | 442399.00   | 442428.38              | 0.007 %               |
## +---+-----------------+-------------------------+-------------+------------------------+-----------------------+
## | 4 | 2021-03-18      | Thursday                | 443296.00   | 443324.17              | 0.006 %               |
## +---+-----------------+-------------------------+-------------+------------------------+-----------------------+
## | 5 | 2021-03-19      | Friday                  | 444193.00   | 444211.60              | 0.004 %               |
## +---+-----------------+-------------------------+-------------+------------------------+-----------------------+
## | 6 | 2021-03-20      | Saturday                | 445078.00   | 445090.87              | 0.003 %               |
## +---+-----------------+-------------------------+-------------+------------------------+-----------------------+
## | 7 | 2021-03-21      | Sunday                  | 445932.00   | 445962.15              | 0.007 %               |
## +---+-----------------+-------------------------+-------------+------------------------+-----------------------+
## | 8 | 2021-03-22      | Monday                  | 446774.00   | 446825.61              | 0.012 %               |
## +---+-----------------+-------------------------+-------------+------------------------+-----------------------+
print(ascii(data.frame(FD,forecating_date=forecasting_data_by_name,forecasting_by_auto.arima=tail(forecasting_auto_arima$mean,N_forecasting_days),Lower=tail(forecasting_auto_arima$lower,N_forecasting_days),Upper=tail(forecasting_auto_arima$upper,N_forecasting_days))), type = "rest")
## 
## +---+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## |   | FD         | forecating_date | forecasting_by_auto.arima | Lower.80. | Lower.95. | Upper.80. | Upper.95. |
## +===+============+=================+===========================+===========+===========+===========+===========+
## | 1 | 2021-03-23 | Tuesday         | 447681.43                 | 446768.08 | 446284.57 | 448594.79 | 449078.30 |
## +---+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 2 | 2021-03-24 | Wednesday       | 448529.78                 | 447449.31 | 446877.35 | 449610.26 | 450182.22 |
## +---+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 3 | 2021-03-25 | Thursday        | 449370.82                 | 448109.42 | 447441.67 | 450632.23 | 451299.98 |
## +---+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 4 | 2021-03-26 | Friday          | 450204.72                 | 448748.46 | 447977.57 | 451660.97 | 452431.87 |
## +---+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 5 | 2021-03-27 | Saturday        | 451031.62                 | 449366.51 | 448485.05 | 452696.73 | 453578.19 |
## +---+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 6 | 2021-03-28 | Sunday          | 451851.69                 | 449963.62 | 448964.14 | 453739.75 | 454739.23 |
## +---+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 7 | 2021-03-29 | Monday          | 452665.06                 | 450539.87 | 449414.86 | 454790.26 | 455915.27 |
## +---+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
plot(forecasting_auto_arima)
x1_test <- ts(testing_data, start =(rows-validation_data_days+1) )
lines(x1_test, col='red',lwd=2)

graph4<-autoplot(forecasting_auto_arima,xlab = paste ("Time in", frequency ,y_lab , sep=" "), ylab=y_lab)
graph4

MAPE_Mean_All.ARIMA
## [1] "0.006 % MAPE  8 days Cumulative Covid 19 Infection cases in Siberian Federal"
# Table for MAPE For counry
best_recommended_model <- min(MAPE_Mean_All_NNAR,MAPE_Mean_All.bats_Model,MAPE_Mean_All.TBATS_Model,MAPE_Mean_All.Holt_Model,MAPE_Mean_All.ARIMA_Model)
paste("System Choose Least Error ==> ( MAPE %) of Forecasting  by using bats model and BATS Model, Holt's Linear Models , and autoarima for  ==> ", y_lab , sep=" ")
## [1] "System Choose Least Error ==> ( MAPE %) of Forecasting  by using bats model and BATS Model, Holt's Linear Models , and autoarima for  ==>  Cumulative Covid 19 Infection cases in Siberian Federal"
best_recommended_model
## [1] 0.004
x1<-if(best_recommended_model >= MAPE_Mean_All.bats_Model) {paste("BATS Model")}
x2<-if(best_recommended_model >= MAPE_Mean_All.TBATS_Model) {paste("TBATS Model")}
x3<-if(best_recommended_model >= MAPE_Mean_All.Holt_Model) {paste("Holt Model")}
x4<-if(best_recommended_model >= MAPE_Mean_All.ARIMA_Model) {paste("ARIMA Model")}
x5<-if(best_recommended_model >= MAPE_Mean_All_NNAR) {paste("NNAR Model")}
panderOptions('table.split.table', Inf)
paste("Forecasting by using BATS Model  ==> ", y_lab , sep=" ")
## [1] "Forecasting by using BATS Model  ==>  Cumulative Covid 19 Infection cases in Siberian Federal"
print(ascii(data.frame(FD,forecating_date=forecasting_data_by_name,forecasting_by_bats=tail(forecasting_bats$mean,N_forecasting_days),lower=tail(forecasting_bats$lower,N_forecasting_days),Upper=tail(forecasting_bats$lower,N_forecasting_days))), type = "rest")
## 
## +---+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## |   | FD         | forecating_date | forecasting_by_bats | lower.80. | lower.95. | Upper.80. | Upper.95. |
## +===+============+=================+=====================+===========+===========+===========+===========+
## | 1 | 2021-03-23 | Tuesday         | 447588.68           | 446734.74 | 446282.70 | 446734.74 | 446282.70 |
## +---+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 2 | 2021-03-24 | Wednesday       | 448419.75           | 447412.68 | 446879.56 | 447412.68 | 446879.56 |
## +---+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 3 | 2021-03-25 | Thursday        | 449241.94           | 448068.85 | 447447.85 | 448068.85 | 447447.85 |
## +---+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 4 | 2021-03-26 | Friday          | 450055.35           | 448703.22 | 447987.44 | 448703.22 | 447987.44 |
## +---+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 5 | 2021-03-27 | Saturday        | 450860.06           | 449315.77 | 448498.27 | 449315.77 | 448498.27 |
## +---+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 6 | 2021-03-28 | Sunday          | 451656.18           | 449906.51 | 448980.29 | 449906.51 | 448980.29 |
## +---+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 7 | 2021-03-29 | Monday          | 452443.80           | 450475.45 | 449433.46 | 450475.45 | 449433.46 |
## +---+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
paste("Forecasting by using TBATS Model  ==> ", y_lab , sep=" ")
## [1] "Forecasting by using TBATS Model  ==>  Cumulative Covid 19 Infection cases in Siberian Federal"
print(ascii(data.frame(FD,forecating_date=forecasting_data_by_name,forecasting_by_TBATS=tail(forecasting_tbats$mean,N_forecasting_days),Lower=tail(forecasting_tbats$lower,N_forecasting_days),Upper=tail(forecasting_tbats$upper,N_forecasting_days))), type = "rest")
## 
## +---+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## |   | FD         | forecating_date | forecasting_by_TBATS | Lower.80. | Lower.95. | Upper.80. | Upper.95. |
## +===+============+=================+======================+===========+===========+===========+===========+
## | 1 | 2021-03-23 | Tuesday         | 448090.53            | 447918.39 | 447827.26 | 448262.68 | 448353.81 |
## +---+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 2 | 2021-03-24 | Wednesday       | 449021.08            | 448840.13 | 448744.34 | 449202.04 | 449297.83 |
## +---+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 3 | 2021-03-25 | Thursday        | 449958.95            | 449769.62 | 449669.39 | 450148.28 | 450248.51 |
## +---+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 4 | 2021-03-26 | Friday          | 450893.46            | 450696.15 | 450591.71 | 451090.76 | 451195.21 |
## +---+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 5 | 2021-03-27 | Saturday        | 451821.49            | 451616.59 | 451508.12 | 452026.38 | 452134.85 |
## +---+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 6 | 2021-03-28 | Sunday          | 452757.00            | 452544.83 | 452432.52 | 452969.17 | 453081.49 |
## +---+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 7 | 2021-03-29 | Monday          | 453691.68            | 453472.50 | 453356.47 | 453910.87 | 454026.90 |
## +---+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
paste("Forecasting by using Holt's Linear Trend Model  ==> ", y_lab , sep=" ")
## [1] "Forecasting by using Holt's Linear Trend Model  ==>  Cumulative Covid 19 Infection cases in Siberian Federal"
print(ascii(data.frame(FD,forecating_date=forecasting_data_by_name,forecasting_by_holt=tail(forecasting_holt$mean,N_forecasting_days),Lower=tail(forecasting_holt$lower,N_forecasting_days),Upper=tail(forecasting_holt$upper,N_forecasting_days))), type = "rest")
## 
## +---+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## |   | FD         | forecating_date | forecasting_by_holt | Lower.80. | Lower.95. | Upper.80. | Upper.95. |
## +===+============+=================+=====================+===========+===========+===========+===========+
## | 1 | 2021-03-23 | Tuesday         | 448198.15           | 444848.09 | 443079.48 | 451560.08 | 453344.58 |
## +---+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 2 | 2021-03-24 | Wednesday       | 449148.26           | 445307.49 | 443280.62 | 453004.60 | 455052.33 |
## +---+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 3 | 2021-03-25 | Thursday        | 450099.32           | 445745.68 | 443449.09 | 454472.95 | 456796.30 |
## +---+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 4 | 2021-03-26 | Friday          | 451051.34           | 446163.46 | 443586.16 | 455964.35 | 458575.32 |
## +---+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 5 | 2021-03-27 | Saturday        | 452004.30           | 446561.59 | 443693.00 | 457478.14 | 460388.41 |
## +---+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 6 | 2021-03-28 | Sunday          | 452958.21           | 446940.72 | 443770.64 | 459013.70 | 462234.65 |
## +---+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 7 | 2021-03-29 | Monday          | 453913.07           | 447301.44 | 443820.01 | 460570.50 | 464113.26 |
## +---+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
paste("Forecasting by using ARIMA Model  ==> ", y_lab , sep=" ")
## [1] "Forecasting by using ARIMA Model  ==>  Cumulative Covid 19 Infection cases in Siberian Federal"
print(ascii(data.frame(FD,forecating_date=forecasting_data_by_name,forecasting_by_auto.arima=tail(forecasting_auto_arima$mean,N_forecasting_days),Lower=tail(forecasting_auto_arima$lower,N_forecasting_days),Upper=tail(forecasting_auto_arima$upper,N_forecasting_days))), type = "rest")
## 
## +---+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## |   | FD         | forecating_date | forecasting_by_auto.arima | Lower.80. | Lower.95. | Upper.80. | Upper.95. |
## +===+============+=================+===========================+===========+===========+===========+===========+
## | 1 | 2021-03-23 | Tuesday         | 447681.43                 | 446768.08 | 446284.57 | 448594.79 | 449078.30 |
## +---+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 2 | 2021-03-24 | Wednesday       | 448529.78                 | 447449.31 | 446877.35 | 449610.26 | 450182.22 |
## +---+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 3 | 2021-03-25 | Thursday        | 449370.82                 | 448109.42 | 447441.67 | 450632.23 | 451299.98 |
## +---+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 4 | 2021-03-26 | Friday          | 450204.72                 | 448748.46 | 447977.57 | 451660.97 | 452431.87 |
## +---+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 5 | 2021-03-27 | Saturday        | 451031.62                 | 449366.51 | 448485.05 | 452696.73 | 453578.19 |
## +---+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 6 | 2021-03-28 | Sunday          | 451851.69                 | 449963.62 | 448964.14 | 453739.75 | 454739.23 |
## +---+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 7 | 2021-03-29 | Monday          | 452665.06                 | 450539.87 | 449414.86 | 454790.26 | 455915.27 |
## +---+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
paste("Forecasting by using NNAR Model  ==> ", y_lab , sep=" ")
## [1] "Forecasting by using NNAR Model  ==>  Cumulative Covid 19 Infection cases in Siberian Federal"
print(ascii(data.frame(FD,forecating_date=forecasting_data_by_name,forecasting_by_NNAR=tail(forecasting_NNAR$mean,N_forecasting_days))), type = "rest")
## 
## +---+------------+-----------------+---------------------+
## |   | FD         | forecating_date | forecasting_by_NNAR |
## +===+============+=================+=====================+
## | 1 | 2021-03-23 | Tuesday         | 445756.45           |
## +---+------------+-----------------+---------------------+
## | 2 | 2021-03-24 | Wednesday       | 446342.16           |
## +---+------------+-----------------+---------------------+
## | 3 | 2021-03-25 | Thursday        | 446911.36           |
## +---+------------+-----------------+---------------------+
## | 4 | 2021-03-26 | Friday          | 447464.36           |
## +---+------------+-----------------+---------------------+
## | 5 | 2021-03-27 | Saturday        | 448001.50           |
## +---+------------+-----------------+---------------------+
## | 6 | 2021-03-28 | Sunday          | 448523.10           |
## +---+------------+-----------------+---------------------+
## | 7 | 2021-03-29 | Monday          | 449029.50           |
## +---+------------+-----------------+---------------------+
result<-c(x1,x2,x3,x4,x5)
table.error<-data.frame(country.name,NNAR.model=MAPE_Mean_All_NNAR, BATS.Model=MAPE_Mean_All.bats_Model,TBATS.Model=MAPE_Mean_All.TBATS_Model,Holt.Model=MAPE_Mean_All.Holt_Model,ARIMA.Model=MAPE_Mean_All.ARIMA_Model,Best.Model=result)
library(ascii)
print(ascii(table(table.error)), type = "rest")
## 
## +---+------------------+------------+------------+-------------+------------+-------------+------------+------+
## |   | country.name     | NNAR.model | BATS.Model | TBATS.Model | Holt.Model | ARIMA.Model | Best.Model | Freq |
## +===+==================+============+============+=============+============+=============+============+======+
## | 1 | Siberian Federal | 0.193      | 0.004      | 0.038       | 0.049      | 0.006       | BATS Model | 1.00 |
## +---+------------------+------------+------------+-------------+------------+-------------+------------+------+
MAPE.Value<-c(MAPE_Mean_All_NNAR,MAPE_Mean_All.bats_Model,MAPE_Mean_All.TBATS_Model,MAPE_Mean_All.Holt_Model,MAPE_Mean_All.ARIMA_Model)
Model<-c("NNAR.model","BATS.Model","TBATS.Model","Holt.Model","ARIMA.Model")
channel_data<-data.frame(Model,MAPE.Value)
# Normally, the entire expression below would be assigned to an object, but we're
# going bare bones here.
ggplot(channel_data, aes(x = Model, y = MAPE.Value)) +
  geom_bar(stat = "identity") +
  geom_text(aes(label = MAPE.Value)) +  # x AND y INHERITED. WE JUST NEED TO SPECIFY "label"
  coord_flip() +
  scale_y_continuous(expand = c(0, 0))

message("System finished Modelling and Forecasting  by using BATS, TBATS, Holt's Linear Trend, and ARIMA Model ==>",y_lab, sep=" ")
## System finished Modelling and Forecasting  by using BATS, TBATS, Holt's Linear Trend, and ARIMA Model ==>Cumulative Covid 19 Infection cases in Siberian Federal
message(" Thank you for using our System For Modelling and Forecasting ==> ",y_lab, sep=" ")
##  Thank you for using our System For Modelling and Forecasting ==> Cumulative Covid 19 Infection cases in Siberian Federal