Smaller dataset for infected cases only

infected <- covid_pl[covid_pl$Infection.Death.Recovery=="I",]

infected$Infection.Death.Recovery<- 1
#number of unique cities
cities_number <- length(unique(infected$City))
# how many infections in each province each time stamp
infected_city <- aggregate(Infection.Death.Recovery~Province+Timestamp, data = infected, FUN = sum)
# mean daily infections for each province
infected_city_mean <- aggregate(Infection.Death.Recovery~Province, data = infected_city, FUN = mean)

Preparing large dataset

#This time we will extract data (daily updated) using covid19.analytics package
data <- covid19.data("ts-confirmed")

## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
## --------------------------------------------------------------------------------

d<-growth.rate(data,geo.loc="Poland",stride=1)

## Processing...  POLAND

d2<-as.data.frame(t(d$Changes),row.names=FALSE)
dat <- tail(d2, -1) #to keep all rows except the first
zmiany<-as.numeric(dat$V1)

Time period: 01.04.2020 - 31.05.2020

Aggregated daily infections for each province of Poland presented on the box plot.
It is visible that in Śląskie there is the highest median value of daily infections as well as interquartile range. Moreover in this voivodeship the largest daily change (492 new cases) was observed.
In every province outliers over the upper quartile are present which probably indicates on the pandemic rapid growth. However, the most significant outliers are observed in Śląskie, Dolnośląskie, Łódzkie, Mazowieckie and Wielkopolskie.

ggboxplot(infected_city, x = "Province", y = "Infection.Death.Recovery", 
          color = "Province", ylab = "Number of infected", xlab = "Province") + theme(axis.text.x = element_text(angle = -45, hjust = 0)) + 
  ggtitle("Daily infections in Poland", subtitle = paste(min(infected_city$Timestamp), " -- ", max(infected_city$Timestamp)) )

Daily infections count distribution for each province.

ggpubr::gghistogram(infected_city, x = "Infection.Death.Recovery", add = "mean", rug = FALSE, fill = "Province", add_density = TRUE, bins = 100) + 
  ggtitle("Daily infections distribution in Poland", subtitle = paste(min(infected_city$Timestamp), " -- ", max(infected_city$Timestamp)) )

Time period: 23.01.2020 – today (21.06.2021)

lower <- min(zmiany)
upper <- round(max(zmiany) + 500, digits = -3)
interval <- 1000
classes <- upper/interval
limits<- cut(zmiany,seq(lower,upper,by=interval), dig.lab = 10) #setting limits of classes
table1<-table(limits) 
 kbl(table1, format = "html", col.names=c("Daily infections range", "Frequency"), caption="Infections frequency table", format.args = list(scientific = FALSE)) %>% kable_styling(bootstrap_options = c( "bordered", "condensed")) %>% row_spec(1, background = "gold") %>% row_spec(c(2,3,4,5,6,7,8,9,10,13,14,22), background = "ivory")

TAI

Tabular accuracy index is really satisfactory for 36 classes.

tab2<-classIntervals(zmiany, n = classes, style="fixed", fixedBreaks=seq(lower, upper, by=interval) )
tab2 <-jenks.tests(tab2)
kable(tab2) %>% kable_material_dark() %>% kable_styling( font_size = 24)

	x
# classes	36.0000000
Goodness of fit	0.9989104
Tabular accuracy	0.9651643

Distribution of daily infections in Poland seems extremely right skewed - obviously that’s better for the country and citizens. Below is the summary of daily infections in Poland. In fact, skewness is strong and positive which is caused by median being way lower than mean.
Moreover distribution is leptokurtic - caused by “long tail”.

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     0.0   313.0   783.5  5579.1  8617.2 35253.0

Distribution shape measures

Skewness	1.648967
Kurtosis	2.013898

Box-Cox transformation

The first transformation will be taking a logarithm of raw data because the range is wide. However, it is not possible for now since there are 0 values in the time series and operation log(0) would give an infinite result. For that reason, value of 1 is added to each observation in the time series.
Then, Box_cox transformation with lambda = 0 is applied. It is a power transformation which takes log() for a lambda = 0.
Time series looks better for the further work but still needs to be differenced - now double differencing is suggested.

