Przygotowania wstępne

Biblioteki

Na początku mojego zadania załaduję kilka bibliotek przydatnych w czasie jego realizacji. Podstawową biblioteką będzie w tym przypadku biblioteka forecast w której znajduje się bardzo dużo funkcji ogromnie ułatwiających analizę, modelowanie oraz prognozowanie szeregów czasowych. Przygotuję także dwie funkcje kable1 oraz kable2 ułatwiające formatowanie danych prezentowanych w raporcie.

library(kableExtra)
library(readxl)
library(tidyverse)
library(forecast)
library(lubridate)

kable1 = function(data, digits=3) data %>% kable(digits = digits) %>% 
  kable_styling(bootstrap_options = 
                  c("striped", "hover", "condensed", "responsive"))

kable2 = function(data, digits=3, height="320px") data %>% 
  kable1(digits = digits) %>% 
  scroll_box(width = "100%", height = height)

Wczytanie danych

Dane pobiorę z przygotowanego wcześniej arkusza Excel w którym umieściłem dane z dokumentu zawierającego treść zadania. Odczytane dane zostaną także od razu przekształcone w odpowiednie szeregi czasowe tsWM oraz tsPol.

dfZad1 <- read_excel("Zad1.xlsx", col_types = c("date", "numeric", "numeric"))
tsPol = ts(dfZad1$Pol, frequency = 7)
tsWM = ts(dfZad1$WM, frequency = 7)
dfZad1 %>% kable2

Data	Pol	WM
2021-02-13	6586	503
2021-02-14	5334	221
2021-02-15	2543	140
2021-02-16	5178	616
2021-02-17	8694	795
2021-02-18	9073	774
2021-02-19	8777	732
2021-02-20	8510	736
2021-02-21	7038	364
2021-02-22	3890	158
2021-02-23	6310	725
2021-02-24	12146	970
2021-02-25	12142	810
2021-02-26	11539	960
2021-02-27	12100	943
2021-02-28	10099	593
2021-03-01	4786	197
2021-03-02	7937	924
2021-03-03	15698	1132
2021-03-04	15250	831
2021-03-05	15829	1083
2021-03-06	14857	1129
2021-03-07	13574	706
2021-03-08	6170	202
2021-03-09	9954	835
2021-03-10	17260	1272
2021-03-11	21045	999
2021-03-12	18775	969

Warmińsko-Mazurskie

a)

Wykreślić szereg empiryczny

Do wykreślenia szeregu czasowego posłużę się bardzo wygodną funkcją autoplot z pakietu forecast. Funkcja ta tworzy wykres klasy ggplot. Tworzenie szablonowego wykresu zrealizuję w własnej funkcji plotWM. Została ona tak przygotowana aby oś X wyraźnie wskazywała na datę wraz z dniem tygodnia. W tym przypadku przyjąłem poniedziałek jako pierwszy dzień tygodnia.

Tak powstały wykres będzie szablonem w dalszych punktach zadnia.

plotWM = function(ts) suppressMessages(autoplot(ts)+
  ggtitle("Liczba zachorowań na koronawirusa w woj. Warmińsko-Mazurskim \nw dniach 13.02.2021-12.03.2021")+
  xlab("Data") +
  ylab("Liczba zachorowań")+
  theme(legend.position="bottom",
        legend.title = element_blank())+
  scale_x_continuous(
    breaks = seq(1+2/7,9,1),
    minor_breaks = NULL,
    labels = paste("Pn", seq(ymd("21-02-13"), ymd("21-04-04"), by=7))))

plotWM(tsWM)

Spoglądając na ten wykres można od razu wyciągnąć kilka ważnych wniosków. Po pierwsze ewidentnie mamy do czynienia z szeregiem z trendem oraz bardzo wyraźną, tygodniową sezonowością. Ponadto nie można pominąć faktu rosnącej wariancji w kolejnych tygodniach. To prowadzi mnie do wniosku, że dla poprawnego prognozowania, a szczególnie dla poprawnego oszacowania przedziałów predykcji należało by skorzystać z jakiejś transformacji matematycznej. Zdecydowałem, że w tym przypadku najlepiej będzie stosować transformację Boxa-Coxa¹ która zależy od parametru $\lambda$. Transformacja ta zdefiniowana jest następująco:

\[ w_t= \left\{ \begin{array}{ll} \log(y_t) & \textrm{jeżeli $\lambda =0$};\\ (y^\lambda - 1)/\lambda & \textrm{w przeciwnym przypadku}\\ \end{array} \right. \]

Do wyznaczenia odpowiedniej wartości współczynnika $\lambda$ skorzystam z funkcji Box.Cox.lambda która to zadanie wykona w pełni automatycznie.

lambda = BoxCox.lambda(tsWM)
cat("Wartość współczynnika lambda:", round(lambda, 3))

## Wartość współczynnika lambda: 0.473

W przypadku danych z województwa warmińsko-mazurskiego parametr $\lambda$ wynosi 0.473. Zobaczmy teraz jak będzie wyglądał przebieg badanego szeregu czasowego po zastosowaniu tej transformacji z wyznaczonym parametrem $\lambda$.

plotWM(BoxCox(tsWM, lambda))+
  ylab("Liczba zachorowań po zastosowaniu transformacji \nBoxa-Coxa z parametrem lambda = 0.473")

Z powyższego wykresu można odczytać, że zastosowanie transformacji Boxa-Coxa doprowadziło do wyraźnego wyrównania wariancji w poszczególnych tygodniach.

b) - c)

Wyznaczyć tendencję rozwojową za pomocą średniej ruchomej 7-dniowej
Zestawić obok siebie dane z zadania oraz dane ze średnich ruchomych 7-dniowych Zastanów się, czy różnice są duże i postaraj się wyjaśnić dlaczego. Dane zaokrąglić do całości

Do wyznaczenia tendencji rozwojowej za pomocą średniej ruchomej użyję kolejnej, bardzo przydatnej funkcji z pakietu forecast. Jest nią funkcja ma. Dane od razu zestawię w jednej tabeli wraz z różnicami względnymi.

dfZad1 %>% select(-Pol) %>% 
  mutate(
    WM.ma7 = round(ma(tsWM,7)),
    WM.wzgl = round((WM-WM.ma7)/WM, 2)) %>% 
  kable2()

Data	WM	WM.ma7	WM.wzgl
2021-02-13	503	NA	NA
2021-02-14	221	NA	NA
2021-02-15	140	NA	NA
2021-02-16	616	540	0.12
2021-02-17	795	573	0.28
2021-02-18	774	594	0.23
2021-02-19	732	596	0.19
2021-02-20	736	612	0.17
2021-02-21	364	637	-0.75
2021-02-22	158	642	-3.06
2021-02-23	725	675	0.07
2021-02-24	970	704	0.27
2021-02-25	810	737	0.09
2021-02-26	960	743	0.23
2021-02-27	943	771	0.18
2021-02-28	593	794	-0.34
2021-03-01	197	797	-3.05
2021-03-02	924	815	0.12
2021-03-03	1132	841	0.26
2021-03-04	831	857	-0.03
2021-03-05	1083	858	0.21
2021-03-06	1129	845	0.25
2021-03-07	706	865	-0.23
2021-03-08	202	889	-3.40
2021-03-09	835	873	-0.05
2021-03-10	1272	NA	NA
2021-03-11	999	NA	NA
2021-03-12	969	NA	NA

Oczywiście bardzo duże różnice (dochodzące w wartościach absolutnych do 340%) pomiędzy wartościami szeregu czasowego a wartościami średniej ruchomej wynikają z dużych wahań sezonowych, co nie było specjalnym zaskoczeniem.

d)

Dane przedstawić na jednym zbiorczym wykresie (średnia ruchoma 7-dniowa na czerwono)

Do uzupełnienia pierwotnego wykresu o kolejne serie wykorzystam funkcję autolayer pozwalającą na wygodne dodawanie kolejnych serii do pierwotnego wykresu bez konieczności zestawiania ich w jedną tabelę danych.

plotWM(tsWM) + 
  autolayer(ma(tsWM, 7), series="Średnia ruchoma 7 dniowa")

e)

Wyznaczyć trend liniowy

Kolejna przydatna funkcja z pakietu forecast będzie pomocna do wyznaczenia trendu liniowego. Jest nią funkcja tslm. Funkcja ta, w przeciwieństwie do standardowej funkcji lm jest zoptymalizowana do używania szeregów czasowych a ponadto ma tą zaletę, że można przy jej użyciu stosować opisaną wcześniej transformację Boxa-Coxa. Ponadto w nieco odmienny sposób rozwiązuje problem sezonowości, co będzie opisane w dalszej części zadania. Model z tej funkcji zapisany zostanie do zmiennej lmtWM.

lmtWM = tslm(tsWM ~ trend)
summary(lmtWM)

## 
## Call:
## tslm(formula = tsWM ~ trend)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -709.50 -168.00   81.74  218.48  321.38 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  442.056    110.664   3.995 0.000474 ***
## trend         19.560      6.667   2.934 0.006908 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 285 on 26 degrees of freedom
## Multiple R-squared:  0.2487, Adjusted R-squared:  0.2198 
## F-statistic: 8.607 on 1 and 26 DF,  p-value: 0.006908

f)

Zestawić obok siebie dane z zadania i teoretyczne uzyskane z trendu

dfZad1 %>% select(-Pol) %>% 
  mutate(
    WM.lmt = round(fitted(lmtWM)),
    dWM.wzgl = round((WM-WM.lmt)/WM, 2)) %>% 
  kable2()

Data	WM	WM.lmt	dWM.wzgl
2021-02-13	503	462	0.08
2021-02-14	221	481	-1.18
2021-02-15	140	501	-2.58
2021-02-16	616	520	0.16
2021-02-17	795	540	0.32
2021-02-18	774	559	0.28
2021-02-19	732	579	0.21
2021-02-20	736	599	0.19
2021-02-21	364	618	-0.70
2021-02-22	158	638	-3.04
2021-02-23	725	657	0.09
2021-02-24	970	677	0.30
2021-02-25	810	696	0.14
2021-02-26	960	716	0.25
2021-02-27	943	735	0.22
2021-02-28	593	755	-0.27
2021-03-01	197	775	-2.93
2021-03-02	924	794	0.14
2021-03-03	1132	814	0.28
2021-03-04	831	833	0.00
2021-03-05	1083	853	0.21
2021-03-06	1129	872	0.23
2021-03-07	706	892	-0.26
2021-03-08	202	912	-3.51
2021-03-09	835	931	-0.11
2021-03-10	1272	951	0.25
2021-03-11	999	970	0.03
2021-03-12	969	990	-0.02

