Разработка алгоритмов выбора лучших моделей временных рядов и нейронных сетей для прогнозирования случаев COVID-19

Cumulative Covid 19 Infections cases In Chelyabinsk

Makarovskikh Tatyana Anatolyevna “Макаровских Татьяна Анатольевна”

Abotaleb mostafa“Аботалеб Мостафа”

Faculty of Electrical Engineering and Computer Science

Department of system programming

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.csv("data.csv") # path of your data ( time series data)
original_data<-Full_original_data$Infections
y_lab <- "Forecasting Cumulative Covid 19 infection cases in Chelyabinsk"   # input name of data
Actual_date_interval <- c("2020/03/12","2021/04/09")
Forecast_date_interval <- c("2021/04/10","2021/04/30")
validation_data_days <-7
frequency<-"day"
Number_Neural<-5   # Number of Neural For model NNAR Model
NNAR_Model<- FALSE     #create new model (TRUE/FALSE)
frequency<-"days"
Population <-1130319 # population in England for SIR Model
country.name <- "Chelyabinsk"

# Data Preparation & calculate some of statistics measures
summary(original_data) # Summary your time series

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       0    5166   15331   20702   34932   56513

# calculate standard deviation 
data.frame(kurtosis=kurtosis(original_data))   # calculate Cofficient of kurtosis

##   kurtosis
## 1 2.090858

data.frame(skewness=skewness(original_data))  # calculate Cofficient of skewness

##    skewness
## 1 0.6291764

data.frame(Standard.deviation =sd(original_data))

##   Standard.deviation
## 1           17959.76

#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
}
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
}

## 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 916.1

# 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  7 days by using NNAR Model for  ==>  Forecasting Cumulative Covid 19 infection cases in Chelyabinsk"

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  7  days in NNAR Model for  ==>  Forecasting Cumulative Covid 19 infection cases in Chelyabinsk"

paste(MAPE_Mean_All,"%")

## [1] "0.458 % MAPE  7 days Forecasting Cumulative Covid 19 infection cases in Chelyabinsk %"

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  7  days in NNAR Model for  ==>  Forecasting Cumulative Covid 19 infection cases in Chelyabinsk"

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-04-03 | Saturday                | 55783.00    | 55719.97         | 0.113 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 2 | 2021-04-04 | Sunday                  | 55908.00    | 55782.10         | 0.225 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 3 | 2021-04-05 | Monday                  | 56031.00    | 55841.47         | 0.338 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 4 | 2021-04-06 | Tuesday                 | 56153.00    | 55898.19         | 0.454 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 5 | 2021-04-07 | Wednesday               | 56274.00    | 55952.34         | 0.572 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 6 | 2021-04-08 | Thursday                | 56394.00    | 56004.02         | 0.692 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 7 | 2021-04-09 | Friday                  | 56513.00    | 56053.32         | 0.813 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+

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-04-10 | Saturday        | 56100.34            |
## +----+------------+-----------------+---------------------+
## | 2  | 2021-04-11 | Sunday          | 56145.15            |
## +----+------------+-----------------+---------------------+
## | 3  | 2021-04-12 | Monday          | 56187.85            |
## +----+------------+-----------------+---------------------+
## | 4  | 2021-04-13 | Tuesday         | 56228.52            |
## +----+------------+-----------------+---------------------+
## | 5  | 2021-04-14 | Wednesday       | 56267.25            |
## +----+------------+-----------------+---------------------+
## | 6  | 2021-04-15 | Thursday        | 56304.11            |
## +----+------------+-----------------+---------------------+
## | 7  | 2021-04-16 | Friday          | 56339.19            |
## +----+------------+-----------------+---------------------+
## | 8  | 2021-04-17 | Saturday        | 56372.56            |
## +----+------------+-----------------+---------------------+
## | 9  | 2021-04-18 | Sunday          | 56404.30            |
## +----+------------+-----------------+---------------------+
## | 10 | 2021-04-19 | Monday          | 56434.47            |
## +----+------------+-----------------+---------------------+
## | 11 | 2021-04-20 | Tuesday         | 56463.14            |
## +----+------------+-----------------+---------------------+
## | 12 | 2021-04-21 | Wednesday       | 56490.39            |
## +----+------------+-----------------+---------------------+
## | 13 | 2021-04-22 | Thursday        | 56516.28            |
## +----+------------+-----------------+---------------------+
## | 14 | 2021-04-23 | Friday          | 56540.87            |
## +----+------------+-----------------+---------------------+
## | 15 | 2021-04-24 | Saturday        | 56564.22            |
## +----+------------+-----------------+---------------------+
## | 16 | 2021-04-25 | Sunday          | 56586.39            |
## +----+------------+-----------------+---------------------+
## | 17 | 2021-04-26 | Monday          | 56607.44            |
## +----+------------+-----------------+---------------------+
## | 18 | 2021-04-27 | Tuesday         | 56627.41            |
## +----+------------+-----------------+---------------------+
## | 19 | 2021-04-28 | Wednesday       | 56646.36            |
## +----+------------+-----------------+---------------------+
## | 20 | 2021-04-29 | Thursday        | 56664.34            |
## +----+------------+-----------------+---------------------+
## | 21 | 2021-04-30 | Friday          | 56681.39            |
## +----+------------+-----------------+---------------------+

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 0.4316178 11.44078 6.540191 NaN  Inf 0.04536005 -0.002938951

# Print Model Parameters
model_bats

## BATS(1, {0,0}, 1, -)
## 
## Call: bats(y = data_series)
## 
## Parameters
##   Alpha: 0.985641
##   Beta: 0.8084883
##   Damping Parameter: 1
## 
## Seed States:
##           [,1]
## [1,] -6.177343
## [2,] -5.825196
## 
## Sigma: 11.44078
## AIC: 4200.291

#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  7 days by using bats Model for  ==>  Forecasting Cumulative Covid 19 infection cases in Chelyabinsk"

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  7  days in bats Model for  ==>  Forecasting Cumulative Covid 19 infection cases in Chelyabinsk"

paste(MAPE_Mean_All.bats,"%")