zmiany_ts1 <- zmiany_ts + 1
zmiany_bc <- BoxCox(zmiany_ts1, 0)
ndiffs(zmiany_bc)

## [1] 2

tsdisplay(zmiany_bc)

First differencing

zmiany_bc_diff1 <- diff(zmiany_bc, lag = 1)
tsdisplay(zmiany_bc_diff1)

ndiffs(zmiany_bc_diff1)

## [1] 1

Acf(zmiany_bc_diff1, lag.max = 40)

On the autocorrelation plot it can be observed that observations exceed confidence interval with lag = 7 For that reason second differencing will be conducted with a lag of 7.

Second differencing

zmiany_bc_diff2 <- diff(zmiany_bc_diff1, lag = 7)
tsdisplay(zmiany_bc_diff2)

ndiffs(zmiany_bc_diff2)

## [1] 0

Acf(zmiany_bc_diff2, lag.max = 40)

No more differencing operation is required.
Testing the stationarity using the unit root test by Kwiatkowski–Phillips–Schmidt–Shin (KPSS).

library(urca)
zmiany_bc_diff1 %>% ur.kpss() %>% summary()

## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: mu with 5 lags. 
## 
## Value of test-statistic is: 1.1844 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.347 0.463  0.574 0.739

zmiany_bc_diff2 %>% ur.kpss() %>% summary()

## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: mu with 5 lags. 
## 
## Value of test-statistic is: 0.0388 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.347 0.463  0.574 0.739

Second test indicates that the statistical value is smaller than 1% critical value, which concludes that the data is stationary.

ACF, PACF

According to autocorrelation and partial autocorrelation plots, for further work models of order 7 also will be considered. It is because for the lag of 7 the largest value exceeding the confidence interval is observed.
Basing on the ACF and PACF plots analysis, appropriate ARIMA model is selected. PACF becomes insignificant after the 2nd lag, therefore the part of the model is more likely to be MA(2). ACF also becomes insignificant after 2nd lag, therefore the part of the model is more likely to be AR(2). I the arima(p,d,q), parameter ‘d’ will be set to 2 since time series was differenced twice.

Model selection

We use the principle of parsimony to decide which model is best: that is, we assume that the model with the fewest parameters is best.
For that reason three models will be considered: - ARMA(2, 0) - two parameters - ARMA(0, 2) - two parameters - ARMA(2, 2) - the most complex one with four parameters

arima0.2.2 <- Arima(zmiany_ts1, order =c(0,2,2), lambda= 0)
summary(arima0.2.2)

## Series: zmiany_ts1 
## ARIMA(0,2,2) 
## Box Cox transformation: lambda= 0 
## 
## Coefficients:
##           ma1     ma2
##       -1.6455  0.6769
## s.e.      NaN     NaN
## 
## sigma^2 estimated as 0.1203:  log likelihood=-171.56
## AIC=349.12   AICc=349.17   BIC=361.6
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE      MASE      ACF1
## Training set -69.7191 3065.353 1569.744 -8.937411 28.52469 0.1545703 0.4705795

arima2.2.0 <- Arima(zmiany_ts1, order =c(2,2,0), lambda= 0)
summary(arima2.2.0)

## Series: zmiany_ts1 
## ARIMA(2,2,0) 
## Box Cox transformation: lambda= 0 
## 
## Coefficients:
##           ar1      ar2
##       -0.7533  -0.3526
## s.e.   0.0433   0.0449
## 
## sigma^2 estimated as 0.1914:  log likelihood=-279.46
## AIC=564.93   AICc=564.98   BIC=577.4
## 
## Training set error measures:
##                     ME     RMSE      MAE       MPE     MAPE      MASE      ACF1
## Training set -385.3305 3774.596 1927.881 -9.543319 33.58852 0.1898355 0.2391352

arima2.2.2 <- Arima(zmiany_ts1, order =c(2,2,2), lambda= 0)
summary(arima2.2.2)