Wyniki są zgodne z oczekiwaniami. Różnice względne są na podobnym poziomie jak dla średniej ruchomej.

g)

Sporządzić wykres z trendem

plotWM(tsWM) + 
  autolayer(ma(tsWM, 7), series="Średnia ruchoma 7 dniowa", )+
  autolayer(fitted(lmtWM), series="Trend liniowy", )

To czego należało się spodziewać to dość dobre pokrycie siedmiodniowej średniej ruchomej trendem liniowym. Oczywiście, średnia ruchoma ma tutaj wartość jedynie “wizualną”. Nie nadaje się ona bowiem ani do modelowania ani do prognozowania. A już na pewno nie można szacować żadnych współczynników ani statystyk na jej podstawie.

h)

Obliczyć odchylenie standardowe składnika resztowego

Wartość odchylenia standardowego rezyduów można oczywiście odczytać bezpośrednio z podsumowania modelu lub też wyekstrahować z modelu funkcją sigma. Można go także wyznaczyć korzystając z poniższego równania

\[ s^2 \left (z_i\right) = \frac{\sum_{t=1}^n \left (y_t-\bar{y_t}\right)^2}{\left (n-k\right)}= \frac{\sum_{i=1}^n z_i^2}{n-k}\ \]

n=length(tsWM)
s2.zt = 1/(n-2)*sum(residuals(lmtWM)^2)
cat("Odchylenie standardowe rezyduów wyliczone samodzielnie:", round(sqrt(s2.zt)))

## Odchylenie standardowe rezyduów wyliczone samodzielnie: 285

cat("\nOdchylenie standardowe rezyduów wyekstrahowane z modelu:", round(sigma(lmtWM)))

## 
## Odchylenie standardowe rezyduów wyekstrahowane z modelu: 285

W obu przypadkach wynik jest oczywiście taki sam co tylko potwierdziło poprawność przeprowadzonych “ręcznie” obliczeń.

i)

Obliczyć średni błąd szacunku parametrów $\alpha_0$ i $\alpha_1$ liniowej funkcji trendu

Również i w tym przypadku wartości te można odczytać bezpośrednio z podsumowania modelu, jak i obliczyć korzystając ze znanych równań

\[ D({\alpha_0})=\sqrt{\frac{s^2 \left (z_i\right) \sum_{t=1}^n t^2}{n(\sum_{t=1}^n t^2 -n \bar{t}^2 )}} \]

\[ D({\alpha_1})=\sqrt{\frac{s^2\left(z_i\right)}{\sum_{t=1}^n t^2 -n \bar{t}^2 }} \]

t=1:n
D.a0 = sqrt((s2.zt*sum(t^2))/(n*(sum(t^2)-n*mean(t)^2)))
D.a1 = sqrt(s2.zt/(sum(t^2)-n*mean(t)^2))
cat("Błąd standardowy współczynników regresji wyliczony:", 
    round(c(D.a0, D.a1), 3))

## Błąd standardowy współczynników regresji wyliczony: 110.664 6.667

cat("\nBłąd standardowy współczynników regresji z modelu:", 
    round(coef(summary(lmtWM))[, "Std. Error"], 3))

## 
## Błąd standardowy współczynników regresji z modelu: 110.664 6.667

j)

Obliczyć współczynnik zbieżności i współczynnik determinacji

Współczynnika zbieżności nie można odczytać wprost z podsumowania modelu. Można go jednak wyliczyć odczytując współczynnik determinacji. Możemy go także wyliczyć we własnym zakresie.

Współczynnik zbieżności liczy się korzystając z poniższego równania

\[ \varphi^2=\frac{\sum_{t=1}^n \left (y_t-\hat{y_t}\right)^2}{\sum_{t=1}^n \left (y_t-\bar{y_t}\right)^2} = \frac{(n-2)s^2 \left (z_i\right)}{\sum_{t=1}^n \left (y_t-\bar{y_t}\right)^2} \]

Natomiast współczynnik determinacji można wyznaczyć następująco

\[ R^2 = 1- \varphi^2 \]

phi2 = (n-2)*s2.zt/(sum((tsWM-mean(tsWM))^2))
cat("Współczynnik zbieżności wyliczony:", round(phi2, 3))

## Współczynnik zbieżności wyliczony: 0.751

cat("\nWspółczynnik zbieżności z modelu:", round(1-summary(lmtWM)$r.squared, 3))

## 
## Współczynnik zbieżności z modelu: 0.751

cat("\nWspółczynnik determinacji wyliczony:", round(1 - phi2, 3))

## 
## Współczynnik determinacji wyliczony: 0.249

cat("\nWspółczynnik determinacji z modelu:", round(summary(lmtWM)$r.squared, 3))

## 
## Współczynnik determinacji z modelu: 0.249

Bardzo niska wartość współczynnika determinacji nie jest żadnym zaskoczeniem. Przyczyną są tutaj bardzo duże wahania sezonowe.

k)

Zweryfikować hipotezę o istotności współczynników $\alpha_0$ i $\alpha_1$ trendu liniowego

Do weryfikacji hipotezy o istotności współczynników $\alpha_0$ i $\alpha_1$ posłużymy się statystyką $t$

\[ t=\frac{\alpha}{D(\alpha)} \]

dla której obszar krytyczny wynosi

\[ K=(-\infty; -t_{\alpha,n-2})\cup(t_{\alpha,n-2};\infty) \]

alpha = 0.05
K = qt(1-alpha/2, n-2)
cat("Obszar krytyczny: ", round(K, 3))

## Obszar krytyczny:  2.056

t = round(coefficients(lmtWM)/c(D.a0, D.a1), 3)
cat("\nWartość statystyki t a0 i a1:", t)

## 
## Wartość statystyki t a0 i a1: 3.995 2.934

W obu przypadkach wartość statystyki $t$ znalazła się w obszarze krytycznym więc na poziomie istotności $\alpha=0.05$ możemy odrzucić hipotezę zerową mówiącą o zerowej wartości współczynników regresji.

l) - o)

Obliczyć wskaźniki sezonowości
Sprawdzić, ile wynosi suma wskaźników sezonowości
Sporządzić wykres wskaźników sezonowości
Obliczyć absolutne poziomy wahań sezonowych

Aby wygodnie obliczać wskaźniki sezonowości opracowałem własną funkcję seasind. Funkcja ta automatycznie rozpoznaje jaki jest okres szeregu czasowego. Potrafi także rozpoznać szeregi które nie rozpoczynają się w pierwszym sezonie lub nie kończą się na ostatnim sezonie odpowiednio ograniczając dane uwzględniane w obliczeniach. Funkcja ta zwraca listę z wszystkimi obliczonymi wartościami wraz z absolutnymi poziomami wahań sezonowych.

Moja funkcja wyznacza wskaźniki sezonowości według następującego równania

\[ S_i=\frac{\bar{y_i}d}{\sum_{i=1}^d \bar{y_i}} \]

gdzie

$S_i$ - wskaźnik sezonowości dla i-tego okresu (Si),

$\bar{y_i}$ - średnia arytmetyczna wielkości badanego zjawiska w jednoimiennych podokresach (mean.yi),

$d$ - liczba podokresów (d),

Ponadto zwraca początkowy (start) oraz końcowy (end) indeks danych branych pod uwagę w obliczeniach, co jest istotne w przypadku szeregów rozpoczynających się nie w pierwszym lub kończących nie na ostatnim okresie sezonowym.

seasind = function(y){
  if(class(y) != "ts") {
    stop("argument y nie jest szeregiem czasowym")}
  out = list()
  out[["y"]] = y
  out[["d"]] = frequency(y)
  out[["n"]] = length(y)
  out[["start"]] = ifelse(start(y)[2]==1, 1, out$d-start(y)[2]+2)
  out[["end"]] = ifelse(end(y)[2]==out$d, out$n, out$n-end(y)[2])
  if((out$end - out$start + 1) < out$d) {
    stop("argument y jest zbyt krótki aby wyznaczyć wskaźniki sezonowości")}
  out[["mean.yi"]] = vector("double", out$d)
  for(i in sequence(out$d)) {
    out$mean.yi[i] = mean(y[seq(out$start, out$end, out$d)+(i-1)])}
  out[["Si"]] = out$mean.yi*out$d/sum(out$mean.yi)
  out[["sum.Si"]] = sum(out$Si)
  out[["gi"]] = mean(y)*(out$Si-1)
  out
}

Funkcja ta zwraca również sumę wskaźników sezonowości (sum.Si) daną równaniem

\[ \sum_{i=1}^d S_i \]

oraz absolutne poziomy wahań sezonowych (gi) obliczanych następująco

\[ g_i = \bar{y}(S_i -1) \]

Sprawdźmy jak funkcja ta poradzi sobie z kilkoma testowymi szeregami czasowymi.

Najpierw sprawdźmy ją z prostym szeregiem z sezonowością miesięczną.

seasind(ts(1:24, frequency = 12, start = c(2012, 1)))

## $y
##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 2012   1   2   3   4   5   6   7   8   9  10  11  12
## 2013  13  14  15  16  17  18  19  20  21  22  23  24
## 
## $d
## [1] 12
## 
## $n
## [1] 24
## 
## $start
## [1] 1
## 
## $end
## [1] 24
## 
## $mean.yi
##  [1]  7  8  9 10 11 12 13 14 15 16 17 18
## 
## $Si
##  [1] 0.56 0.64 0.72 0.80 0.88 0.96 1.04 1.12 1.20 1.28 1.36 1.44
## 
## $sum.Si
## [1] 12
## 
## $gi
##  [1] -5.5 -4.5 -3.5 -2.5 -1.5 -0.5  0.5  1.5  2.5  3.5  4.5  5.5

Teraz sprawdźmy ją z szeregiem z sezonowością kwartalną.