## [1] "0.037 % MAPE  7 days Forecasting Cumulative Covid 19 infection cases in Chelyabinsk %"

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  7  days in bats Model for  ==>  Forecasting Cumulative Covid 19 infection cases in Chelyabinsk"

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-04-03 | Saturday                | 55783.00    | 55784.24         | 0.002 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 2 | 2021-04-04 | Sunday                  | 55908.00    | 55913.46         | 0.01 %          |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 3 | 2021-04-05 | Monday                  | 56031.00    | 56042.68         | 0.021 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 4 | 2021-04-06 | Tuesday                 | 56153.00    | 56171.90         | 0.034 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 5 | 2021-04-07 | Wednesday               | 56274.00    | 56301.13         | 0.048 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 6 | 2021-04-08 | Thursday                | 56394.00    | 56430.35         | 0.064 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 7 | 2021-04-09 | Friday                  | 56513.00    | 56559.57         | 0.082 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+

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-04-10 | Saturday        | 56688.79            | 56512.89  | 56419.77  | 56512.89  | 56419.77  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 2  | 2021-04-11 | Sunday          | 56818.01            | 56610.92  | 56501.30  | 56610.92  | 56501.30  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 3  | 2021-04-12 | Monday          | 56947.23            | 56707.32  | 56580.31  | 56707.32  | 56580.31  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 4  | 2021-04-13 | Tuesday         | 57076.45            | 56802.14  | 56656.93  | 56802.14  | 56656.93  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 5  | 2021-04-14 | Wednesday       | 57205.68            | 56895.47  | 56731.26  | 56895.47  | 56731.26  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 6  | 2021-04-15 | Thursday        | 57334.90            | 56987.36  | 56803.39  | 56987.36  | 56803.39  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 7  | 2021-04-16 | Friday          | 57464.12            | 57077.87  | 56873.40  | 57077.87  | 56873.40  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 8  | 2021-04-17 | Saturday        | 57593.34            | 57167.03  | 56941.36  | 57167.03  | 56941.36  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 9  | 2021-04-18 | Sunday          | 57722.56            | 57254.91  | 57007.34  | 57254.91  | 57007.34  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 10 | 2021-04-19 | Monday          | 57851.78            | 57341.52  | 57071.41  | 57341.52  | 57071.41  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 11 | 2021-04-20 | Tuesday         | 57981.01            | 57426.92  | 57133.61  | 57426.92  | 57133.61  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 12 | 2021-04-21 | Wednesday       | 58110.23            | 57511.13  | 57193.99  | 57511.13  | 57193.99  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 13 | 2021-04-22 | Thursday        | 58239.45            | 57594.19  | 57252.61  | 57594.19  | 57252.61  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 14 | 2021-04-23 | Friday          | 58368.67            | 57676.12  | 57309.51  | 57676.12  | 57309.51  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 15 | 2021-04-24 | Saturday        | 58497.89            | 57756.95  | 57364.72  | 57756.95  | 57364.72  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 16 | 2021-04-25 | Sunday          | 58627.11            | 57836.70  | 57418.28  | 57836.70  | 57418.28  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 17 | 2021-04-26 | Monday          | 58756.33            | 57915.40  | 57470.23  | 57915.40  | 57470.23  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 18 | 2021-04-27 | Tuesday         | 58885.56            | 57993.06  | 57520.60  | 57993.06  | 57520.60  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 19 | 2021-04-28 | Wednesday       | 59014.78            | 58069.71  | 57569.43  | 58069.71  | 57569.43  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 20 | 2021-04-29 | Thursday        | 59144.00            | 58145.37  | 57616.73  | 58145.37  | 57616.73  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 21 | 2021-04-30 | Friday          | 59273.22            | 58220.06  | 57662.55  | 58220.06  | 57662.55  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+

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 0.4294133 11.34378 6.647381 NaN  Inf 0.04610347 -0.002703582

# Print Model Parameters
model_TBATS

## TBATS(1, {0,0}, 1, {<6,2>})
## 
## Call: NULL
## 
## Parameters
##   Alpha: 0.9794355
##   Beta: 0.8160431
##   Damping Parameter: 1
##   Gamma-1 Values: -0.001024159
##   Gamma-2 Values: -0.0002236461
## 
## Seed States:
##            [,1]
## [1,] -6.3221139
## [2,] -5.7677411
## [3,]  0.5029448
## [4,] -0.3633024
## [5,]  1.0302208
## [6,]  0.3373406
## 
## Sigma: 11.34378
## AIC: 4205.7

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  7 days by using TBATS Model for  ==>  Forecasting Cumulative Covid 19 infection cases in Chelyabinsk"

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  7  days in TBATS Model for  ==>  Forecasting Cumulative Covid 19 infection cases in Chelyabinsk"

paste(MAPE_Mean_All.TBATS,"%")

## [1] "0.041 % MAPE  7 days Forecasting Cumulative Covid 19 infection cases in Chelyabinsk %"

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  7  days in TBATS Model for  ==>  Forecasting Cumulative Covid 19 infection cases in Chelyabinsk"

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-04-03 | Saturday                | 55783.00    | 55783.40          | 0.001 %          |
## +---+------------+-------------------------+-------------+-------------------+------------------+
## | 2 | 2021-04-04 | Sunday                  | 55908.00    | 55913.63          | 0.01 %           |
## +---+------------+-------------------------+-------------+-------------------+------------------+
## | 3 | 2021-04-05 | Monday                  | 56031.00    | 56043.19          | 0.022 %          |
## +---+------------+-------------------------+-------------+-------------------+------------------+
## | 4 | 2021-04-06 | Tuesday                 | 56153.00    | 56173.80          | 0.037 %          |
## +---+------------+-------------------------+-------------+-------------------+------------------+
## | 5 | 2021-04-07 | Wednesday               | 56274.00    | 56305.35          | 0.056 %          |
## +---+------------+-------------------------+-------------+-------------------+------------------+
## | 6 | 2021-04-08 | Thursday                | 56394.00    | 56434.04          | 0.071 %          |
## +---+------------+-------------------------+-------------+-------------------+------------------+
## | 7 | 2021-04-09 | Friday                  | 56513.00    | 56562.46          | 0.088 %          |
## +---+------------+-------------------------+-------------+-------------------+------------------+

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-04-10 | Saturday        | 56692.70             | 56652.47  | 56631.18  | 56732.93  | 56754.22  |
## +----+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 2  | 2021-04-11 | Sunday          | 56822.26             | 56779.63  | 56757.06  | 56864.89  | 56887.45  |
## +----+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 3  | 2021-04-12 | Monday          | 56952.87             | 56907.96  | 56884.19  | 56997.77  | 57021.54  |
## +----+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 4  | 2021-04-13 | Tuesday         | 57084.42             | 57037.35  | 57012.44  | 57131.48  | 57156.40  |
## +----+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 5  | 2021-04-14 | Wednesday       | 57213.11             | 57163.98  | 57137.97  | 57262.24  | 57288.25  |
## +----+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 6  | 2021-04-15 | Thursday        | 57341.53             | 57290.43  | 57263.38  | 57392.63  | 57419.68  |
## +----+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 7  | 2021-04-16 | Friday          | 57471.77             | 57418.78  | 57390.72  | 57524.76  | 57552.81  |
## +----+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 8  | 2021-04-17 | Saturday        | 57601.33             | 57546.51  | 57517.49  | 57656.14  | 57685.16  |
## +----+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 9  | 2021-04-18 | Sunday          | 57731.93             | 57675.36  | 57645.40  | 57788.51  | 57818.47  |
## +----+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 10 | 2021-04-19 | Monday          | 57863.49             | 57805.20  | 57774.34  | 57921.78  | 57952.63  |
## +----+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 11 | 2021-04-20 | Tuesday         | 57992.18             | 57932.23  | 57900.50  | 58052.12  | 58083.86  |
## +----+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 12 | 2021-04-21 | Wednesday       | 58120.60             | 58059.05  | 58026.47  | 58182.15  | 58214.74  |
## +----+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 13 | 2021-04-22 | Thursday        | 58250.84             | 58187.73  | 58154.32  | 58313.95  | 58347.35  |
## +----+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 14 | 2021-04-23 | Friday          | 58380.39             | 58315.77  | 58281.55  | 58445.02  | 58479.24  |
## +----+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 15 | 2021-04-24 | Saturday        | 58511.00             | 58444.89  | 58409.89  | 58577.12  | 58612.12  |
## +----+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 16 | 2021-04-25 | Sunday          | 58642.55             | 58574.99  | 58539.22  | 58710.12  | 58745.89  |
## +----+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 17 | 2021-04-26 | Monday          | 58771.25             | 58702.26  | 58665.75  | 58840.23  | 58876.74  |
## +----+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 18 | 2021-04-27 | Tuesday         | 58899.67             | 58829.31  | 58792.06  | 58970.03  | 59007.28  |
## +----+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 19 | 2021-04-28 | Wednesday       | 59029.90             | 58958.19  | 58920.23  | 59101.62  | 59139.58  |
## +----+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 20 | 2021-04-29 | Thursday        | 59159.46             | 59086.43  | 59047.76  | 59232.50  | 59271.16  |
## +----+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 21 | 2021-04-30 | Friday          | 59290.07             | 59215.73  | 59176.38  | 59364.41  | 59403.76  |
## +----+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+

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 -0.5620706 11.73011 7.066356 Inf  Inf 0.04900932 0.1667545