## Series: zmiany_ts1 
## ARIMA(2,2,2) 
## Box Cox transformation: lambda= 0 
## 
## Coefficients:
##          ar1      ar2     ma1     ma2
##       0.4477  -0.1078  -1.874  0.9070
## s.e.  0.0541   0.0525   0.032  0.0299
## 
## sigma^2 estimated as 0.1026:  log likelihood=-133.4
## AIC=276.79   AICc=276.92   BIC=297.59
## 
## Training set error measures:
##                     ME     RMSE      MAE       MPE     MAPE     MASE     ACF1
## Training set -22.50683 2471.207 1219.294 -6.797896 24.92259 0.120062 0.313234

Summary of three models indicated that model with at least 3 parameters provides the best results so far. However, during ACF and PACF plots examination it was concluded that it will be worth to consider models of order 7. For that reason additional two models are examined.

arima2.2.7 <- Arima(zmiany_ts1, order =c(2,2,7), lambda= 0)
summary(arima2.2.7)

## Series: zmiany_ts1 
## ARIMA(2,2,7) 
## Box Cox transformation: lambda= 0 
## 
## Coefficients:
##           ar1      ar2      ma1     ma2      ma3     ma4      ma5     ma6
##       -0.4417  -0.9839  -0.9295  0.8121  -1.5020  0.5391  -0.2278  0.3311
## s.e.   0.0115   0.0132   0.0509  0.0665   0.0858  0.1190   0.0898  0.0766
##          ma7
##       0.1412
## s.e.  0.0574
## 
## sigma^2 estimated as 0.08458:  log likelihood=-86.42
## AIC=192.84   AICc=193.31   BIC=234.43
## 
## Training set error measures:
##                     ME     RMSE      MAE       MPE     MAPE      MASE      ACF1
## Training set -62.17826 2342.477 1190.141 -5.336486 22.91458 0.1171913 0.3116076

arima7.2.2 <- Arima(zmiany_ts1, order =c(7,2,2), lambda= 0)
summary(arima7.2.2)

## Series: zmiany_ts1 
## ARIMA(7,2,2) 
## Box Cox transformation: lambda= 0 
## 
## Coefficients:
##           ar1      ar2      ar3      ar4      ar5      ar6      ar7      ma1
##       -1.3464  -1.2947  -1.2785  -1.2721  -1.2688  -1.0813  -0.4761  -0.2578
## s.e.   0.1069   0.0995   0.0955   0.0942   0.0943   0.0918   0.0686   0.1113
##           ma2
##       -0.2011
## s.e.   0.0966
## 
## sigma^2 estimated as 0.06886:  log likelihood=-37.36
## AIC=94.72   AICc=95.19   BIC=136.31
## 
## Training set error measures:
##                     ME     RMSE      MAE       MPE     MAPE       MASE
## Training set -90.94902 2009.586 926.9229 -4.940974 19.86228 0.09127267
##                   ACF1
## Training set 0.3593049

All Information Criteria and error measures indicate that the most appropriate model is ARIMA(7,2,2)

ARIMA(7,2,2) performs even better than automatically selected ARIMA model. It has 10 times lower AIC and BIC and forecast values are more reasonable.

Residuals

Time plot of the forecast errors shows that residuals have roughly constant variance around the mean which is close to zero.

resid7.2.2 <- arima7.2.2$residuals
checkresiduals(resid7.2.2) 

On the histogram of residuals and the quantile comparison plot it is observed that errors are almost normally distributed.

par ( mfrow =c(1 ,2) )
plotForecastErrors(resid7.2.2)
qqPlot(resid7.2.2, col.lines = "red")

## [1] 15  5

Residuals diagnostics

Shapiro - Wilk test for residuals normality:
The null hypothesis for this test is that the data are normally distributed. Hypothesis is rejected since p-value < 0.05, but generally level of normality is accepted.

shapiro.test (resid7.2.2)

## 
##  Shapiro-Wilk normality test
## 
## data:  resid7.2.2
## W = 0.93844, p-value = 4.11e-13

Ljung-Box test for residuals randomness:

Box.test(resid7.2.2, lag = 1, type = "Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  resid7.2.2
## X-squared = 0.36625, df = 1, p-value = 0.5451

Box.test(resid7.2.2, lag = 7, type = "Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  resid7.2.2
## X-squared = 13.947, df = 7, p-value = 0.05213

All p-values in the above test are greater than 0.05, therefore there is no reason to reject the null hypothesis of errors independence.