seasind(ts(1:12, frequency = 4, start = c(2012, 1)))

## $y
##      Qtr1 Qtr2 Qtr3 Qtr4
## 2012    1    2    3    4
## 2013    5    6    7    8
## 2014    9   10   11   12
## 
## $d
## [1] 4
## 
## $n
## [1] 12
## 
## $start
## [1] 1
## 
## $end
## [1] 12
## 
## $mean.yi
## [1] 5 6 7 8
## 
## $Si
## [1] 0.7692308 0.9230769 1.0769231 1.2307692
## 
## $sum.Si
## [1] 4
## 
## $gi
## [1] -1.5 -0.5  0.5  1.5

Na koniec zobaczmy jak funkcja ta poradzi sobie z szeregiem nie rozpoczynającym na pierwszym sezonie i nie kończącym się na sezonie ostatnim.

seasind(ts(1:20, frequency = 12, start = c(2012, 8)))

## $y
##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 2012                               1   2   3   4   5
## 2013   6   7   8   9  10  11  12  13  14  15  16  17
## 2014  18  19  20                                    
## 
## $d
## [1] 12
## 
## $n
## [1] 20
## 
## $start
## [1] 6
## 
## $end
## [1] 17
## 
## $mean.yi
##  [1]  6  7  8  9 10 11 12 13 14 15 16 17
## 
## $Si
##  [1] 0.5217391 0.6086957 0.6956522 0.7826087 0.8695652 0.9565217 1.0434783
##  [8] 1.1304348 1.2173913 1.3043478 1.3913043 1.4782609
## 
## $sum.Si
## [1] 12
## 
## $gi
##  [1] -5.0217391 -4.1086957 -3.1956522 -2.2826087 -1.3695652 -0.4565217
##  [7]  0.4565217  1.3695652  2.2826087  3.1956522  4.1086957  5.0217391

Jak można się przekonać funkcja seasind całkiem sprawnie wylicza wszelkie wskaźniki a w przypadku niepełnych okresów początkowych i/lub końcowych do obliczeń bierze dane tylko z całkowitych okresów sezonowych.

Po tym przygotowaniu mogłem już wyliczyć wskaźniki sezonowości dla badanego szeregu czasowego.

siWM = seasind(tsWM)
siWM

## $y
## Time Series:
## Start = c(1, 1) 
## End = c(4, 7) 
## Frequency = 7 
##  [1]  503  221  140  616  795  774  732  736  364  158  725  970  810  960  943
## [16]  593  197  924 1132  831 1083 1129  706  202  835 1272  999  969
## 
## $d
## [1] 7
## 
## $n
## [1] 28
## 
## $start
## [1] 1
## 
## $end
## [1] 28
## 
## $mean.yi
## [1]  827.75  471.00  174.25  775.00 1042.25  853.50  936.00
## 
## $Si
## [1] 1.1406565 0.6490477 0.2401201 1.0679659 1.4362419 1.1761406 1.2898273
## 
## $sum.Si
## [1] 7
## 
## $gi
## [1]  102.07143 -254.67857 -551.42857   49.32143  316.57143  127.82143  210.32143

Suma otrzymanych wskaźników sezonowości wynosi 7 czego należało się spodziewać w przypadku szeregu z tygodniową (siedmiodniową) sezonowością.

n)

Sporządzić wykres wskaźników sezonowości

Wykres zostanie przygotowane z użyciem funkcji pakietu ggplot.

tibble(Si = siWM$Si, d = 1:siWM$d) %>% 
  ggplot(aes(d, Si))+
  geom_point()+
  geom_line()+
  geom_hline(yintercept = 1, color="red")+
  ggtitle("Wskaźniki sezonowości zachorowań na koronawirusa \nw woj. Warmińsko-Mazurskim dniach 13.02.2021-12.03.2021")+
  ylab("Wskaźnik sezonowości")+
  xlab("Podokres d")+
  geom_text(aes(label = paste(round(Si*100), "%")), nudge_x = .4)+
  scale_x_continuous(
    breaks = 1:7,
    minor_breaks = NULL,
    labels = paste(1:7, c("So", "Ni", "Pn", "Wt", "Śr", "Cz", "Pt")))

p)

Zestawić dla poszczególnych miesięcy średnie miesięczne, wskaźniki sezonowości oraz absolutne poziomy wahań sezonowych.

Choć polecenie zadania wskazuje, iż należy zestawić średnie miesięczne wraz z wskaźnikami sezonowości trzeba zauważyć, że dla miesiąca marca mamy do dyspozycji tylko jeden pełny okres. Trudno więc będzie wyznaczyć wskaźniki sezonowości mając do dyspozycji tak ograniczone dane. Dlatego zdecydowałem aby zestawić te dane dla całego zakresu dostępnych danych, bez dzielenia ich na poszczególne miesiące.

tibble(
  d = paste(1:7, c("So", "Ni", "Pn", "Wt", "Śr", "Cz", "Pt")),
  mean = round(mean(tsWM)),
  Si = round(seasind(tsWM)$Si, 3),
  gi = round(seasind(tsWM)$gi)
) %>% kable2()

d	mean	Si	gi
1 So	726	1.141	102
2 Ni	726	0.649	-255
3 Pn	726	0.240	-551
4 Wt	726	1.068	49
5 Śr	726	1.436	317
6 Cz	726	1.176	128
7 Pt	726	1.290	210

q)

Utworzyć nowy model liniowy z uwzględnieniem wahań sezonowych

To będzie najciekawsza część zadania. Najpierw stworzę nowy model, który powstanie po dodaniu do wartości dopasowanych poprzedniego modelu poziomów wahań sezonowych.

fitt.y = fitted(lmtWM)
fitt.seas.y = fitt.y + seasind(tsWM)$gi

Niestety w przypadku takiego postępowania dodajemy wahania sezonowe niejako post factum, to znaczy po utworzeniu modelu. Z ciekawości badacza postanowiłem więc sprawdzić kilka nieco innych , alternatywnych, i jak się okazało bardzo ciekawych metod tworzenia modelu liniowego uwzględniającego wahania sezonowe.

Model lm

Na początku chciałem sprawdzić jaki otrzymam model, jeżeli do standardowej funkcji lm użyję danych z dodatkową zmienną utworzoną z wahań sezonowych.

seasWM = seasind(tsWM)
dfWM1 = tibble(
  t = 1:n,
  y = tsWM,
  d = rep(seasWM$gi, n/seasWM$d)
)

lmWM = lm(y ~ t + d, dfWM1)
summary(lmWM)

## 
## Call:
## lm(formula = y ~ t + d, data = dfWM1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -164.664  -71.718    6.858   61.908  190.610 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 500.7273    35.5745  14.075 2.19e-13 ***
## t            15.5139     2.1472   7.225 1.43e-07 ***
## d             0.9440     0.0623  15.153 4.17e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 91.07 on 25 degrees of freedom
## Multiple R-squared:  0.9262, Adjusted R-squared:  0.9203 
## F-statistic:   157 on 2 and 25 DF,  p-value: 7.052e-15

Jak widać, jest to całkiem skuteczny sposób, który dodaje wahania sezonowe niejako pre factum, to znaczy przed utworzeniem modelu. Współczynnik regresji dla zmiennej $d$ jest bardzo bliski jedności czego należało się spodziewać. Od razu także można ocenić bardzo dobre dopasowanie tego modelu o czym świadczy bardzo bliska jedności wartość współczynnika determinacji $R^2$ wynosząca w tym przypadku $0.926$.

Model tslm + season

Druga metoda którą przetestowałem polegała na skorzystaniu z funkcji tslm z pakietu forecast. W przypadku tej funkcji dodanie sezonowości robi się bardzo prosto i intuicyjnie dodając do formuły “zmienną” season. Formułą ma więc postać szereg.czasowy ~ trend + season. W tym przypadku model zbudowany jest jednak nieco inaczej i ma następującą postać

\[ y_t = \beta_0 + \beta_1t+\beta_2d_{2,t}+\dots+\beta_{k+1}d_{k,t}+\varepsilon_t \]

gdzie $d_{2,t}=1$ jeśli $t$ pochodzi z drugiego sezonu oraz 0 dla pozostałych sezonów itd. Pierwszy sezon, jako zmienna zależna liniowo od wszystkich pozostałych sezonów jest tu z oczywistych względów pominięty.

lmtsWM = tslm(tsWM ~ trend + season)
summary(lmtsWM)

## 
## Call:
## tslm(formula = tsWM ~ trend + season)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -153.90  -56.64    9.40   64.38  136.60 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   640.63      55.18  11.609 2.44e-10 ***
## trend          16.27       2.33   6.984 8.89e-07 ***
## season2      -373.02      68.27  -5.464 2.39e-05 ***
## season3      -686.04      68.38 -10.032 3.00e-09 ***
## season4      -101.56      68.58  -1.481   0.1542    
## season5       149.41      68.86   2.170   0.0422 *  
## season6       -55.61      69.21  -0.803   0.4312    
## season7        10.62      69.64   0.153   0.8803    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 96.49 on 20 degrees of freedom
## Multiple R-squared:  0.9338, Adjusted R-squared:  0.9106 
## F-statistic: 40.27 on 7 and 20 DF,  p-value: 1.977e-10

Porównując otrzymane wyniki z poprzednim modelem można zauważyć, że aktualny model jest minimalnie lepiej dopasowany o czym świadczy nieco bliższa jedności wartość współczynnika determinacji $R^2$ wtóry wynosi $0.934$. Dodatkowo można zauważyć, że bardzo istotne są tylko wahania z sezonu 2 (u nas niedziela) oraz z sezonu 3 (u nas poniedziałek). Istotne jest także wahanie dla sezonu 5 (czyli z środy) choć tylko na poziomie istotności $\alpha=0.05$.

Model tslm + lambda

Zapoznając się z pomocą dotycząca funkcji tslm zauważyłem, że pozwala ona na zastosowanie transformacji Boxa-Coxa poprzez proste podanie parametru $\lambda$. Postanowiłem więc przekonać się jak taka transformacja wpłynie na dopasowanie modelu, ale także jak później wpłynie to na przedziały predykcji dla ekstrapolowanych wartości.