# 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.6407 
## 
##   Smoothing parameters:
##     alpha = 0.9999 
##     beta  = 0.6195 
## 
##   Initial states:
##     l = -1.503 
##     b = -0.1886 
## 
##   sigma:  0.885
## 
##      AIC     AICc      BIC 
## 2217.329 2217.487 2237.121 
## 
## Training set error measures:
##                      ME     RMSE      MAE MPE MAPE       MASE      ACF1
## Training set -0.5620706 11.73011 7.066356 Inf  Inf 0.04900932 0.1667545

# 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  7 days by using holt Model for  ==>  Forecasting Cumulative Covid 19 infection cases in Chelyabinsk"

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  7  days in holt Model for  ==>  Forecasting Cumulative Covid 19 infection cases in Chelyabinsk"

paste(MAPE_Mean_All.Holt,"%")

## [1] "0.043 % MAPE  7 days Forecasting Cumulative Covid 19 infection cases in Chelyabinsk %"

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  7  days in holt Model for  ==>  Forecasting Cumulative Covid 19 infection cases in Chelyabinsk"

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-04-03 | Saturday                | 55783.00    | 55784.84         | 0.003 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 2 | 2021-04-04 | Sunday                  | 55908.00    | 55914.80         | 0.012 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 3 | 2021-04-05 | Monday                  | 56031.00    | 56044.86         | 0.025 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 4 | 2021-04-06 | Tuesday                 | 56153.00    | 56175.03         | 0.039 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 5 | 2021-04-07 | Wednesday               | 56274.00    | 56305.30         | 0.056 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 6 | 2021-04-08 | Thursday                | 56394.00    | 56435.69         | 0.074 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+
## | 7 | 2021-04-09 | Friday                  | 56513.00    | 56566.18         | 0.094 %         |
## +---+------------+-------------------------+-------------+------------------+-----------------+

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-04-10 | Saturday        | 56696.79            | 56129.25  | 55829.65  | 57266.37  | 57568.72  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 2  | 2021-04-11 | Sunday          | 56827.50            | 56163.37  | 55812.94  | 57494.42  | 57848.60  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 3  | 2021-04-12 | Monday          | 56958.32            | 56192.65  | 55788.84  | 57727.70  | 58136.48  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 4  | 2021-04-13 | Tuesday         | 57089.24            | 56217.30  | 55757.67  | 57966.00  | 58432.07  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 5  | 2021-04-14 | Wednesday       | 57220.28            | 56237.51  | 55719.74  | 58209.15  | 58735.09  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 6  | 2021-04-15 | Thursday        | 57351.42            | 56253.45  | 55675.32  | 58456.99  | 59045.31  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 7  | 2021-04-16 | Friday          | 57482.67            | 56265.28  | 55624.63  | 58709.40  | 59362.55  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 8  | 2021-04-17 | Saturday        | 57614.03            | 56273.13  | 55567.90  | 58966.24  | 59686.61  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 9  | 2021-04-18 | Sunday          | 57745.50            | 56277.12  | 55505.32  | 59227.41  | 60017.34  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 10 | 2021-04-19 | Monday          | 57877.07            | 56277.37  | 55437.08  | 59492.82  | 60354.59  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 11 | 2021-04-20 | Tuesday         | 58008.75            | 56273.99  | 55363.34  | 59762.36  | 60698.24  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 12 | 2021-04-21 | Wednesday       | 58140.54            | 56267.06  | 55284.25  | 60035.97  | 61048.17  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 13 | 2021-04-22 | Thursday        | 58272.44            | 56256.68  | 55199.96  | 60313.57  | 61404.26  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 14 | 2021-04-23 | Friday          | 58404.44            | 56242.93  | 55110.60  | 60595.08  | 61766.43  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 15 | 2021-04-24 | Saturday        | 58536.55            | 56225.89  | 55016.29  | 60880.46  | 62134.58  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 16 | 2021-04-25 | Sunday          | 58668.77            | 56205.64  | 54917.16  | 61169.63  | 62508.63  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 17 | 2021-04-26 | Monday          | 58801.09            | 56182.23  | 54813.32  | 61462.56  | 62888.51  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 18 | 2021-04-27 | Tuesday         | 58933.53            | 56155.74  | 54704.87  | 61759.18  | 63274.15  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 19 | 2021-04-28 | Wednesday       | 59066.07            | 56126.23  | 54591.91  | 62059.45  | 63665.47  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 20 | 2021-04-29 | Thursday        | 59198.71            | 56093.75  | 54474.54  | 62363.33  | 64062.43  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 21 | 2021-04-30 | Friday          | 59331.47            | 56058.36  | 54352.84  | 62670.78  | 64464.96  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+

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  ==>  Forecasting Cumulative Covid 19 infection cases in Chelyabinsk"

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 = 6.1522, 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.1993, Truncation lag parameter = 5, p-value
## = 0.9842
## alternative hypothesis: stationary

adf.test(data_series)  # applay adf test

## 
##  Augmented Dickey-Fuller Test
## 
## data:  data_series
## Dickey-Fuller = -3.4538, Lag order = 7, p-value = 0.04719
## 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  ==>  Forecasting Cumulative Covid 19 infection cases in Chelyabinsk"

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.9696, Truncation lag parameter = 5, 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) = -1.8717, Truncation lag parameter = 5, p-value
## = 0.9724
## alternative hypothesis: stationary

adf.test(diff1_x1)    # applay adf test after taking first differences

## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff1_x1
## Dickey-Fuller = -0.49021, Lag order = 7, p-value = 0.9821
## 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 Forecasting Cumulative Covid 19 infection cases in Chelyabinsk"

kpss.test(diff2_x1)   # applay kpss test after taking Second differences

## 
##  KPSS Test for Level Stationarity
## 
## data:  diff2_x1
## KPSS Level = 0.35539, Truncation lag parameter = 5, p-value = 0.09638

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) = -435.74, Truncation lag parameter = 5, p-value
## = 0.01
## alternative hypothesis: stationary

adf.test(diff2_x1)    # applay adf test after taking Second differences

## Warning in adf.test(diff2_x1): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff2_x1
## Dickey-Fuller = -6.6133, Lag order = 7, p-value = 0.01
## 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)                    : 2988.391
##  ARIMA(0,2,1)                    : 2974.147
##  ARIMA(0,2,2)                    : 2976.106
##  ARIMA(0,2,3)                    : 2973.714
##  ARIMA(0,2,4)                    : 2974.144
##  ARIMA(0,2,5)                    : 2967.241
##  ARIMA(1,2,0)                    : 2973.806
##  ARIMA(1,2,1)                    : 2975.831
##  ARIMA(1,2,2)                    : 2977.821
##  ARIMA(1,2,3)                    : 2975.381
##  ARIMA(1,2,4)                    : 2970.371
##  ARIMA(2,2,0)                    : 2975.833
##  ARIMA(2,2,1)                    : Inf
##  ARIMA(2,2,2)                    : 2978.748
##  ARIMA(2,2,3)                    : 2967.901
##  ARIMA(3,2,0)                    : 2977.293
##  ARIMA(3,2,1)                    : 2977.207
##  ARIMA(3,2,2)                    : 2970.983
##  ARIMA(4,2,0)                    : 2967.396
##  ARIMA(4,2,1)                    : 2964.695
##  ARIMA(5,2,0)                    : 2965.317
## 
## 
## 
##  Best model: ARIMA(4,2,1)

model1 # show the result of autoarima 

## Series: data_series 
## ARIMA(4,2,1) 
## 
## Coefficients:
##           ar1      ar2      ar3      ar4     ma1
##       -0.6737  -0.0999  -0.0780  -0.1985  0.4827
## s.e.   0.1475   0.0667   0.0603   0.0507  0.1456
## 
## sigma^2 estimated as 126.9:  log likelihood=-1476.24
## AIC=2964.47   AICc=2964.7   BIC=2988.19

#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] 4 2 1

strtoi(bestmodel[3])

## [1] 1

#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      ar2      ar3      ar4     ma1
##       -0.6737  -0.0999  -0.0780  -0.1985  0.4827
## s.e.   0.1475   0.0667   0.0603   0.0507  0.1456
## 
## sigma^2 estimated as 125.2:  log likelihood = -1476.24,  aic = 2964.47

paste ("accuracy of autoarima Model For  ==> ",y_lab, sep=" ")

## [1] "accuracy of autoarima Model For  ==>  Forecasting Cumulative Covid 19 infection cases in Chelyabinsk"

accuracy(x1_model1)  # aacuracy of best model from auto arima

##                     ME     RMSE     MAE       MPE     MAPE       MASE
## Training set 0.4635585 11.16239 6.49806 0.3410782 1.864202 0.04506785
##                      ACF1
## Training set 0.0008884164

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(4,2,1)
## Q* = 9.2535, df = 5, p-value = 0.09937
## 
## 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   ==>  Forecasting Cumulative Covid 19 infection cases in Chelyabinsk"

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 = 168.18, df = 20, p-value < 2.2e-16

library(tseries)
jarque.bera.test(x1_model1$residuals)  # Do test jarque.bera.test 

## 
##  Jarque Bera Test
## 
## data:  x1_model1$residuals
## X-squared = 2300.8, 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  7 days by using bats Model for  ==>  Forecasting Cumulative Covid 19 infection cases in Chelyabinsk"

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  7  days in bats Model for  ==>  Forecasting Cumulative Covid 19 infection cases in Chelyabinsk"

paste(MAPE_Mean_All.ARIMA,"%")

## [1] "0.041 % MAPE  7 days Forecasting Cumulative Covid 19 infection cases in Chelyabinsk %"

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  7  days in bats Model for  ==>  Forecasting Cumulative Covid 19 infection cases in Chelyabinsk"

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-04-03      | Saturday                | 55783.00    | 55784.49               | 0.003 %               |
## +---+-----------------+-------------------------+-------------+------------------------+-----------------------+
## | 2 | 2021-04-04      | Sunday                  | 55908.00    | 55914.23               | 0.011 %               |
## +---+-----------------+-------------------------+-------------+------------------------+-----------------------+
## | 3 | 2021-04-05      | Monday                  | 56031.00    | 56044.03               | 0.023 %               |
## +---+-----------------+-------------------------+-------------+------------------------+-----------------------+
## | 4 | 2021-04-06      | Tuesday                 | 56153.00    | 56173.92               | 0.037 %               |
## +---+-----------------+-------------------------+-------------+------------------------+-----------------------+
## | 5 | 2021-04-07      | Wednesday               | 56274.00    | 56303.63               | 0.053 %               |
## +---+-----------------+-------------------------+-------------+------------------------+-----------------------+
## | 6 | 2021-04-08      | Thursday                | 56394.00    | 56433.40               | 0.07 %                |
## +---+-----------------+-------------------------+-------------+------------------------+-----------------------+
## | 7 | 2021-04-09      | Friday                  | 56513.00    | 56563.12               | 0.089 %               |
## +---+-----------------+-------------------------+-------------+------------------------+-----------------------+

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-04-10 | Saturday        | 56692.87                  | 56528.85  | 56442.02  | 56856.89  | 56943.71  |
## +----+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 2  | 2021-04-11 | Sunday          | 56822.64                  | 56630.71  | 56529.11  | 57014.57  | 57116.17  |
## +----+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 3  | 2021-04-12 | Monday          | 56952.38                  | 56731.21  | 56614.12  | 57173.55  | 57290.63  |
## +----+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 4  | 2021-04-13 | Tuesday         | 57082.14                  | 56830.30  | 56696.98  | 57333.99  | 57467.31  |
## +----+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 5  | 2021-04-14 | Wednesday       | 57211.89                  | 56928.17  | 56777.97  | 57495.61  | 57645.80  |
## +----+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 6  | 2021-04-15 | Thursday        | 57341.64                  | 57024.76  | 56857.01  | 57658.52  | 57826.27  |
## +----+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 7  | 2021-04-16 | Friday          | 57471.40                  | 57120.18  | 56934.25  | 57822.62  | 58008.54  |
## +----+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 8  | 2021-04-17 | Saturday        | 57601.15                  | 57214.44  | 57009.72  | 57987.86  | 58192.57  |
## +----+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 9  | 2021-04-18 | Sunday          | 57730.90                  | 57307.57  | 57083.47  | 58154.24  | 58378.34  |
## +----+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 10 | 2021-04-19 | Monday          | 57860.65                  | 57399.63  | 57155.58  | 58321.68  | 58565.73  |
## +----+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 11 | 2021-04-20 | Tuesday         | 57990.41                  | 57490.62  | 57226.05  | 58490.19  | 58754.76  |
## +----+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 12 | 2021-04-21 | Wednesday       | 58120.16                  | 57580.60  | 57294.97  | 58659.72  | 58945.35  |
## +----+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 13 | 2021-04-22 | Thursday        | 58249.91                  | 57669.57  | 57362.35  | 58830.26  | 59137.47  |
## +----+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 14 | 2021-04-23 | Friday          | 58379.66                  | 57757.56  | 57428.23  | 59001.77  | 59331.09  |
## +----+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 15 | 2021-04-24 | Saturday        | 58509.42                  | 57844.60  | 57492.66  | 59174.24  | 59526.17  |
## +----+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 16 | 2021-04-25 | Sunday          | 58639.17                  | 57930.70  | 57555.66  | 59347.64  | 59722.68  |
## +----+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 17 | 2021-04-26 | Monday          | 58768.92                  | 58015.89  | 57617.26  | 59521.96  | 59920.59  |
## +----+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 18 | 2021-04-27 | Tuesday         | 58898.68                  | 58100.18  | 57677.49  | 59697.17  | 60119.86  |
## +----+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 19 | 2021-04-28 | Wednesday       | 59028.43                  | 58183.60  | 57736.37  | 59873.26  | 60320.49  |
## +----+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 20 | 2021-04-29 | Thursday        | 59158.18                  | 58266.15  | 57793.94  | 60050.21  | 60522.43  |
## +----+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 21 | 2021-04-30 | Friday          | 59287.93                  | 58347.85  | 57850.20  | 60228.01  | 60725.66  |
## +----+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+

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.041 % MAPE  7 days Forecasting Cumulative Covid 19 infection cases in Chelyabinsk"

