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

Cumulative vaccination Covid 19 cases In Russian Federation
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("data.xlsx", sheet = "Vaccine") # path of your data ( time series data)
original_data<-Full_original_data$Vaccine
y_lab <- "Cumulative vaccination Covid 19  cases in Russia"   # input name of data
Actual_date_interval <- c("2020/12/16","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 <- "Russia"
# Data Preparation & calculate some of statistics measures
summary(original_data) # Summary your time series
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    3357  854298 2742279 3068390 4797696 8092201
# calculate standard deviation 
data.frame(kurtosis=kurtosis(original_data))   # calculate Cofficient of kurtosis
##   kurtosis
## 1 2.025165
data.frame(skewness=skewness(original_data))  # calculate Cofficient of skewness
##    skewness
## 1 0.4625011
data.frame(Standard.deviation =sd(original_data))
##   Standard.deviation
## 1            2470098
#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 407520513
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 vaccination Covid 19  cases in Russia"
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 vaccination Covid 19  cases in Russia"
paste(MAPE_Mean_All,"%")
## [1] "4.747 % MAPE  8 days Cumulative vaccination Covid 19  cases in Russia %"
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 vaccination Covid 19  cases in Russia"
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                  | 7344639.00  | 7270453.46       | 1.01 %          |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 2 | 2021-03-16 | Tuesday                 | 7499489.00  | 7323941.07       | 2.341 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 3 | 2021-03-17 | Wednesday               | 7609472.00  | 7359186.35       | 3.289 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 4 | 2021-03-18 | Thursday                | 7700258.00  | 7381895.12       | 4.134 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 5 | 2021-03-19 | Friday                  | 7805641.00  | 7396311.22       | 5.244 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 6 | 2021-03-20 | Saturday                | 7900944.00  | 7405375.85       | 6.272 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 7 | 2021-03-21 | Sunday                  | 7995756.00  | 7411041.06       | 7.313 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 8 | 2021-03-22 | Monday                  | 8092201.00  | 7414568.21       | 8.374 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
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         | 7416758.97          |
## +---+------------+-----------------+---------------------+
## | 2 | 2021-03-24 | Wednesday       | 7418117.66          |
## +---+------------+-----------------+---------------------+
## | 3 | 2021-03-25 | Thursday        | 7418959.53          |
## +---+------------+-----------------+---------------------+
## | 4 | 2021-03-26 | Friday          | 7419480.87          |
## +---+------------+-----------------+---------------------+
## | 5 | 2021-03-27 | Saturday        | 7419803.61          |
## +---+------------+-----------------+---------------------+
## | 6 | 2021-03-28 | Sunday          | 7420003.35          |
## +---+------------+-----------------+---------------------+
## | 7 | 2021-03-29 | Monday          | 7420126.96          |
## +---+------------+-----------------+---------------------+
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 1548.089 17746.45 7972.698 16.3128 49.46082 0.09759931 0.005464383
# Print Model Parameters
model_bats
## BATS(1, {0,0}, 1, -)
## 
## Call: bats(y = data_series)
## 
## Parameters
##   Lambda: 1
##   Alpha: 0.8716822
##   Beta: 1.044995
##   Damping Parameter: 1
## 
## Seed States:
##            [,1]
## [1,] -93019.016
## [2,]   3860.791
## attr(,"lambda")
## [1] 1
## 
## Sigma: 17746.44
## AIC: 2151.03
#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 vaccination Covid 19  cases in Russia"
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 vaccination Covid 19  cases in Russia"
paste(MAPE_Mean_All.bats,"%")
## [1] "1.485 % MAPE  8 days Cumulative vaccination Covid 19  cases in Russia %"
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 vaccination Covid 19  cases in Russia"
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                  | 7344639.00  | 7339702.11       | 0.067 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 2 | 2021-03-16 | Tuesday                 | 7499489.00  | 7487542.31       | 0.159 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 3 | 2021-03-17 | Wednesday               | 7609472.00  | 7635382.51       | 0.341 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 4 | 2021-03-18 | Thursday                | 7700258.00  | 7783222.72       | 1.077 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 5 | 2021-03-19 | Friday                  | 7805641.00  | 7931062.92       | 1.607 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 6 | 2021-03-20 | Saturday                | 7900944.00  | 8078903.12       | 2.252 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 7 | 2021-03-21 | Sunday                  | 7995756.00  | 8226743.32       | 2.889 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 8 | 2021-03-22 | Monday                  | 8092201.00  | 8374583.52       | 3.49 %          |
## +---+------------+-------------------------+-------------+------------------+-----------------+
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         | 8522423.73          | 8131512.25 | 7924576.21 | 8131512.25 | 7924576.21 |
## +---+------------+-----------------+---------------------+------------+------------+------------+------------+
## | 2 | 2021-03-24 | Wednesday       | 8670263.93          | 8214810.81 | 7973708.50 | 8214810.81 | 7973708.50 |
## +---+------------+-----------------+---------------------+------------+------------+------------+------------+
## | 3 | 2021-03-25 | Thursday        | 8818104.13          | 8294904.80 | 8017939.81 | 8294904.80 | 8017939.81 |
## +---+------------+-----------------+---------------------+------------+------------+------------+------------+
## | 4 | 2021-03-26 | Friday          | 8965944.33          | 8371939.75 | 8057492.72 | 8371939.75 | 8057492.72 |
## +---+------------+-----------------+---------------------+------------+------------+------------+------------+
## | 5 | 2021-03-27 | Saturday        | 9113784.53          | 8446042.86 | 8092561.77 | 8446042.86 | 8092561.77 |
## +---+------------+-----------------+---------------------+------------+------------+------------+------------+
## | 6 | 2021-03-28 | Sunday          | 9261624.74          | 8517326.60 | 8123318.97 | 8517326.60 | 8123318.97 |
## +---+------------+-----------------+---------------------+------------+------------+------------+------------+
## | 7 | 2021-03-29 | Monday          | 9409464.94          | 8585891.37 | 8149917.85 | 8585891.37 | 8149917.85 |
## +---+------------+-----------------+---------------------+------------+------------+------------+------------+
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
## Training set 3386.069 36603.97 18568.06 50.61896 120.5399 0.2273044
##                      ACF1
## Training set -0.002627535
# Print Model Parameters
model_TBATS
## TBATS(1, {0,0}, 1, {<6,2>})
## 
## Call: NULL
## 
## Parameters
##   Alpha: 0.9842787
##   Beta: 0.5003967
##   Damping Parameter: 1
##   Gamma-1 Values: -0.00336074
##   Gamma-2 Values: 0.0005554326
## 
## Seed States:
##              [,1]
## [1,] -263682.8534
## [2,]   11966.3201
## [3,]     454.0979
## [4,]    3937.1547
## [5,]    9238.1052
## [6,]    2660.4872
## 
## Sigma: 36603.97
## AIC: 2289.897
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 vaccination Covid 19  cases in Russia"
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 vaccination Covid 19  cases in Russia"
paste(MAPE_Mean_All.TBATS,"%")
## [1] "2.375 % MAPE  8 days Cumulative vaccination Covid 19  cases in Russia %"
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 vaccination Covid 19  cases in Russia"
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                  | 7344639.00  | 7350189.59        | 0.076 %          |
## +---+------------+-------------------------+-------------+-------------------+------------------+
## | 2 | 2021-03-16 | Tuesday                 | 7499489.00  | 7529721.12        | 0.403 %          |
## +---+------------+-------------------------+-------------+-------------------+------------------+
## | 3 | 2021-03-17 | Wednesday               | 7609472.00  | 7697515.57        | 1.157 %          |
## +---+------------+-------------------------+-------------+-------------------+------------------+
## | 4 | 2021-03-18 | Thursday                | 7700258.00  | 7854985.79        | 2.009 %          |
## +---+------------+-------------------------+-------------+-------------------+------------------+
## | 5 | 2021-03-19 | Friday                  | 7805641.00  | 8017295.73        | 2.712 %          |
## +---+------------+-------------------------+-------------+-------------------+------------------+
## | 6 | 2021-03-20 | Saturday                | 7900944.00  | 8168591.80        | 3.388 %          |
## +---+------------+-------------------------+-------------+-------------------+------------------+
## | 7 | 2021-03-21 | Sunday                  | 7995756.00  | 8326785.22        | 4.14 %           |
## +---+------------+-------------------------+-------------+-------------------+------------------+
## | 8 | 2021-03-22 | Monday                  | 8092201.00  | 8506316.75        | 5.117 %          |
## +---+------------+-------------------------+-------------+-------------------+------------------+
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         | 8674111.21           | 8537295.42 | 8464869.52 | 8810926.99 | 8883352.90 |
## +---+------------+-----------------+----------------------+------------+------------+------------+------------+
## | 2 | 2021-03-24 | Wednesday       | 8831581.43           | 8687604.50 | 8611387.72 | 8975558.35 | 9051775.13 |
## +---+------------+-----------------+----------------------+------------+------------+------------+------------+
## | 3 | 2021-03-25 | Thursday        | 8993891.36           | 8843111.51 | 8763293.47 | 9144671.22 | 9224489.26 |
## +---+------------+-----------------+----------------------+------------+------------+------------+------------+
## | 4 | 2021-03-26 | Friday          | 9145187.43           | 8987983.22 | 8904764.33 | 9302391.65 | 9385610.53 |
## +---+------------+-----------------+----------------------+------------+------------+------------+------------+
## | 5 | 2021-03-27 | Saturday        | 9303380.86           | 9140090.28 | 9053649.46 | 9466671.44 | 9553112.25 |
## +---+------------+-----------------+----------------------+------------+------------+------------+------------+
## | 6 | 2021-03-28 | Sunday          | 9482912.39           | 9313799.85 | 9224277.08 | 9652024.93 | 9741547.70 |
## +---+------------+-----------------+----------------------+------------+------------+------------+------------+
## | 7 | 2021-03-29 | Monday          | 9650706.84           | 9475990.19 | 9383500.79 | 9825423.49 | 9917912.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 -26.93708 11291.17 6149.15 -0.4359066 1.814672 0.075276 0.5493567
# 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.6147 
## 
##   Smoothing parameters:
##     alpha = 0.9984 
##     beta  = 0.9984 
## 
##   Initial states:
##     l = -17.3592 
##     b = 260.4519 
## 
##   sigma:  44.2008
## 
##      AIC     AICc      BIC 
## 1079.792 1080.515 1092.235 
## 
## Training set error measures:
##                     ME     RMSE     MAE        MPE     MAPE     MASE      ACF1
## Training set -26.93708 11291.17 6149.15 -0.4359066 1.814672 0.075276 0.5493567
# 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 vaccination Covid 19  cases in Russia"
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 vaccination Covid 19  cases in Russia"
paste(MAPE_Mean_All.Holt,"%")
## [1] "1.752 % MAPE  8 days Cumulative vaccination Covid 19  cases in Russia %"
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 vaccination Covid 19  cases in Russia"
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                  | 7344639.00  | 7342149.28       | 0.034 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 2 | 2021-03-16 | Tuesday                 | 7499489.00  | 7493581.89       | 0.079 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 3 | 2021-03-17 | Wednesday               | 7609472.00  | 7646203.02       | 0.483 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 4 | 2021-03-18 | Thursday                | 7700258.00  | 7800007.19       | 1.295 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 5 | 2021-03-19 | Friday                  | 7805641.00  | 7954989.02       | 1.913 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 6 | 2021-03-20 | Saturday                | 7900944.00  | 8111143.19       | 2.66 %          |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 7 | 2021-03-21 | Sunday                  | 7995756.00  | 8268464.50       | 3.411 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 8 | 2021-03-22 | Monday                  | 8092201.00  | 8426947.82       | 4.137 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
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         | 8586588.09          | 8142575.66 | 7911222.34 | 9039629.93  | 9283074.59  |
## +---+------------+-----------------+---------------------+------------+------------+-------------+-------------+
## | 2 | 2021-03-24 | Wednesday       | 8747380.35          | 8228409.84 | 7958658.33 | 9278496.85  | 9564513.64  |
## +---+------------+-----------------+---------------------+------------+------------+-------------+-------------+
## | 3 | 2021-03-25 | Thursday        | 8909319.71          | 8311085.97 | 8000921.64 | 9523449.97  | 9854902.93  |
## +---+------------+-----------------+---------------------+------------+------------+-------------+-------------+
## | 4 | 2021-03-26 | Friday          | 9072401.35          | 8390718.98 | 8038214.64 | 9774419.75  | 10154159.33 |
## +---+------------+-----------------+---------------------+------------+------------+-------------+-------------+
## | 5 | 2021-03-27 | Saturday        | 9236620.54          | 8467408.74 | 8070716.53 | 10031351.07 | 10462221.27 |
## +---+------------+-----------------+---------------------+------------+------------+-------------+-------------+
## | 6 | 2021-03-28 | Sunday          | 9401972.59          | 8541243.08 | 8098587.96 | 10294200.22 | 10779044.12 |
## +---+------------+-----------------+---------------------+------------+------------+-------------+-------------+
## | 7 | 2021-03-29 | Monday          | 9568452.92          | 8612300.06 | 8121974.51 | 10562932.69 | 11104596.88 |
## +---+------------+-----------------+---------------------+------------+------------+-------------+-------------+
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 vaccination Covid 19  cases in Russia"
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 = 2.2785, Truncation lag parameter = 3, p-value = 0.01
pp.test(data_series)   # applay pp test
## Warning in pp.test(data_series): p-value greater than printed p-value
## 
##  Phillips-Perron Unit Root Test
## 
## data:  data_series
## Dickey-Fuller Z(alpha) = -0.65233, Truncation lag parameter = 3,
## p-value = 0.99
## alternative hypothesis: stationary
adf.test(data_series)  # applay adf test
## 
##  Augmented Dickey-Fuller Test
## 
## data:  data_series
## Dickey-Fuller = -0.58768, Lag order = 4, p-value = 0.976
## 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 vaccination Covid 19  cases in Russia"
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 = 1.2711, Truncation lag parameter = 3, p-value = 0.01
pp.test(diff1_x1)     # applay pp test after taking first differences
## 
##  Phillips-Perron Unit Root Test
## 
## data:  diff1_x1
## Dickey-Fuller Z(alpha) = -10.76, Truncation lag parameter = 3, p-value
## = 0.4862
## alternative hypothesis: stationary
adf.test(diff1_x1)    # applay adf test after taking first differences
## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff1_x1
## Dickey-Fuller = -3.5844, Lag order = 4, p-value = 0.03935
## 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 vaccination Covid 19  cases in Russia"
kpss.test(diff2_x1)   # applay kpss test after taking Second differences
## Warning in kpss.test(diff2_x1): p-value greater than printed p-value
## 
##  KPSS Test for Level Stationarity
## 
## data:  diff2_x1
## KPSS Level = 0.049076, Truncation lag parameter = 3, p-value = 0.1
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) = -38.546, Truncation lag parameter = 3, 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.5403, Lag order = 4, p-value = 0.04332
## 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)                    : 1875.29
##  ARIMA(0,2,1)                    : 1852.518
##  ARIMA(0,2,2)                    : 1848.255
##  ARIMA(0,2,3)                    : 1848.912
##  ARIMA(0,2,4)                    : 1844.181
##  ARIMA(0,2,5)                    : Inf
##  ARIMA(1,2,0)                    : 1844.406
##  ARIMA(1,2,1)                    : 1846.544
##  ARIMA(1,2,2)                    : 1848.66
##  ARIMA(1,2,3)                    : 1850.911
##  ARIMA(1,2,4)                    : 1838.57
##  ARIMA(2,2,0)                    : 1846.543
##  ARIMA(2,2,1)                    : Inf
##  ARIMA(2,2,2)                    : Inf
##  ARIMA(2,2,3)                    : 1842.612
##  ARIMA(3,2,0)                    : 1848.644
##  ARIMA(3,2,1)                    : Inf
##  ARIMA(3,2,2)                    : Inf
##  ARIMA(4,2,0)                    : 1850.896
##  ARIMA(4,2,1)                    : Inf
##  ARIMA(5,2,0)                    : 1852.281
## 
## 
## 
##  Best model: ARIMA(1,2,4)
model1 # show the result of autoarima 
## Series: data_series 
## ARIMA(1,2,4) 
## 
## Coefficients:
##           ar1     ma1     ma2     ma3     ma4
##       -0.7740  1.5186  0.8987  0.3513  0.2825
## s.e.   0.1193  0.1511  0.2162  0.1931  0.1060
## 
## sigma^2 estimated as 78032351:  log likelihood=-912.76
## AIC=1837.52   AICc=1838.57   BIC=1852.32
#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 4
strtoi(bestmodel[3])
## [1] 4
#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     ma3     ma4
##       -0.7740  1.5186  0.8987  0.3513  0.2825
## s.e.   0.1193  0.1511  0.2162  0.1931  0.1060
## 
## sigma^2 estimated as 73547733:  log likelihood = -912.76,  aic = 1837.52
paste ("accuracy of autoarima Model For  ==> ",y_lab, sep=" ")
## [1] "accuracy of autoarima Model For  ==>  Cumulative vaccination Covid 19  cases in Russia"
accuracy(x1_model1)  # aacuracy of best model from auto arima
##                    ME    RMSE     MAE       MPE      MAPE      MASE
## Training set 667.1011 8479.09 3850.22 0.4710244 0.6254332 0.0471332
##                      ACF1
## Training set 0.0003793958
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,4)
## Q* = 13.064, df = 5, p-value = 0.02279
## 
## Model df: 5.   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 vaccination Covid 19  cases in Russia"
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 = 31.095, df = 20, p-value = 0.05395
library(tseries)
jarque.bera.test(x1_model1$residuals)  # Do test jarque.bera.test 
## 
##  Jarque Bera Test
## 
## data:  x1_model1$residuals
## X-squared = 588.61, 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 vaccination Covid 19  cases in Russia"
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 vaccination Covid 19  cases in Russia"
paste(MAPE_Mean_All.ARIMA,"%")
## [1] "1.081 % MAPE  8 days Cumulative vaccination Covid 19  cases in Russia %"
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 vaccination Covid 19  cases in Russia"
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                  | 7344639.00  | 7328249.75             | 0.223 %               |
## +---+-----------------+-------------------------+-------------+------------------------+-----------------------+
## | 2 | 2021-03-16      | Tuesday                 | 7499489.00  | 7471398.11             | 0.375 %               |
## +---+-----------------+-------------------------+-------------+------------------------+-----------------------+
## | 3 | 2021-03-17      | Wednesday               | 7609472.00  | 7608716.90             | 0.01 %                |
## +---+-----------------+-------------------------+-------------+------------------------+-----------------------+
## | 4 | 2021-03-18      | Thursday                | 7700258.00  | 7749252.94             | 0.636 %               |
## +---+-----------------+-------------------------+-------------+------------------------+-----------------------+
## | 5 | 2021-03-19      | Friday                  | 7805641.00  | 7887298.91             | 1.046 %               |
## +---+-----------------+-------------------------+-------------+------------------------+-----------------------+
## | 6 | 2021-03-20      | Saturday                | 7900944.00  | 8027272.13             | 1.599 %               |
## +---+-----------------+-------------------------+-------------+------------------------+-----------------------+
## | 7 | 2021-03-21      | Sunday                  | 7995756.00  | 8165753.70             | 2.126 %               |
## +---+-----------------+-------------------------+-------------+------------------------+-----------------------+
## | 8 | 2021-03-22      | Monday                  | 8092201.00  | 8305389.77             | 2.634 %               |
## +---+-----------------+-------------------------+-------------+------------------------+-----------------------+
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         | 8444132.29                | 8080171.51 | 7887502.30 | 8808093.08  | 9000762.28  |
## +---+------------+-----------------+---------------------------+------------+------------+-------------+-------------+
## | 2 | 2021-03-24 | Wednesday       | 8583566.40                | 8154224.73 | 7926944.98 | 9012908.06  | 9240187.81  |
## +---+------------+-----------------+---------------------------+------------+------------+-------------+-------------+
## | 3 | 2021-03-25 | Thursday        | 8722465.23                | 8224098.64 | 7960279.31 | 9220831.82  | 9484651.16  |
## +---+------------+-----------------+---------------------------+------------+------------+-------------+-------------+
## | 4 | 2021-03-26 | Friday          | 8861778.35                | 8291124.16 | 7989038.07 | 9432432.55  | 9734518.63  |
## +---+------------+-----------------+---------------------------+------------+------------+-------------+-------------+
## | 5 | 2021-03-27 | Saturday        | 9000770.83                | 8354563.82 | 8012482.49 | 9646977.84  | 9989059.16  |
## +---+------------+-----------------+---------------------------+------------+------------+-------------+-------------+
## | 6 | 2021-03-28 | Sunday          | 9140011.47                | 8415230.74 | 8031555.00 | 9864792.21  | 10248467.95 |
## +---+------------+-----------------+---------------------------+------------+------------+-------------+-------------+
## | 7 | 2021-03-29 | Monday          | 9279060.04                | 8472716.93 | 8045864.68 | 10085403.15 | 10512255.41 |
## +---+------------+-----------------+---------------------------+------------+------------+-------------+-------------+
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] "1.081 % MAPE  8 days Cumulative vaccination Covid 19  cases in Russia"
# 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 vaccination Covid 19  cases in Russia"
best_recommended_model
## [1] 1.081
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 vaccination Covid 19  cases in Russia"
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         | 8522423.73          | 8131512.25 | 7924576.21 | 8131512.25 | 7924576.21 |
## +---+------------+-----------------+---------------------+------------+------------+------------+------------+
## | 2 | 2021-03-24 | Wednesday       | 8670263.93          | 8214810.81 | 7973708.50 | 8214810.81 | 7973708.50 |
## +---+------------+-----------------+---------------------+------------+------------+------------+------------+
## | 3 | 2021-03-25 | Thursday        | 8818104.13          | 8294904.80 | 8017939.81 | 8294904.80 | 8017939.81 |
## +---+------------+-----------------+---------------------+------------+------------+------------+------------+
## | 4 | 2021-03-26 | Friday          | 8965944.33          | 8371939.75 | 8057492.72 | 8371939.75 | 8057492.72 |
## +---+------------+-----------------+---------------------+------------+------------+------------+------------+
## | 5 | 2021-03-27 | Saturday        | 9113784.53          | 8446042.86 | 8092561.77 | 8446042.86 | 8092561.77 |
## +---+------------+-----------------+---------------------+------------+------------+------------+------------+
## | 6 | 2021-03-28 | Sunday          | 9261624.74          | 8517326.60 | 8123318.97 | 8517326.60 | 8123318.97 |
## +---+------------+-----------------+---------------------+------------+------------+------------+------------+
## | 7 | 2021-03-29 | Monday          | 9409464.94          | 8585891.37 | 8149917.85 | 8585891.37 | 8149917.85 |
## +---+------------+-----------------+---------------------+------------+------------+------------+------------+
paste("Forecasting by using TBATS Model  ==> ", y_lab , sep=" ")
## [1] "Forecasting by using TBATS Model  ==>  Cumulative vaccination Covid 19  cases in Russia"
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         | 8674111.21           | 8537295.42 | 8464869.52 | 8810926.99 | 8883352.90 |
## +---+------------+-----------------+----------------------+------------+------------+------------+------------+
## | 2 | 2021-03-24 | Wednesday       | 8831581.43           | 8687604.50 | 8611387.72 | 8975558.35 | 9051775.13 |
## +---+------------+-----------------+----------------------+------------+------------+------------+------------+
## | 3 | 2021-03-25 | Thursday        | 8993891.36           | 8843111.51 | 8763293.47 | 9144671.22 | 9224489.26 |
## +---+------------+-----------------+----------------------+------------+------------+------------+------------+
## | 4 | 2021-03-26 | Friday          | 9145187.43           | 8987983.22 | 8904764.33 | 9302391.65 | 9385610.53 |
## +---+------------+-----------------+----------------------+------------+------------+------------+------------+
## | 5 | 2021-03-27 | Saturday        | 9303380.86           | 9140090.28 | 9053649.46 | 9466671.44 | 9553112.25 |
## +---+------------+-----------------+----------------------+------------+------------+------------+------------+
## | 6 | 2021-03-28 | Sunday          | 9482912.39           | 9313799.85 | 9224277.08 | 9652024.93 | 9741547.70 |
## +---+------------+-----------------+----------------------+------------+------------+------------+------------+
## | 7 | 2021-03-29 | Monday          | 9650706.84           | 9475990.19 | 9383500.79 | 9825423.49 | 9917912.90 |
## +---+------------+-----------------+----------------------+------------+------------+------------+------------+
paste("Forecasting by using Holt's Linear Trend Model  ==> ", y_lab , sep=" ")
## [1] "Forecasting by using Holt's Linear Trend Model  ==>  Cumulative vaccination Covid 19  cases in Russia"
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         | 8586588.09          | 8142575.66 | 7911222.34 | 9039629.93  | 9283074.59  |
## +---+------------+-----------------+---------------------+------------+------------+-------------+-------------+
## | 2 | 2021-03-24 | Wednesday       | 8747380.35          | 8228409.84 | 7958658.33 | 9278496.85  | 9564513.64  |
## +---+------------+-----------------+---------------------+------------+------------+-------------+-------------+
## | 3 | 2021-03-25 | Thursday        | 8909319.71          | 8311085.97 | 8000921.64 | 9523449.97  | 9854902.93  |
## +---+------------+-----------------+---------------------+------------+------------+-------------+-------------+
## | 4 | 2021-03-26 | Friday          | 9072401.35          | 8390718.98 | 8038214.64 | 9774419.75  | 10154159.33 |
## +---+------------+-----------------+---------------------+------------+------------+-------------+-------------+
## | 5 | 2021-03-27 | Saturday        | 9236620.54          | 8467408.74 | 8070716.53 | 10031351.07 | 10462221.27 |
## +---+------------+-----------------+---------------------+------------+------------+-------------+-------------+
## | 6 | 2021-03-28 | Sunday          | 9401972.59          | 8541243.08 | 8098587.96 | 10294200.22 | 10779044.12 |
## +---+------------+-----------------+---------------------+------------+------------+-------------+-------------+
## | 7 | 2021-03-29 | Monday          | 9568452.92          | 8612300.06 | 8121974.51 | 10562932.69 | 11104596.88 |
## +---+------------+-----------------+---------------------+------------+------------+-------------+-------------+
paste("Forecasting by using ARIMA Model  ==> ", y_lab , sep=" ")
## [1] "Forecasting by using ARIMA Model  ==>  Cumulative vaccination Covid 19  cases in Russia"
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         | 8444132.29                | 8080171.51 | 7887502.30 | 8808093.08  | 9000762.28  |
## +---+------------+-----------------+---------------------------+------------+------------+-------------+-------------+
## | 2 | 2021-03-24 | Wednesday       | 8583566.40                | 8154224.73 | 7926944.98 | 9012908.06  | 9240187.81  |
## +---+------------+-----------------+---------------------------+------------+------------+-------------+-------------+
## | 3 | 2021-03-25 | Thursday        | 8722465.23                | 8224098.64 | 7960279.31 | 9220831.82  | 9484651.16  |
## +---+------------+-----------------+---------------------------+------------+------------+-------------+-------------+
## | 4 | 2021-03-26 | Friday          | 8861778.35                | 8291124.16 | 7989038.07 | 9432432.55  | 9734518.63  |
## +---+------------+-----------------+---------------------------+------------+------------+-------------+-------------+
## | 5 | 2021-03-27 | Saturday        | 9000770.83                | 8354563.82 | 8012482.49 | 9646977.84  | 9989059.16  |
## +---+------------+-----------------+---------------------------+------------+------------+-------------+-------------+
## | 6 | 2021-03-28 | Sunday          | 9140011.47                | 8415230.74 | 8031555.00 | 9864792.21  | 10248467.95 |
## +---+------------+-----------------+---------------------------+------------+------------+-------------+-------------+
## | 7 | 2021-03-29 | Monday          | 9279060.04                | 8472716.93 | 8045864.68 | 10085403.15 | 10512255.41 |
## +---+------------+-----------------+---------------------------+------------+------------+-------------+-------------+
paste("Forecasting by using NNAR Model  ==> ", y_lab , sep=" ")
## [1] "Forecasting by using NNAR Model  ==>  Cumulative vaccination Covid 19  cases in Russia"
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         | 7416758.97          |
## +---+------------+-----------------+---------------------+
## | 2 | 2021-03-24 | Wednesday       | 7418117.66          |
## +---+------------+-----------------+---------------------+
## | 3 | 2021-03-25 | Thursday        | 7418959.53          |
## +---+------------+-----------------+---------------------+
## | 4 | 2021-03-26 | Friday          | 7419480.87          |
## +---+------------+-----------------+---------------------+
## | 5 | 2021-03-27 | Saturday        | 7419803.61          |
## +---+------------+-----------------+---------------------+
## | 6 | 2021-03-28 | Sunday          | 7420003.35          |
## +---+------------+-----------------+---------------------+
## | 7 | 2021-03-29 | Monday          | 7420126.96          |
## +---+------------+-----------------+---------------------+
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 | Russia       | 4.747      | 1.485      | 2.375       | 1.752      | 1.081       | ARIMA 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 vaccination Covid 19  cases in Russia
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 vaccination Covid 19  cases in Russia