lmtslWM = tslm(tsWM ~ trend + season, lambda = lambda, biasadj = T)
summary(lmtslWM)

## 
## Call:
## tslm(formula = tsWM ~ trend + season, lambda = lambda, biasadj = T)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.3369 -1.5290  0.3776  2.2539  3.5822 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  42.11665    1.75376  24.015 3.20e-16 ***
## trend         0.52378    0.07404   7.074 7.38e-07 ***
## season2     -12.75286    2.16949  -5.878 9.47e-06 ***
## season3     -27.07480    2.17328 -12.458 6.99e-11 ***
## season4      -2.73286    2.17958  -1.254    0.224    
## season5       4.08067    2.18836   1.865    0.077 .  
## season6      -1.40398    2.19961  -0.638    0.531    
## season7       0.30381    2.21328   0.137    0.892    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.066 on 20 degrees of freedom
## Multiple R-squared:  0.9491, Adjusted R-squared:  0.9312 
## F-statistic: 53.22 on 7 and 20 DF,  p-value: 1.485e-11

Okazało się, że dodanie transformacji Boxa-Coxa jeszcze bardziej poprawiło uzyskany model podwyższając współczynnik determinacji $R^2$ do wartości $0.949$. co oznacza, że w aż 95% model wyjaśnia zmienność badanej cechy.

Warto w tym miejscu zwrócić uwagę na jeszcze jedną ciekawą cechę modeli uzyskiwanych z funkcji tslm z pakietu forecast. Po pierwsze, zarówno wartości dopasowane jak i rezydua zwracane przez te modele są szeregami czasowymi. Po drugie, choć modele z zastosowaną transformacją Boxa-Coxa zwracają współczynniki regresji w wartościach transformowanych, to już wartości dopasowane są z poddawane odwrotnej transformacji aby uzyskać wartości w oryginalnej skali.

Można się o tym bardzo łatwo przekonać analizując poniższe wyniki.

cat("Wartości dopasowne modelu lmtsWM\n\n")

## Wartości dopasowne modelu lmtsWM

round(fitted(lmtsWM))

## Time Series:
## Start = c(1, 1) 
## End = c(4, 7) 
## Frequency = 7 
##    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15   16 
##  657  300    3  604  871  683  765  771  414  117  718  985  797  879  885  528 
##   17   18   19   20   21   22   23   24   25   26   27   28 
##  231  832 1099  910  993  999  642  345  946 1213 1024 1107

cat("\nWartości dopasowne modelu lmtslWM\n\n")

## 
## Wartości dopasowne modelu lmtslWM

round(fitted(lmtslWM))

## Time Series:
## Start = c(1, 1) 
## End = c(4, 7) 
## Frequency = 7 
##  [1]  636  325  103  602  835  673  744  751  407  149  713  967  791  868  876
## [16]  499  205  835 1109  919 1002 1011  600  269  967 1261 1058 1147

r)

Zestawić obok siebie dane, trend liniowy oraz trend liniowy z uwzględnieniem wahań sezonowych

Aby wykonać to polecenie zestawiłem badany szereg czsowy tsWM oraz wartości dopasowane z wszystkich pięciu modeli czyli:

wartości z modelu liniowego z trendem lmtWM,
wartości z modelu liniowego z dodanymi wahaniami sezonowymi lmt.d,
wartości z klasycznego modelu liniowego (lm) uwzględniającego wahania sezonowe lmd,
wartości z modelu z efektem sezonowości lmts (modelu tslm)
oraz na końcu wartości z modelu z dodatkowa transformacją Boxa-Coxa lmtsl.

tibble(
  tsWM = round(tsWM),
  lmt = round(fitted(lmtWM)),
  lmt.d = round(fitted(lmtWM) + seasind(tsWM)$gi),
  lmd = round(fitted(lmWM)),
  lmts = round(fitted(lmtsWM)),
  lmtsl = round(fitted(lmtslWM))
) %>% kable2()

tsWM	lmt	lmt.d	lmd	lmts	lmtsl
503	462	564	613	657	636
221	481	226	291	300	325
140	501	-51	27	3	103
616	520	570	609	604	602
795	540	856	877	871	835
774	559	687	714	683	673
732	579	789	808	765	744
736	599	701	721	771	751
364	618	363	400	414	407
158	638	86	135	117	149
725	657	707	718	718	713
970	677	993	986	985	967
810	696	824	823	797	791
960	716	926	916	879	868
943	735	838	830	885	876
593	755	500	509	528	499
197	775	223	244	231	205
924	794	843	827	832	835
1132	814	1130	1094	1099	1109
831	833	961	932	910	919
1083	853	1063	1025	993	1002
1129	872	974	938	999	1011
706	892	637	617	642	600
202	912	360	352	345	269
835	931	980	935	946	967
1272	951	1267	1203	1213	1261
999	970	1098	1040	1024	1058
969	990	1200	1134	1107	1147

s)

Zilustrować powyższe na wykresie

Jest oczywistym, że porównanie wizualne liczb jest znacznie trudniejsze od porównania wizualnego tych samych wartości zaprezentowanych na wykresie. Dlatego zgromadziłem wszystkie powyższe wyniki na jednym, zbiorczym wykresie gdzie o wiele łatwiej można je porównać.

plotWM(tsWM) + 
  autolayer(fitted(lmtWM), series="Trend liniowy")+ 
  autolayer(fitted(lmtWM) + seasind(tsWM)$gi, series="Trend liniowy + d")+ 
  autolayer(ts(fitted(lmWM), frequency = 7), series="Model lm")+ 
  autolayer(fitted(lmtsWM), series="Model lmts")+ 
  autolayer(fitted(lmtslWM), series="Model lmtsl")

Pierwszy wniosek jaki się nasuwa jest taki, że każdy model w ten czy inny sposób uwzględniający sezonowość jest znacznie lepiej dopasowany do danych niż model uwzględniający jedynie trend. Ponadto trzeba dodać, że z wszystkich modeli uwzględniających sezonowość tylko model ostatni lmtsl najlepiej oddaje zmienną wariancję szeregu, co jest oczywiście efektem zastosowania transformacji Boxa-Coxa.

t)

Obliczyć odchylenie standardowe składnika resztowego dla nowego modelu uwzględniającego wahania okresowe

Ten punkt zadania postanowiłem nieco rozszerzyć przeprowadzając dokładniejszą analizę reszt poszczególnych modeli. Pominąłem tu jednak model który powstał poprzez dodanie do dopasowanych wartości wahań sezonowych, skupiając się jedynie na tych modelach które uzyskałem bezpośrednio z funkcji lm oraz tslm.

W dokładnej analizie rezyduów bardzo pomocny będzie rozbudowany wykres, a raczej kompilacja wykresów dostępna poprzez funkcję checkresiduals również pochodzącej z pakietu forecast.

Warto zwrócić uwagę, że funkcja ta oprócz odpowiednich wykresów zwraca także wyniki z testu Ljung-Boxa, jednego z tak zwanych testów portmanteau (od francuskiego słowa opisującego walizkę lub wieszak na płaszcze, na którym znajduje się kilka elementów garderoby), opartego na następującej statystyce

\[ Q^*=T(T+2)\sum_{k=1}^h(T-k)^{-1}r^2_k \]

gdzie $h$ jest opóźnieniem a $T$ to liczba obserwacji. Duże wartości $Q^*$ sugerują, że autokorelacje w badanych rezyduach nie pochodzą z szumu białego. Dodatkowo, w celach porównawczych zbadałem zgodność rozkładu rezyduów z rozkładem normalnym posługując się znanym testem Shapiro-Wilka.

Model tslm

Pierwszym badanym modelem będzie prosty model z trendem liniowym.

checkresiduals(lmtWM, test="LB")

## 
##  Ljung-Box test
## 
## data:  Residuals from Linear regression model
## Q* = 17.923, df = 4, p-value = 0.001278
## 
## Model df: 2.   Total lags used: 6

shapiro.test(residuals(lmtWM))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(lmtWM)
## W = 0.89725, p-value = 0.009868

Mode ten jest bardzo słabo dopasowany. Świadczy o tym wyraźny wzór w resztach widoczny na pierwszym podwykresie a także brak zgodności rozkładu reszt z rozkładem normalnym, co wyraźnie widać na histogramie. O słabym dopasowaniu świadczy również wykres autokorelacji ACF. Dwa duże skoki dla opóźnień 2 oraz 7 informują nas o wyraźnej sezonowości której okres, nawet bez znajomości badanego szeregu czasowego można od razu oszacować na 7 dniowy (największy skok na wykresie ACF).

Dodatkowo można zauważyć, że p-value testu Ljung-Boxa jest bardzo niskie a co za tym idzie zmuszony jestem odrzucić hipotezę zerową o pochodzeniu reszt z szumu białego. Również dodatkowy test Shapiro-Wilka każe mi odrzucić hipotezę o pochodzeniu reszt z rozkładu normalnego.

Model lm

checkresiduals(lmWM, test="LB")

## 
##  Ljung-Box test
## 
## data:  Residuals
## Q* = 11.834, df = 3, p-value = 0.007974
## 
## Model df: 3.   Total lags used: 6

shapiro.test(residuals(lmWM))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(lmWM)
## W = 0.98415, p-value = 0.9354

Choć w tym przypadku mam znacznie lepsze dopasowanie reszt do rozkładu normalnego (p-value dla testu Shapiro-Wilka wynosi aż $0.94$), to nadal mogę obserwować znaczący skok w autokorelacji dla opóźnienia 2. Co prawda p-value testu Ljung-Boxa jest już wyższe niż w przypadku poprzedniego modelu, ale nadal jest na tyle małe, że zmuszony jestem odrzuć hipotezę zerową o pochodzeniu z szumu białego.

Model tslm + season

checkresiduals(lmtsWM, test="LB")

## 
##  Ljung-Box test
## 
## data:  Residuals from Linear regression model
## Q* = 15.295, df = 3, p-value = 0.001581
## 
## Model df: 8.   Total lags used: 11

shapiro.test(residuals(lmtsWM))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(lmtsWM)
## W = 0.96177, p-value = 0.3838