# SIR Model 
#install.packages("dplyr")
library(deSolve)
first<-rows-(validation_data_days+N_forecasting_days-1)
secondr<-rows-N_forecasting_days
vector_SIR<-original_data[first:secondr]
Infected <- c(vector_SIR)
Day <- 1:(length(Infected))
N <- Population # population of the us
SIR <- function(time, state, parameters) {
  par <- as.list(c(state, parameters))
  with(par, {
    dS <- -beta/N * I * S
    dI <- beta/N * I * S - gamma * I
    dR <- gamma * I
    list(c(dS, dI, dR))
  })
}

init <- c(S = N-Infected[1], I = Infected[1], R = 0)
RSS <- function(parameters) {
  names(parameters) <- c("beta", "gamma")
  out <- ode(y = init, times = Day, func = SIR, parms = parameters)
  fit <- out[ , 3]
  sum((Infected - fit)^2)
}

# optimize with some sensible conditions
Opt <- optim(c(0.5, 0.5), RSS, method = "L-BFGS-B", 
             lower = c(0, 0), upper = c(10, 10))
Opt$message

## [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

Opt_par <- setNames(Opt$par, c("beta", "gamma"))
Opt_par

##       beta      gamma 
## 0.06607126 0.05991285

# beta     gamma 
# 0.6512503 0.4920399 

out <- ode(y = init, times = Day, func = SIR, parms = Opt_par)

plot(out)
plot(out, obs=data.frame(time=Day, I=Infected))

result_SIR<-data.frame(out)
validation_forecast<-result_SIR$I
MAPE_Mean_SIR<-round(mean(abs(((testing_data-validation_forecast)/testing_data)*100)),3)


## forecasting by SIR model

Infected <- c(tail(original_data,N_forecasting_days))
Day <- 1:(length(Infected))
N <- Population # population of the us

SIR <- function(time, state, parameters) {
  par <- as.list(c(state, parameters))
  with(par, {
    dS <- -beta/N * I * S
    dI <- beta/N * I * S - gamma * I
    dR <- gamma * I
    list(c(dS, dI, dR))
  })
}

init <- c(S = N-Infected[1], I = Infected[1], R = 0)
RSS <- function(parameters) {
  names(parameters) <- c("beta", "gamma")
  out <- ode(y = init, times = Day, func = SIR, parms = parameters)
  fit <- out[ , 3]
  sum((Infected - fit)^2)
}

# optimize with some sensible conditions
Opt <- optim(c(0.5, 0.5), RSS, method = "L-BFGS-B", 
             lower = c(0, 0), upper = c(10, 10))
Opt$message

## [1] "ERROR: ABNORMAL_TERMINATION_IN_LNSRCH"

Opt_par <- setNames(Opt$par, c("beta", "gamma"))
Opt_par

##       beta      gamma 
## 0.05065253 0.04491823

# beta     gamma 
# 0.6512503 0.4920399 

out <- ode(y = init, times = Day, func = SIR, parms = Opt_par)
plot(out)
plot(out, obs=data.frame(time=Day, I=Infected))

result_SIR <-data.frame(out)
print(ascii(data.frame(FD,forecating_date=forecasting_data_by_name,forecasting_by_SIR=result_SIR$I)), type = "rest")

## 
## +----+------------+-----------------+--------------------+
## |    | FD         | forecating_date | forecasting_by_SIR |
## +====+============+=================+====================+
## | 1  | 2021-04-10 | Saturday        | 53902.00           |
## +----+------------+-----------------+--------------------+
## | 2  | 2021-04-11 | Sunday          | 54078.04           |
## +----+------------+-----------------+--------------------+
## | 3  | 2021-04-12 | Monday          | 54248.32           |
## +----+------------+-----------------+--------------------+
## | 4  | 2021-04-13 | Tuesday         | 54412.79           |
## +----+------------+-----------------+--------------------+
## | 5  | 2021-04-14 | Wednesday       | 54571.39           |
## +----+------------+-----------------+--------------------+
## | 6  | 2021-04-15 | Thursday        | 54724.05           |
## +----+------------+-----------------+--------------------+
## | 7  | 2021-04-16 | Friday          | 54870.73           |
## +----+------------+-----------------+--------------------+
## | 8  | 2021-04-17 | Saturday        | 55011.38           |
## +----+------------+-----------------+--------------------+
## | 9  | 2021-04-18 | Sunday          | 55145.93           |
## +----+------------+-----------------+--------------------+
## | 10 | 2021-04-19 | Monday          | 55274.36           |
## +----+------------+-----------------+--------------------+
## | 11 | 2021-04-20 | Tuesday         | 55396.61           |
## +----+------------+-----------------+--------------------+
## | 12 | 2021-04-21 | Wednesday       | 55512.64           |
## +----+------------+-----------------+--------------------+
## | 13 | 2021-04-22 | Thursday        | 55622.42           |
## +----+------------+-----------------+--------------------+
## | 14 | 2021-04-23 | Friday          | 55725.91           |
## +----+------------+-----------------+--------------------+
## | 15 | 2021-04-24 | Saturday        | 55823.07           |
## +----+------------+-----------------+--------------------+
## | 16 | 2021-04-25 | Sunday          | 55913.89           |
## +----+------------+-----------------+--------------------+
## | 17 | 2021-04-26 | Monday          | 55998.34           |
## +----+------------+-----------------+--------------------+
## | 18 | 2021-04-27 | Tuesday         | 56076.38           |
## +----+------------+-----------------+--------------------+
## | 19 | 2021-04-28 | Wednesday       | 56148.01           |
## +----+------------+-----------------+--------------------+
## | 20 | 2021-04-29 | Thursday        | 56213.21           |
## +----+------------+-----------------+--------------------+
## | 21 | 2021-04-30 | Friday          | 56271.95           |
## +----+------------+-----------------+--------------------+

# 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,MAPE_Mean_SIR)
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  ==>  Forecasting Cumulative Covid 19 infection cases in Chelyabinsk"

best_recommended_model

## [1] 0.037

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")}
x6<-if(best_recommended_model >= MAPE_Mean_SIR) {paste("SIR Model")}
panderOptions('table.split.table', Inf)
paste("Forecasting by using BATS Model  ==> ", y_lab , sep=" ")

## [1] "Forecasting by using BATS Model  ==>  Forecasting Cumulative Covid 19 infection cases in Chelyabinsk"

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-04-10 | Saturday        | 56688.79            | 56512.89  | 56419.77  | 56512.89  | 56419.77  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 2  | 2021-04-11 | Sunday          | 56818.01            | 56610.92  | 56501.30  | 56610.92  | 56501.30  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 3  | 2021-04-12 | Monday          | 56947.23            | 56707.32  | 56580.31  | 56707.32  | 56580.31  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 4  | 2021-04-13 | Tuesday         | 57076.45            | 56802.14  | 56656.93  | 56802.14  | 56656.93  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 5  | 2021-04-14 | Wednesday       | 57205.68            | 56895.47  | 56731.26  | 56895.47  | 56731.26  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 6  | 2021-04-15 | Thursday        | 57334.90            | 56987.36  | 56803.39  | 56987.36  | 56803.39  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 7  | 2021-04-16 | Friday          | 57464.12            | 57077.87  | 56873.40  | 57077.87  | 56873.40  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 8  | 2021-04-17 | Saturday        | 57593.34            | 57167.03  | 56941.36  | 57167.03  | 56941.36  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 9  | 2021-04-18 | Sunday          | 57722.56            | 57254.91  | 57007.34  | 57254.91  | 57007.34  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 10 | 2021-04-19 | Monday          | 57851.78            | 57341.52  | 57071.41  | 57341.52  | 57071.41  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 11 | 2021-04-20 | Tuesday         | 57981.01            | 57426.92  | 57133.61  | 57426.92  | 57133.61  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 12 | 2021-04-21 | Wednesday       | 58110.23            | 57511.13  | 57193.99  | 57511.13  | 57193.99  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 13 | 2021-04-22 | Thursday        | 58239.45            | 57594.19  | 57252.61  | 57594.19  | 57252.61  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 14 | 2021-04-23 | Friday          | 58368.67            | 57676.12  | 57309.51  | 57676.12  | 57309.51  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 15 | 2021-04-24 | Saturday        | 58497.89            | 57756.95  | 57364.72  | 57756.95  | 57364.72  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 16 | 2021-04-25 | Sunday          | 58627.11            | 57836.70  | 57418.28  | 57836.70  | 57418.28  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 17 | 2021-04-26 | Monday          | 58756.33            | 57915.40  | 57470.23  | 57915.40  | 57470.23  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 18 | 2021-04-27 | Tuesday         | 58885.56            | 57993.06  | 57520.60  | 57993.06  | 57520.60  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 19 | 2021-04-28 | Wednesday       | 59014.78            | 58069.71  | 57569.43  | 58069.71  | 57569.43  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 20 | 2021-04-29 | Thursday        | 59144.00            | 58145.37  | 57616.73  | 58145.37  | 57616.73  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 21 | 2021-04-30 | Friday          | 59273.22            | 58220.06  | 57662.55  | 58220.06  | 57662.55  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+

paste("Forecasting by using TBATS Model  ==> ", y_lab , sep=" ")

## [1] "Forecasting by using TBATS Model  ==>  Forecasting Cumulative Covid 19 infection cases in Chelyabinsk"

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-04-10 | Saturday        | 56692.70             | 56652.47  | 56631.18  | 56732.93  | 56754.22  |
## +----+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 2  | 2021-04-11 | Sunday          | 56822.26             | 56779.63  | 56757.06  | 56864.89  | 56887.45  |
## +----+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 3  | 2021-04-12 | Monday          | 56952.87             | 56907.96  | 56884.19  | 56997.77  | 57021.54  |
## +----+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 4  | 2021-04-13 | Tuesday         | 57084.42             | 57037.35  | 57012.44  | 57131.48  | 57156.40  |
## +----+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 5  | 2021-04-14 | Wednesday       | 57213.11             | 57163.98  | 57137.97  | 57262.24  | 57288.25  |
## +----+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 6  | 2021-04-15 | Thursday        | 57341.53             | 57290.43  | 57263.38  | 57392.63  | 57419.68  |
## +----+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 7  | 2021-04-16 | Friday          | 57471.77             | 57418.78  | 57390.72  | 57524.76  | 57552.81  |
## +----+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 8  | 2021-04-17 | Saturday        | 57601.33             | 57546.51  | 57517.49  | 57656.14  | 57685.16  |
## +----+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 9  | 2021-04-18 | Sunday          | 57731.93             | 57675.36  | 57645.40  | 57788.51  | 57818.47  |
## +----+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 10 | 2021-04-19 | Monday          | 57863.49             | 57805.20  | 57774.34  | 57921.78  | 57952.63  |
## +----+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 11 | 2021-04-20 | Tuesday         | 57992.18             | 57932.23  | 57900.50  | 58052.12  | 58083.86  |
## +----+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 12 | 2021-04-21 | Wednesday       | 58120.60             | 58059.05  | 58026.47  | 58182.15  | 58214.74  |
## +----+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 13 | 2021-04-22 | Thursday        | 58250.84             | 58187.73  | 58154.32  | 58313.95  | 58347.35  |
## +----+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 14 | 2021-04-23 | Friday          | 58380.39             | 58315.77  | 58281.55  | 58445.02  | 58479.24  |
## +----+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 15 | 2021-04-24 | Saturday        | 58511.00             | 58444.89  | 58409.89  | 58577.12  | 58612.12  |
## +----+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 16 | 2021-04-25 | Sunday          | 58642.55             | 58574.99  | 58539.22  | 58710.12  | 58745.89  |
## +----+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 17 | 2021-04-26 | Monday          | 58771.25             | 58702.26  | 58665.75  | 58840.23  | 58876.74  |
## +----+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 18 | 2021-04-27 | Tuesday         | 58899.67             | 58829.31  | 58792.06  | 58970.03  | 59007.28  |
## +----+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 19 | 2021-04-28 | Wednesday       | 59029.90             | 58958.19  | 58920.23  | 59101.62  | 59139.58  |
## +----+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 20 | 2021-04-29 | Thursday        | 59159.46             | 59086.43  | 59047.76  | 59232.50  | 59271.16  |
## +----+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+
## | 21 | 2021-04-30 | Friday          | 59290.07             | 59215.73  | 59176.38  | 59364.41  | 59403.76  |
## +----+------------+-----------------+----------------------+-----------+-----------+-----------+-----------+