Okazuje się, że z punktu widzenia rozkładu rezyduów model ten jest minimalnie gorszy od modelu uzyskanego przy zastosowaniu standardowej funkcji lm oraz dodaniu zmiennej w postaci wahań sezonowych. Świadczy o tym nieznacznie gorsze dopasowanie rozkładu reszt do rozkładu normalnego widoczne na hitogramie oraz znacznie mniejsza wartość p-value testu Ljung-Boxa, która jest tylko nieznacznie wyższa niż uzyskana z pierwszego modelu trendu liniowego. Co prawda nie widać już skoku dla opóźnienia 7 na wykresie ACF, ale nadal pozostał znaczący skok dla opóźnienia 2.

Model tslm + lambda

checkresiduals(lmtslWM, test="LB")

## 
##  Ljung-Box test
## 
## data:  Residuals from Linear regression model
## Q* = 5.7232, df = 3, p-value = 0.1259
## 
## Model df: 8.   Total lags used: 11

shapiro.test(residuals(lmtslWM))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(lmtslWM)
## W = 0.93755, p-value = 0.09574

Dopiero po zastosowaniu transformacji Boxa-Coxa udało się uzyskać model którego rezydua nie zawierają żadnej istotnej autokorelacji, a więc są nieodróżnialne od szumu białego. Świadczy o tym bezpośrednio wartość statystki $Q^*$ testu Ljung-Boxa która jest najmniejsza dla reszt z wszystkich czterech modeli, a co za tym idzie brak jest przesłanek aby na poziomie istotności $\alpha=0.05$ odrzucić hipotezę zerową mówiąca o pochodzeniu reszt z szumu białego. Warto tu także zauważyć, że w przypadku szeregów czasowych większe znaczenie będzie miał brak znaczących autokorelacji reszt niż ich zgodność z rozkładem normalnym. Co prawda w tym przypadku możemy pozostać przy hipotezie zerowej zerowej mówiącej o pochodzeniu reszt z rozkładu normalnego jednak gdybyśmy operali się jedynie na tym teście, za najlepszy należało by uznać model drugi dla którego p-value testu Shapiro-Wilka wyniosła aż $0.94$ a mimo to rezydua posiadały znaczące korelacje poszczególnych szeregów czasowych.

t) - w)

Obliczyć odchylenie standardowe składnika resztowego dla nowego modelu uwzględniającego wahania okresowe
Porównać odchylenia standardowe składników resztkowych dla modeli (starego i nowego)
Obliczyć współczynnik zmienności i determinacji dla nowego modelu
Czy nowy model jest lepiej dopasowany?

Aby łatwo można było porównać wszystkie wartości zarówno błędów standardowych jak i współczynników zmienności i determinacji zebrałem je w przejrzystej formie w jednej tabeli.

tibble(
  model = c("tslm", "lm", "tslm + season", "tslm + lambda"),
  s2.zt = c(round(sigma(lmtWM), 3),
         round(sigma(lmWM), 3),
         round(sigma(lmtsWM), 3),
         round(sigma(lmtslWM), 3)),
  phi2 = c(round(1-summary(lmtWM)$r.squared, 3),
           round(1-summary(lmWM)$r.squared, 3),
           round(1-summary(lmtsWM)$r.squared, 3),
           round(1-summary(lmtslWM)$r.squared, 3)),
  R2 = c(round(summary(lmtWM)$r.squared, 3),
           round(summary(lmWM)$r.squared, 3),
           round(summary(lmtsWM)$r.squared, 3),
           round(summary(lmtslWM)$r.squared, 3))
) %>% kable1()

model	s2.zt	phi2	R2
tslm	284.979	0.751	0.249
lm	91.067	0.074	0.926
tslm + season	96.485	0.066	0.934
tslm + lambda	3.066	0.051	0.949

Analizując powyższą tablę jeszcze raz mogłem się przekonać, że najlepiej dopasowany jest model uwzględniający sezonowość z zastosowaniem transformacji Baxa-Coxa. W przypadku tego modelu wartość współczynnika determinacji jest najbliższa jedności. Zastanawiać nieco może bardzo niska wartość odchylenia standardowego rezyduów, jednak należy pamiętać iż ten model został skalkulowany dla wartości poddanych transformacji. Z czystej ciekawości wyliczyłem ile wynosi odchylenie liczone bezpośrednio z podanego wcześniej wzoru

s2.zt = 1/(n-8)*sum((tsWM-fitted(lmtslWM))^2)
cat("Odchylenie standardowe rezyduów wyliczone samodzielnie:", round(sqrt(s2.zt),3))

## Odchylenie standardowe rezyduów wyliczone samodzielnie: 92.704

Trzeba tu jednak wziąć pod uwagę mniejszą liczbę stopni swobody jako, że model wyznacza aż osiem współczynników regresji. Tak uzyskana wartość jest już bliższa wartościom uzyskanym dla pozostałych modeli z sezonowością choć nie jest to najniższa wartość.

Płynie stąd dla mnie ważny wniosek, że nie można bezpośrednio porównywać wartości odchylenia standardowego rezyduów uzyskanych z modeli z transformacją z wartościami uzyskanymi z modeli w których nie stosowano transformacji.

Kolejny raz przekonałem się także o dużym znaczeniu odpowiedniej diagnostyki rezyduów (wraz z analizą autokorelacji) i stosowaniu odpowiednich testów portmanteau. Wyciąganie bowiem wosków na podstawie porównywania samych tylko wartości współczynników determinacji oraz błędów standardowych rezyduów może prowadzić do niewpopranych wniosków.

x)

Ekstrapolować liczbę zakażeń na 13.03.2021 (prawdziwe wartości 13.03.2021 to 21063 dla Polski i 1115 dla warmińsko-mazurskiego oraz 17259 i 759 dla 14.03.2021)

Ostatni punkt zadania postanowiłem również nieco rozwinąć przeprowadzając prognozowanie w nieco dłuższym okresie. Do takiego prognozowania użyłem funkcji forecast z pakietu o tej samej nazwie. Funkcja ta oprócz wyznaczenia prognozy na zadany okres podaje także przedziały predykcji które bardzo łatwo można dodać do wykresu. Aby jednak nie zaciemniać wykresu każdy z modeli zostanie przedstawiony na oddzielnym wykresie.

Model tslm

dfRealWM = tibble(
  t = c(28,29)/7+1,
  y = c(1115, 759))

plotWM(tsWM) + 
  autolayer(fitted(lmtWM), series="Trend liniowy")+
  autolayer(forecast(lmtWM, h=14), series="Prognoza na kolejne 2 tyg.")+
  geom_point(aes(t, y), dfRealWM, size=3, fill='red', shape=21)

Ja można zauważyć model uwzględniający jedynie trend dał w rezultacie bardzo szerokie przedziały predykcji które odzwierciedlają dużą wariancję badanego szeregu czasowego. Rzeczywiste ilości zachorowań odbiegają znacznie od wartości prognozowanych jednak bez problemu mieszą się zarówno w 95% jak i w 80% przedziale predykcji.

Model lm

W przypadku modelu opartego na klasycznej funkcji lm wyliczenie odpowiedniej prognozy wraz z przedziałami predykcji oraz wygenerowanie wykresu będzie minimalnie trudniejsze niż było to w przypadku kiedy korzystałem z funkcji pakietu forecast. Nadal jednak można to wykonać.

newWM = tibble(
  t = 1:(n+14),
  d =  rep(seasWM$gi, 6)) %>% 
  bind_cols(predict(lmWM, newdata = ., interval = "confidence")) %>% 
  tail(14)


plotWM(tsWM) +
  geom_line(aes(t/7+6/7, fit), data=newWM, color="blue")+
  geom_ribbon(aes(x=t/7+6/7, y=fit, ymin=lwr, ymax=upr), data=newWM, alpha =.2, fill="blue")+
  autolayer(ts(fitted(lmWM), frequency = 7), series="Model liniowy lm")+
  geom_point(aes(t, y), dfRealWM, size=3, fill='red', shape=21)

Po dodaniu wahań sezonowych model ma o wile węższe przedziały predykcji. Można także zauważyć, że rzeczywista wartość zachorowań wynosząca $1115$ znalazła się dokładnie wewnątrz przedziału predykcji który w dla 13 marca wynosił $[974, 1120]$. W dniu 14 marca prognoza była jeszcze bliższa rzeczywistej liczbie zachorowań.

model tslm + season

plotWM(tsWM) + 
  autolayer(fitted(lmtsWM), series="Model tslm + season")+
  autolayer(forecast(lmtsWM, h=14), series="Prognoza na kolejne 2 tyg.")+
  geom_point(aes(t, y), dfRealWM, size=3, fill='red', shape=21)

W przypadku modelu pochodzącego z funkcji tslm obie prognozowane wartości wynoszące $1113, 756$ były bardzo bliskie rzeczywistej liczbie zachorowań występujących 13 i 14 marca 2021.

model tslm + lambda

plotWM(tsWM) + 
  autolayer(fitted(lmtslWM), series="Model tslm + lambda")+
  autolayer(forecast(lmtslWM, h=14), series="Prognoza na kolejne 2 tyg.")+
  geom_point(aes(t, y), dfRealWM, size=3, fill='red', shape=21)

W przypadku ostatniego modelu wyznaczonego z użyciem transformacji Boxa-Coxa różnica pomiędzy prognozą wynoszącą $1160$ a rzeczywistą wartością byłą nieco większa jednak i w tym przypadku obserwowana wartość rzeczywista mieści się bardzo dobrze w oszacowanym 80% przedziale predykcji. Warto także zauważyć, że tylko w tym przypadku wariancja prognozowanych wartości wyraźnie rośnie w kolejnych tygodniach co oczywiście jest efektem zastosowanej transformacji Boxa-Coxa.

Polska

Zgodnie z treścią zadani całość obliczeń została powtórzona dla danych dla danych pochodzących z całego kraju. Tym razem jednak ograniczyłem się do podstawowych wniosków pomijając wcześniejsze szczegółowe opisy stosowanych metod.

a)

Wykreślić szereg empiryczny

plotPol = function(ts) suppressMessages(autoplot(ts)+
  ggtitle("Liczba zachorowań na koronawirusa w Polsce w dniach 13.02.2021-12.03.2021")+
  xlab("Data") +
  ylab("Liczba zachorowań")+
  theme(legend.position="bottom",
        legend.title = element_blank())+
  scale_x_continuous(
    breaks = seq(1+2/7,9,1),
    minor_breaks = NULL,
    labels = paste("Pn", seq(ymd("21-02-13"), ymd("21-04-04"), by=7))))

plotPol(tsPol)

Nie jest wielkim zaskoczeniem, że w przypadku danych z całego kraju badany szereg czasowy odznacza się zarówno trendem jak i bardzo silną sezonowością. Jednak tym razem znacznie bardziej widać efekt narastającej w poszczególnych tygodniach wariancji. Należy się więc spodziewać, że tym razem transformacja Boxa-Coxa będzie miała o wiele większe znaczenie.

Podobnie jak poprzednio wartości parametru $\lambda$ wyznaczyłem całkowicie automatycznie.

lambda = BoxCox.lambda(tsPol)
cat("Wartość współczynnika lambda:", round(lambda, 3))

## Wartość współczynnika lambda: 0.038

W przypadku danych z całego kraju wartość parametru $\lambda$ wynosi 0.038 i jest znacznie niższa niż w przypadku danych z województwa warmińsko-mazurskiego.

plotPol(BoxCox(tsPol, lambda))

Tym razem zastosowanie transformacji Boxa-Coxa doprowadziło do o wiele bardziej widocznego wyrównania wariancji w poszczególnych tygodniach.

b) - c)

Wyznaczyć tendencję rozwojową za pomocą średniej ruchomej 7-dniowej
Zestawić obok siebie dane z zadania oraz dane ze średnich ruchomych 7-dniowych Zastanów się, czy różnice są duże i postaraj się wyjaśnić dlaczego. Dane zaokrąglić do całości

dfZad1 %>% select(-WM) %>% 
  mutate(
    Pol.ma7 = round(ma(tsPol,7)),
    dPol.wzgl = round((Pol-Pol.ma7)/Pol, 2)) %>% 
  kable2()

Data	Pol	Pol.ma7	dPol.wzgl
2021-02-13	6586	NA	NA
2021-02-14	5334	NA	NA
2021-02-15	2543	NA	NA
2021-02-16	5178	6598	-0.27
2021-02-17	8694	6873	0.21
2021-02-18	9073	7116	0.22
2021-02-19	8777	7309	0.17
2021-02-20	8510	7470	0.12
2021-02-21	7038	7963	-0.13
2021-02-22	3890	8402	-1.16
2021-02-23	6310	8796	-0.39
2021-02-24	12146	9309	0.23
2021-02-25	12142	9747	0.20
2021-02-26	11539	9875	0.14
2021-02-27	12100	10107	0.16
2021-02-28	10099	10614	-0.05
2021-03-01	4786	11058	-1.31
2021-03-02	7937	11671	-0.47
2021-03-03	15698	12065	0.23
2021-03-04	15250	12562	0.18
2021-03-05	15829	12759	0.19
2021-03-06	14857	13047	0.12
2021-03-07	13574	13271	0.02
2021-03-08	6170	14098	-1.28
2021-03-09	9954	14519	-0.46
2021-03-10	17260	NA	NA
2021-03-11	21045	NA	NA
2021-03-12	18775	NA	NA

Tak jak dla danych z województwa warmińsko-mazurskieko również i dla całego kraju względne różnice pomiędzy wartościami szeregu czasowego a wartościami średniej ruchomej są równie wysokie (a nawet wyższe). Jednak i w tym przypadku nie jest to żadnym zaskoczeniem i wynika wprost z silnej sezonowości szeregu.

d)

Dane przedstawić na jednym zbiorczym wykresie (średnia ruchoma 7-dniowa na czerwono)

plotPol(tsPol) + 
  autolayer(ma(tsPol, 7), series="Średnia ruchoma 7 dniowa")

e)

Wyznaczyć trend liniowy

lmtPol = tslm(tsPol ~ trend)
summary(lmtPol)

## 
## Call:
## tslm(formula = tsPol ~ trend)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##  -8233  -1431   1326   2323   5376 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4280.35    1303.65   3.283  0.00293 ** 
## trend         421.78      78.54   5.370 1.27e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3357 on 26 degrees of freedom
## Multiple R-squared:  0.5259, Adjusted R-squared:  0.5076 
## F-statistic: 28.84 on 1 and 26 DF,  p-value: 1.268e-05

f)

Zestawić obok siebie dane z zadania i teoretyczne uzyskane z trendu

dfZad1 %>% select(-WM) %>% 
  mutate(
    Pol.lmt = round(fitted(lmtPol)),
    dPol.wzgl = round((Pol-Pol.lmt)/Pol, 2)) %>% 
  kable2()

Data	Pol	Pol.lmt	dPol.wzgl
2021-02-13	6586	4702	0.29
2021-02-14	5334	5124	0.04
2021-02-15	2543	5546	-1.18
2021-02-16	5178	5967	-0.15
2021-02-17	8694	6389	0.27
2021-02-18	9073	6811	0.25
2021-02-19	8777	7233	0.18
2021-02-20	8510	7655	0.10
2021-02-21	7038	8076	-0.15
2021-02-22	3890	8498	-1.18
2021-02-23	6310	8920	-0.41
2021-02-24	12146	9342	0.23
2021-02-25	12142	9764	0.20
2021-02-26	11539	10185	0.12
2021-02-27	12100	10607	0.12
2021-02-28	10099	11029	-0.09
2021-03-01	4786	11451	-1.39
2021-03-02	7937	11872	-0.50
2021-03-03	15698	12294	0.22
2021-03-04	15250	12716	0.17
2021-03-05	15829	13138	0.17
2021-03-06	14857	13560	0.09
2021-03-07	13574	13981	-0.03
2021-03-08	6170	14403	-1.33
2021-03-09	9954	14825	-0.49
2021-03-10	17260	15247	0.12
2021-03-11	21045	15669	0.26
2021-03-12	18775	16090	0.14

Wyniki są zgodne z oczekiwaniami. Różnice względne są na podobnym poziomie jak dla średniej ruchomej.

g)

Sporządzić wykres z trendem

plotPol(tsPol) + 
  autolayer(ma(tsPol, 7), series="Średnia ruchoma 7 dniowa", )+
  autolayer(fitted(lmtPol), series="Trend liniowy", )

h)

Obliczyć odchylenie standardowe składnika resztowego

cat("\nOdchylenie standardowe rezyduów:", round(sigma(lmtPol)))

## 
## Odchylenie standardowe rezyduów: 3357

i)

Obliczyć średni błąd szacunku parametrów $\alpha_0$ i $\alpha_1$ liniowej funkcji trendu

cat("\nBłąd standardowy współczynników regresji:", 
    round(coef(summary(lmtPol))[, "Std. Error"], 3))

## 
## Błąd standardowy współczynników regresji: 1303.654 78.542

j)

Obliczyć współczynnik zbieżności i współczynnik determinacji

cat("\nWspółczynnik zbieżności:", round(1-summary(lmtPol)$r.squared, 3))

## 
## Współczynnik zbieżności: 0.474

cat("\nWspółczynnik determinacji:", round(summary(lmtPol)$r.squared, 3))

## 
## Współczynnik determinacji: 0.526

Dość zaskakujące jest, że w przypadku danych z całego kraju, mimo iż odznaczały się one podobną sezonowością jak dane w województwa warmińsko-mazurskiego, współczynnik determinacji jest już znacznie “lepszy” (ponad dwukrotnie większy).

k)

Zweryfikować hipotezę o istotności współczynników $\alpha_0$ i $\alpha_1$ trendu liniowego

alpha = 0.05
K = qt(1-alpha/2, n-2)
cat("Obszar krytyczny: ", round(K, 3))

## Obszar krytyczny:  2.056

t = round(coefficients(lmtPol)/c(D.a0, D.a1), 3)
cat("\nWartość statystyki t a0 i a1:", t)

## 
## Wartość statystyki t a0 i a1: 38.679 63.262

W obu przypadkach wartość statystyki $t$ znalazły się głęboko w obszarze krytycznym więc na poziomie istotności $\alpha=0.05$ możemy odrzucić hipotezę zerową mówiącą o zerowej wartości współczynników regresji.

l) - o)

Obliczyć wskaźniki sezonowości
Sprawdzić, ile wynosi suma wskaźników sezonowości
Sporządzić wykres wskaźników sezonowości
Obliczyć absolutne poziomy wahań sezonowych

Oczywiście w tym przypadku posłużę się tą samą, przygotowaną wcześniej funkcją seasind.

siPol = seasind(tsPol)
siPol

## $y
## Time Series:
## Start = c(1, 1) 
## End = c(4, 7) 
## Frequency = 7 
##  [1]  6586  5334  2543  5178  8694  9073  8777  8510  7038  3890  6310 12146
## [13] 12142 11539 12100 10099  4786  7937 15698 15250 15829 14857 13574  6170
## [25]  9954 17260 21045 18775
## 
## $d
## [1] 7
## 
## $n
## [1] 28
## 
## $start
## [1] 1
## 
## $end
## [1] 28
## 
## $mean.yi
## [1] 10513.25  9011.25  4347.25  7344.75 13449.50 14377.50 13730.00
## 
## $Si
## [1] 1.0112575 0.8667819 0.4181570 0.7064831 1.2936921 1.3829553 1.3206730
## 
## $sum.Si
## [1] 7
## 
## $gi
## [1]   117.0357 -1384.9643 -6048.9643 -3051.4643  3053.2857  3981.2857  3333.7857

Tu również suma wskaźników sezonowości wynosi 7 co potwierdza poprawność obliczeń.

n)

Sporządzić wykres wskaźników sezonowości

tibble(Si = siPol$Si, d = 1:siPol$d) %>% 
  ggplot(aes(d, Si))+
  geom_point()+
  geom_line()+
  geom_hline(yintercept = 1, color="red")+
  ggtitle("Wskaźniki sezonowości zachorowań na koronawirusa \nw woj. Warmińsko-Mazurskim dniach 13.02.2021-12.03.2021")+
  ylab("Wskaźnik sezonowości")+
  xlab("Podokres d")+
  geom_text(aes(label = paste(round(Si*100), "%")), nudge_x = .4)+
  scale_x_continuous(
    breaks = 1:7,
    minor_breaks = NULL,
    labels = paste(1:7, c("So", "Ni", "Pn", "Wt", "Śr", "Cz", "Pt")))