paste("Forecasting by using Holt's Linear Trend Model  ==> ", y_lab , sep=" ")

## [1] "Forecasting by using Holt's Linear Trend Model  ==>  Forecasting Cumulative Covid 19 infection cases in Chelyabinsk"

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-04-10 | Saturday        | 56696.79            | 56129.25  | 55829.65  | 57266.37  | 57568.72  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 2  | 2021-04-11 | Sunday          | 56827.50            | 56163.37  | 55812.94  | 57494.42  | 57848.60  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 3  | 2021-04-12 | Monday          | 56958.32            | 56192.65  | 55788.84  | 57727.70  | 58136.48  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 4  | 2021-04-13 | Tuesday         | 57089.24            | 56217.30  | 55757.67  | 57966.00  | 58432.07  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 5  | 2021-04-14 | Wednesday       | 57220.28            | 56237.51  | 55719.74  | 58209.15  | 58735.09  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 6  | 2021-04-15 | Thursday        | 57351.42            | 56253.45  | 55675.32  | 58456.99  | 59045.31  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 7  | 2021-04-16 | Friday          | 57482.67            | 56265.28  | 55624.63  | 58709.40  | 59362.55  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 8  | 2021-04-17 | Saturday        | 57614.03            | 56273.13  | 55567.90  | 58966.24  | 59686.61  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 9  | 2021-04-18 | Sunday          | 57745.50            | 56277.12  | 55505.32  | 59227.41  | 60017.34  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 10 | 2021-04-19 | Monday          | 57877.07            | 56277.37  | 55437.08  | 59492.82  | 60354.59  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 11 | 2021-04-20 | Tuesday         | 58008.75            | 56273.99  | 55363.34  | 59762.36  | 60698.24  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 12 | 2021-04-21 | Wednesday       | 58140.54            | 56267.06  | 55284.25  | 60035.97  | 61048.17  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 13 | 2021-04-22 | Thursday        | 58272.44            | 56256.68  | 55199.96  | 60313.57  | 61404.26  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 14 | 2021-04-23 | Friday          | 58404.44            | 56242.93  | 55110.60  | 60595.08  | 61766.43  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 15 | 2021-04-24 | Saturday        | 58536.55            | 56225.89  | 55016.29  | 60880.46  | 62134.58  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 16 | 2021-04-25 | Sunday          | 58668.77            | 56205.64  | 54917.16  | 61169.63  | 62508.63  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 17 | 2021-04-26 | Monday          | 58801.09            | 56182.23  | 54813.32  | 61462.56  | 62888.51  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 18 | 2021-04-27 | Tuesday         | 58933.53            | 56155.74  | 54704.87  | 61759.18  | 63274.15  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 19 | 2021-04-28 | Wednesday       | 59066.07            | 56126.23  | 54591.91  | 62059.45  | 63665.47  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 20 | 2021-04-29 | Thursday        | 59198.71            | 56093.75  | 54474.54  | 62363.33  | 64062.43  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+
## | 21 | 2021-04-30 | Friday          | 59331.47            | 56058.36  | 54352.84  | 62670.78  | 64464.96  |
## +----+------------+-----------------+---------------------+-----------+-----------+-----------+-----------+

paste("Forecasting by using ARIMA Model  ==> ", y_lab , sep=" ")

## [1] "Forecasting by using ARIMA Model  ==>  Forecasting Cumulative Covid 19 infection cases in Chelyabinsk"

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-04-10 | Saturday        | 56692.87                  | 56528.85  | 56442.02  | 56856.89  | 56943.71  |
## +----+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 2  | 2021-04-11 | Sunday          | 56822.64                  | 56630.71  | 56529.11  | 57014.57  | 57116.17  |
## +----+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 3  | 2021-04-12 | Monday          | 56952.38                  | 56731.21  | 56614.12  | 57173.55  | 57290.63  |
## +----+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 4  | 2021-04-13 | Tuesday         | 57082.14                  | 56830.30  | 56696.98  | 57333.99  | 57467.31  |
## +----+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 5  | 2021-04-14 | Wednesday       | 57211.89                  | 56928.17  | 56777.97  | 57495.61  | 57645.80  |
## +----+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 6  | 2021-04-15 | Thursday        | 57341.64                  | 57024.76  | 56857.01  | 57658.52  | 57826.27  |
## +----+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 7  | 2021-04-16 | Friday          | 57471.40                  | 57120.18  | 56934.25  | 57822.62  | 58008.54  |
## +----+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 8  | 2021-04-17 | Saturday        | 57601.15                  | 57214.44  | 57009.72  | 57987.86  | 58192.57  |
## +----+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 9  | 2021-04-18 | Sunday          | 57730.90                  | 57307.57  | 57083.47  | 58154.24  | 58378.34  |
## +----+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 10 | 2021-04-19 | Monday          | 57860.65                  | 57399.63  | 57155.58  | 58321.68  | 58565.73  |
## +----+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 11 | 2021-04-20 | Tuesday         | 57990.41                  | 57490.62  | 57226.05  | 58490.19  | 58754.76  |
## +----+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 12 | 2021-04-21 | Wednesday       | 58120.16                  | 57580.60  | 57294.97  | 58659.72  | 58945.35  |
## +----+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 13 | 2021-04-22 | Thursday        | 58249.91                  | 57669.57  | 57362.35  | 58830.26  | 59137.47  |
## +----+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 14 | 2021-04-23 | Friday          | 58379.66                  | 57757.56  | 57428.23  | 59001.77  | 59331.09  |
## +----+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 15 | 2021-04-24 | Saturday        | 58509.42                  | 57844.60  | 57492.66  | 59174.24  | 59526.17  |
## +----+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 16 | 2021-04-25 | Sunday          | 58639.17                  | 57930.70  | 57555.66  | 59347.64  | 59722.68  |
## +----+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 17 | 2021-04-26 | Monday          | 58768.92                  | 58015.89  | 57617.26  | 59521.96  | 59920.59  |
## +----+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 18 | 2021-04-27 | Tuesday         | 58898.68                  | 58100.18  | 57677.49  | 59697.17  | 60119.86  |
## +----+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 19 | 2021-04-28 | Wednesday       | 59028.43                  | 58183.60  | 57736.37  | 59873.26  | 60320.49  |
## +----+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 20 | 2021-04-29 | Thursday        | 59158.18                  | 58266.15  | 57793.94  | 60050.21  | 60522.43  |
## +----+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+
## | 21 | 2021-04-30 | Friday          | 59287.93                  | 58347.85  | 57850.20  | 60228.01  | 60725.66  |
## +----+------------+-----------------+---------------------------+-----------+-----------+-----------+-----------+