p)

Zestawić dla poszczególnych miesięcy średnie miesięczne, wskaźniki sezonowości oraz absolutne poziomy wahań sezonowych.

Tak samo jak w przypadku danych z województwa warmińsko-mazurskiego trudno będzie wyznaczyć wskaźniki sezonowości mając w marcu do dyspozycji tylko jeden pełny okres. Tu również więc zdecydowałem aby zestawić te dane dla całego zakresu dostępnych danych bez dzielenia ich na poszczególne miesiące.

tibble(
  d = paste(1:7, c("So", "Ni", "Pn", "Wt", "Śr", "Cz", "Pt")),
  mean = round(mean(tsPol)),
  Si = round(seasind(tsPol)$Si, 3),
  gi = round(seasind(tsPol)$gi)
) %>% kable2()

d	mean	Si	gi
1 So	10396	1.011	117
2 Ni	10396	0.867	-1385
3 Pn	10396	0.418	-6049
4 Wt	10396	0.706	-3051
5 Śr	10396	1.294	3053
6 Cz	10396	1.383	3981
7 Pt	10396	1.321	3334

q)

Utworzyć nowy model liniowy z uwzględnieniem wahań sezonowych

Tym razem całkowicie pominę model z dodawanymi wartościami wahań sezonowych od razu przechodząc do modeli uzyskiwanych z funkcji lm oraz tslm.

Model lm

seasPol = seasind(tsPol)
dfPol1 = tibble(
  t = 1:n,
  y = tsPol,
  d = rep(seasPol$gi, n/seasPol$d)
)

lmPol = lm(y ~ t + d, dfPol1)
summary(lmPol)

## 
## Call:
## lm(formula = y ~ t + d, data = dfPol1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2424.44 -1102.75    75.88  1001.68  2608.80 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 5.096e+03  5.139e+02   9.917 3.80e-10 ***
## t           3.655e+02  3.105e+01  11.773 1.08e-11 ***
## d           8.719e-01  7.234e-02  12.052 6.55e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1312 on 25 degrees of freedom
## Multiple R-squared:  0.9304, Adjusted R-squared:  0.9248 
## F-statistic:   167 on 2 and 25 DF,  p-value: 3.421e-15

Tym razem współczynnik regresji dla wahań sezonowych $d$ nie jest już tak bliski jedności jak w przypadku modelu dla danych z województwa warmińsko-mazurskiego, choć nie powinno to stanowić większego problemu. Model jest bardzo dobrze dopasowany o czym świadczy bardzo bliska jedności wartość współczynnika determinacji $R^2$ zbliżona do wartości uzyskanych dla modelu dla dalnych z województwa warmińsko-mazurskiego.

Model tslm + season

lmtsPol = tslm(tsPol ~ trend + season)
summary(lmtsPol)

## 
## Call:
## tslm(formula = tsPol ~ trend + season)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2173.12  -784.48   -78.46   616.92  2671.63 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6136.82     736.90   8.328 6.23e-08 ***
## trend         380.56      31.11  12.232 9.69e-11 ***
## season2     -1882.56     911.58  -2.065 0.052114 .  
## season3     -6927.12     913.17  -7.586 2.62e-07 ***
## season4     -4310.18     915.82  -4.706 0.000135 ***
## season5      1414.01     919.51   1.538 0.139771    
## season6      1961.45     924.24   2.122 0.046497 *  
## season7       933.39     929.98   1.004 0.327525    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1288 on 20 degrees of freedom
## Multiple R-squared:  0.9463, Adjusted R-squared:  0.9275 
## F-statistic: 50.33 on 7 and 20 DF,  p-value: 2.504e-11

W przypadku modelu tslm z sezonowością, nieco inaczej rozłożyła się istotność współczynników regresji dla poszczególnych sezonów. Przypomnijmy, że w poprzednio najbardziej istotne były współczynniki regresji dla sezonu 2 i 3 a nieco mniej dla sezonu 5.

Tym razem najbardziej istotne są współczynniki dla sezonów 3 i 4 i nieco mniej istotne dla sezonu 6 choć i dla 2 sezonu współczynnik jest istotny na poziomie $\alpha=0.1$.

Model tslm + lambda

Najbardziej ciekawiły mnie wyniki uzyskane z modelu z zastosowaniem transformacji Boxa-Coxa.

lmtslPol = tslm(tsPol ~ trend + season, lambda = lambda, biasadj = T)
summary(lmtslPol)

## 
## Call:
## tslm(formula = tsPol ~ trend + season, lambda = lambda, biasadj = T)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.160492 -0.040884 -0.007078  0.061303  0.119736 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 10.433583   0.047073 221.648  < 2e-16 ***
## trend        0.054095   0.001987  27.219  < 2e-16 ***
## season2     -0.291890   0.058231  -5.013 6.68e-05 ***
## season3     -1.347951   0.058333 -23.108 6.74e-16 ***
## season4     -0.645213   0.058502 -11.029 5.95e-10 ***
## season5      0.156464   0.058738   2.664  0.01491 *  
## season6      0.179274   0.059040   3.036  0.00652 ** 
## season7      0.066922   0.059407   1.127  0.27328    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0823 on 20 degrees of freedom
## Multiple R-squared:  0.9904, Adjusted R-squared:  0.987 
## F-statistic: 294.9 on 7 and 20 DF,  p-value: < 2.2e-16

No i nie zawiodłem się w najmniejszym stopniu. Jak się okazało, dodanie transformacji Boxa-Coxa dało model z współczynnikiem determinacji $R^2$ o wartości $0.99$ (sic!).

r)

Zestawić obok siebie dane, trend liniowy oraz trend liniowy z uwzględnieniem wahań sezonowych

Podobnie ja w pierwszej części zadania zestawiłem dane dopasowane uzyskane z wszystkich czterech modeli.

tibble(
  tsPol = round(tsPol),
  lmt = round(fitted(lmtPol)),
  lmd = round(fitted(lmPol)),
  lmts = round(fitted(lmtsPol)),
  lmtsl = round(fitted(lmtslPol))
) %>% kable2()

tsPol	lmt	lmd	lmts	lmtsl
6586	4702	5564	6517	6725
5334	5124	4620	5015	5673
2543	5546	919	351	2737
5178	5967	3898	3349	4756
8694	6389	9586	9454	8760
9073	6811	10761	10382	9248
8777	7233	10562	9734	8876
8510	7655	8123	9181	8796
7038	8076	7178	7679	7434
3890	8498	3477	3015	3614
6310	8920	6456	6013	6243
12146	9342	12145	12118	11427
12142	9764	13319	13046	12057
11539	10185	13120	12398	11577
12100	10607	10681	11845	11474
10099	11029	9737	10343	9713
4786	11451	6036	5679	4757
7937	11872	9015	8677	8172
15698	12294	14703	14781	14867
15250	12716	15878	15709	15678
15829	13138	15679	15062	15060
14857	13560	13240	14509	14927
13574	13981	12295	13007	12658
6170	14403	8594	8343	6245
9954	14825	11573	11341	10667
17260	15247	17262	17445	19291
21045	15669	18436	18373	20332
18775	16090	18237	17726	19538

s)

Zilustrować powyższe na wykresie

plotPol(tsPol) + 
  autolayer(fitted(lmtPol), series="Trend liniowy")+ 
  autolayer(ts(fitted(lmPol), frequency = 7), series="Model lm")+ 
  autolayer(fitted(lmtsPol), series="Model lmts")+ 
  autolayer(fitted(lmtslPol), series="Model lmtsl")

To, co przede wszystkim widać na tym wykresie to doskonałe odwzorowanie zmiennej wariancji przez model z transformacją Boxa-Coxa. Teraz jeszcze lepiej można zrozumieć dlaczego współczynnik determinacji był tak bliski jedności.

t)

Obliczyć odchylenie standardowe składnika resztowego dla nowego modelu uwzględniającego wahania okresowe

Znów, tak jak w pierwszej części zadania przeprowadzę pełną diagnostykę rezyduów.

Model tslm

checkresiduals(lmtPol, test="LB")

## 
##  Ljung-Box test
## 
## data:  Residuals from Linear regression model
## Q* = 31.552, df = 4, p-value = 2.361e-06
## 
## Model df: 2.   Total lags used: 6

shapiro.test(residuals(lmtPol))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(lmtPol)
## W = 0.90992, p-value = 0.01968

Uzyskane wyniki nie są w najmniejszym stopni żadnym zaskoczeniem. Jest oczywistym, że rezydua pochodzące z prostego modelu z trendem nie będą pochodziły z szumu białego i będą zawierały jeszcze sporo informacji. Co jest dość charakterystyczne to kształt skoków autokorelacji układający się w wyraźny sinusoidalny wzór charakterystyczny dla danych z silną sezonowością i zmienną wariancją. Oczywiście wysoka wartość statystyki $Q^*$, a co za tym idzie bardzo niskie p-value powoduje, że musimy odrzucić hipotezę zerową i stwierdzić, że rezydua z tego modelu ni pochodzą z szumu białego.

Model lm

checkresiduals(lmPol, test="LB")

## 
##  Ljung-Box test
## 
## data:  Residuals
## Q* = 9.7138, df = 3, p-value = 0.02116
## 
## Model df: 3.   Total lags used: 6

shapiro.test(residuals(lmPol))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(lmPol)
## W = 0.97446, p-value = 0.7035

W przypadku modelu lm rozkład rezyduów możemy uznać za pochodzący z rozkładu normalnego. Dodatkowo brak jest jakichkolwiek znaczących skoków autokorelacji, chociaż zachowały one sinusoidalny wzór zauważony wcześniej. To może być powodem dość wysokiej wartości statystyki $Q^*$ a co za tym idzie niskiej wartości p-value która na poziomie istotności $\alpha=0.05$ pozwala nam na odrzucenie hipotezy zerowej o pochodzeniu reszt z białego szumu.

Model tslm + season

checkresiduals(lmtsPol, test="LB")

## 
##  Ljung-Box test
## 
## data:  Residuals from Linear regression model
## Q* = 16.2, df = 3, p-value = 0.001032
## 
## Model df: 8.   Total lags used: 11

shapiro.test(residuals(lmtsPol))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(lmtsPol)
## W = 0.96615, p-value = 0.4821

Także w przypadku modelu tslm z sezonowością sinusoidalny wzór na wykresie ACF się utrzymuje. Widać też jeden znaczący skok dla opóźnienia 3. I chociaż nie mam podstaw do odrzucenia hipotezy o pochodzeniu rezyduów z rozkładu normalnego, to niestety zmuszony jestem odrzucić hipotezę o pochodzeniu reszt z szumu białego.

Model tslm + lambda

checkresiduals(lmtslPol, test="LB")

## 
##  Ljung-Box test
## 
## data:  Residuals from Linear regression model
## Q* = 7.1248, df = 3, p-value = 0.06803
## 
## Model df: 8.   Total lags used: 11

shapiro.test(residuals(lmtslPol))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(lmtslPol)
## W = 0.97029, p-value = 0.5884

Dopiero po zastosowaniu transformacji Boxa-Coxa udało się uzyskać model którego rezydua nie zawierają żadnej istotnej autokorelacji, a przy okazji wykres ACF utracił widoczny wcześniej wyraźny, sinusoidalny wzór. Nie jest wiec żadnym zaskoczeniem, że wartość statystki $Q^*$ testu Ljung-Boxa jest najmniejsza dla reszt z wszystkich badanych modeli. Możemy więc bez przeszkód, na poziomie istotności $\alpha=0.05$ pozostać przy hipotezie zerowej mówiącej o pochodzeniu reszt z szumu białego.

t) - w)

Obliczyć odchylenie standardowe składnika resztowego dla nowego modelu uwzględniającego wahania okresowe
Porównać odchylenia standardowe składników resztkowych dla modeli (starego i nowego)
Obliczyć współczynnik zmienności i determinacji dla nowego modelu
Czy nowy model jest lepiej dopasowany?

tibble(
  model = c("tslm", "lm", "tslm + season", "tslm + lambda"),
  s2.zt = c(round(sigma(lmtPol), 3),
         round(sigma(lmPol), 3),
         round(sigma(lmtsPol), 3),
         round(sigma(lmtslPol), 3)),
  phi2 = c(round(1-summary(lmtPol)$r.squared, 3),
           round(1-summary(lmPol)$r.squared, 3),
           round(1-summary(lmtsPol)$r.squared, 3),
           round(1-summary(lmtslPol)$r.squared, 3)),
  R2 = c(round(summary(lmtPol)$r.squared, 3),
           round(summary(lmPol)$r.squared, 3),
           round(summary(lmtsPol)$r.squared, 3),
           round(summary(lmtslPol)$r.squared, 3))
) %>% kable1()

model	s2.zt	phi2	R2
tslm	3357.150	0.474	0.526
lm	1311.930	0.070	0.930
tslm + season	1288.419	0.054	0.946
tslm + lambda	0.082	0.010	0.990

Analizując powyższą tablę, jeszcze raz mogłem się przekonać, że najlepiej dopasowany jest model uwzględniający sezonowość z zastosowaniem transformacji Baxa-Coxa. W przypadku tego modelu wartość współczynnika determinacji jest prawie równa 1. Jednak tak jak poprzednio należało pamiętać, że odchylenia standardowego reszt w przypadku modeli z transformacją nie będzie można porównywać z odchyleniem standardowym reszt z modeli bez transformacji.

Poniżej wyliczyłem to odchylenie w podobny sposób jak poprzednio.

s2.zt = 1/(n-8)*sum((tsPol-fitted(lmtslPol))^2)
cat("Odchylenie standardowe rezyduów wyliczone samodzielnie:", round(sqrt(s2.zt),3))

## Odchylenie standardowe rezyduów wyliczone samodzielnie: 702.663

Tym razem jest to najniższa wartość spośród wszystkich błędów standardowych reszt z pozostałych modeli.

x)

Ekstrapolować liczbę zakażeń na 13.03.2021 (prawdziwe wartości 13.03.2021 to 21063 dla Polski i 1115 dla warmińsko-mazurskiego oraz 17259 i 759 dla 14.03.2021)

Model tslm

dfRealPol = tibble(
  t = c(28,29)/7+1,
  y = c(21063, 17259))

plotPol(tsPol) + 
  autolayer(fitted(lmtPol), series="Trend liniowy")+
  autolayer(forecast(lmtPol, h=14), series="Prognoza na kolejne 2 tyg.")+
  geom_point(aes(t, y), dfRealPol, size=3, fill='red', shape=21)

Podobnie jak poprzednio tai i dla danych z całego krajumodel uwzględniający jedynie trend dał w rezultacie bardzo szerokie przedziały predykcji. Rzeczywista wartość występująca w dniu 13 marca różniła się już znacznie bardziej od wartości prognozowanej i znalazła się na końcu 80% przedziału predykcji. W przypadku wartości z kolejnego dnia prognoza była już bardzo zbliżona do rzeczywistej wartości.

Model lm

newPol = tibble(
  t = 1:(n+14),
  d =  rep(seasPol$gi, 6)) %>% 
  bind_cols(predict(lmPol, newdata = ., interval = "confidence")) %>% 
  tail(14)


plotPol(tsPol) +
  geom_line(aes(t/7+6/7, fit), data=newPol, color="blue")+
  geom_ribbon(aes(x=t/7+6/7, y=fit, ymin=lwr, ymax=upr), data=newPol, alpha =.2, fill="blue")+
  autolayer(ts(fitted(lmPol), frequency = 7), series="Model liniowy lm")+
  geom_point(aes(t, y), dfRealPol, size=3, fill='red', shape=21)

Niestety w przypadku klasycznego modelu liniowego z uwzględnieniem zmiennej w postaci wahań sezonowych prognoza zawiodła w obu przypadkach. Rzeczywista liczba zachorowań $(21063, 17259)$ znacznie odbiegała od prognozowanych wartości $(15798, 14853)$. Rzeczywiste wartości nie zmieściły się także w oszacowanych przedziłach predykcji które 13 marca wynosiły $[14742, 13693]$ a dnia następnego wynosiły $[13693, 16015]$. Niewątpliwie jest to spowodowane narastającą wariancją badanego szeregu czasowego, czego ten model w żadnej mierze nie uwzględniał.

model tslm + season

plotPol(tsPol) + 
  autolayer(fitted(lmtsPol), series="Model tslm + season")+
  autolayer(forecast(lmtsPol, h=14), series="Prognoza na kolejne 2 tyg.")+
  geom_point(aes(t, y), dfRealPol, size=3, fill='red', shape=21)

Równie słabo sprawdził się model pochodzący z funkcji tslm. Co prawda rzeczywista liczba zachorowań z dnia 14 marca zmieściła się w oszacowanym 80% przedziale predykcji ($[13630, 17712]$) jednak dla dnia poprzedniego wartość rzeczywista była poza 95% przedziałem predykcji ($[13961, 20385]$).

model tslm + lambda

plotPol(tsPol) + 
  autolayer(fitted(lmtslPol), series="Model tslm + lambda")+
  autolayer(forecast(lmtslPol, h=14), series="Prognoza na kolejne 2 tyg.")+
  geom_point(aes(t, y), dfRealPol, size=3, fill='red', shape=21)

Jedynym modelem który którego prognozy doskonale pokryły się z rzeczywistymi wartościami był model w którym zastosowano transformację Boxa-Coxa. W przypadku tego modelu prognozy na 13 i 14 marca były następujące $(19393, 16473)$ i nie różniły się zbytnio od wartości rzeczywistych które, przypomnijmy wynosiły $(21063, 17259)$. Warto także zwrócić uwagę na silnie zmienną szerokość przedziałów predykcji które najwęższe są dla niedziel, poniedziałków i wtorków a znacznie szersze dla pozostałych dni tygodnia. Taki efekt występuje jedynie dla tego modelu i jest niewątpliwie związany z budową samego modelu (osobne współczynniki regresji dla każdego sezonu) jak i oczywiście z zastosowaną transformacją Boxa-Coxa.

Wnioski

Przeprowadzone zadanie pozwoliło mi na zapoznanie się z niektórymi z metod prognozowania z wykorzystaniem modeli regresji zastosowanych dla szeregów z trendem oraz sezonowością, wraz z całą metodologią z tym związaną w połączeniu odpowiednią diagnostyką rezyduów. Największym zaskoczeniem był dla mnie bardzo znaczący efekt zastosowania transformacji Boxa-Coxa. Jak mogłem się “naocznie” przekonać, ma ona wpływ nie tylko na szerokość przedziałów predykcji ale także bardzo istotnie wpływa na samą jakość predykcyjną modeli. Nie bez znaczenia są także wnioski płynące z diagnostyki rezyduów. W przypadku modelowania szeregów czasowych diagnostyka ta jest nieco odmienna od diagnostyki stosowanej do klasycznych modeli regresji. Szczególne znaczenie ma w tym przypadku analiza autokorelacji wraz z odpowiednim testem portmanteau.

Box, G. E. P., & Cox, D. R. (1964). An analysis of transformations. Journal of the Royal Statistical Society. Series B, Statistical Methodology, 26(2), 211–252.↩︎

ADPM - Zadanie korona wirus

Marek Fiołka

25 02 2022

Przygotowania wstępne

Biblioteki

Wczytanie danych

Warmińsko-Mazurskie

a)

b) - c)

d)

e)

f)

g)

h)

i)

j)

k)

l) - o)

n)

p)

q)

Model lm

Model tslm + season

Model tslm + lambda

r)

s)

t)

Model tslm

Model lm

Model tslm + season

Model tslm + lambda

t) - w)

x)

Model tslm

Model lm

model tslm + season

model tslm + lambda

Polska

a)

b) - c)

d)

e)

f)

g)

h)

i)

j)

k)

l) - o)

n)

p)

q)

Model lm

Model tslm + season

Model tslm + lambda

r)

s)

t)

Model tslm

Model lm

Model tslm + season

Model tslm + lambda

t) - w)

x)

Model tslm

Model lm

model tslm + season

model tslm + lambda

Wnioski