paste("Forecasting by using NNAR Model  ==> ", y_lab , sep=" ")

## [1] "Forecasting by using NNAR Model  ==>  Forecasting Cumulative Covid 19 infection cases in Chelyabinsk"

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-04-10 | Saturday        | 56100.34            |
## +----+------------+-----------------+---------------------+
## | 2  | 2021-04-11 | Sunday          | 56145.15            |
## +----+------------+-----------------+---------------------+
## | 3  | 2021-04-12 | Monday          | 56187.85            |
## +----+------------+-----------------+---------------------+
## | 4  | 2021-04-13 | Tuesday         | 56228.52            |
## +----+------------+-----------------+---------------------+
## | 5  | 2021-04-14 | Wednesday       | 56267.25            |
## +----+------------+-----------------+---------------------+
## | 6  | 2021-04-15 | Thursday        | 56304.11            |
## +----+------------+-----------------+---------------------+
## | 7  | 2021-04-16 | Friday          | 56339.19            |
## +----+------------+-----------------+---------------------+
## | 8  | 2021-04-17 | Saturday        | 56372.56            |
## +----+------------+-----------------+---------------------+
## | 9  | 2021-04-18 | Sunday          | 56404.30            |
## +----+------------+-----------------+---------------------+
## | 10 | 2021-04-19 | Monday          | 56434.47            |
## +----+------------+-----------------+---------------------+
## | 11 | 2021-04-20 | Tuesday         | 56463.14            |
## +----+------------+-----------------+---------------------+
## | 12 | 2021-04-21 | Wednesday       | 56490.39            |
## +----+------------+-----------------+---------------------+
## | 13 | 2021-04-22 | Thursday        | 56516.28            |
## +----+------------+-----------------+---------------------+
## | 14 | 2021-04-23 | Friday          | 56540.87            |
## +----+------------+-----------------+---------------------+
## | 15 | 2021-04-24 | Saturday        | 56564.22            |
## +----+------------+-----------------+---------------------+
## | 16 | 2021-04-25 | Sunday          | 56586.39            |
## +----+------------+-----------------+---------------------+
## | 17 | 2021-04-26 | Monday          | 56607.44            |
## +----+------------+-----------------+---------------------+
## | 18 | 2021-04-27 | Tuesday         | 56627.41            |
## +----+------------+-----------------+---------------------+
## | 19 | 2021-04-28 | Wednesday       | 56646.36            |
## +----+------------+-----------------+---------------------+
## | 20 | 2021-04-29 | Thursday        | 56664.34            |
## +----+------------+-----------------+---------------------+
## | 21 | 2021-04-30 | Friday          | 56681.39            |
## +----+------------+-----------------+---------------------+

paste("Forecasting by using SIR Model  ==> ", y_lab , sep=" ")

## [1] "Forecasting by using SIR Model  ==>  Forecasting Cumulative Covid 19 infection cases in Chelyabinsk"

print(ascii(data.frame(FD,forecating_date=forecasting_data_by_name,forecasting_by_SIR=result_SIR$I)), type = "rest")

## 
## +----+------------+-----------------+--------------------+
## |    | FD         | forecating_date | forecasting_by_SIR |
## +====+============+=================+====================+
## | 1  | 2021-04-10 | Saturday        | 53902.00           |
## +----+------------+-----------------+--------------------+
## | 2  | 2021-04-11 | Sunday          | 54078.04           |
## +----+------------+-----------------+--------------------+
## | 3  | 2021-04-12 | Monday          | 54248.32           |
## +----+------------+-----------------+--------------------+
## | 4  | 2021-04-13 | Tuesday         | 54412.79           |
## +----+------------+-----------------+--------------------+
## | 5  | 2021-04-14 | Wednesday       | 54571.39           |
## +----+------------+-----------------+--------------------+
## | 6  | 2021-04-15 | Thursday        | 54724.05           |
## +----+------------+-----------------+--------------------+
## | 7  | 2021-04-16 | Friday          | 54870.73           |
## +----+------------+-----------------+--------------------+
## | 8  | 2021-04-17 | Saturday        | 55011.38           |
## +----+------------+-----------------+--------------------+
## | 9  | 2021-04-18 | Sunday          | 55145.93           |
## +----+------------+-----------------+--------------------+
## | 10 | 2021-04-19 | Monday          | 55274.36           |
## +----+------------+-----------------+--------------------+
## | 11 | 2021-04-20 | Tuesday         | 55396.61           |
## +----+------------+-----------------+--------------------+
## | 12 | 2021-04-21 | Wednesday       | 55512.64           |
## +----+------------+-----------------+--------------------+
## | 13 | 2021-04-22 | Thursday        | 55622.42           |
## +----+------------+-----------------+--------------------+
## | 14 | 2021-04-23 | Friday          | 55725.91           |
## +----+------------+-----------------+--------------------+
## | 15 | 2021-04-24 | Saturday        | 55823.07           |
## +----+------------+-----------------+--------------------+
## | 16 | 2021-04-25 | Sunday          | 55913.89           |
## +----+------------+-----------------+--------------------+
## | 17 | 2021-04-26 | Monday          | 55998.34           |
## +----+------------+-----------------+--------------------+
## | 18 | 2021-04-27 | Tuesday         | 56076.38           |
## +----+------------+-----------------+--------------------+
## | 19 | 2021-04-28 | Wednesday       | 56148.01           |
## +----+------------+-----------------+--------------------+
## | 20 | 2021-04-29 | Thursday        | 56213.21           |
## +----+------------+-----------------+--------------------+
## | 21 | 2021-04-30 | Friday          | 56271.95           |
## +----+------------+-----------------+--------------------+

result<-c(x1,x2,x3,x4,x5,x6)
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,SIR.Model=MAPE_Mean_SIR,Best.Model=result)

library(ascii)
print(ascii(table(table.error)), type = "rest")

## 
## +---+--------------+------------+------------+-------------+------------+-------------+-----------+------------+------+
## |   | country.name | NNAR.model | BATS.Model | TBATS.Model | Holt.Model | ARIMA.Model | SIR.Model | Best.Model | Freq |
## +===+==============+============+============+=============+============+=============+===========+============+======+
## | 1 | Chelyabinsk  | 0.458      | 0.037      | 0.041       | 0.043      | 0.041       | 4.98      | 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,MAPE_Mean_SIR)
Model<-c("NNAR.model","BATS.Model","TBATS.Model","Holt.Model","ARIMA.Model" ,"SIR 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,ARIMA Model, and SIR Model ==>",y_lab, sep=" ")

## System finished Modelling and Forecasting  by using BATS, TBATS, Holt's Linear Trend,ARIMA Model, and SIR Model ==>Forecasting Cumulative Covid 19 infection cases in Chelyabinsk

message(" Thank you for using our System For Modelling and Forecasting ==> ",y_lab, sep=" ")

##  Thank you for using our System For Modelling and Forecasting ==> Forecasting Cumulative Covid 19 infection cases in Chelyabinsk